• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 //===-- APFloat.cpp - Implement APFloat class -----------------------------===//
2 //
3 //                     The LLVM Compiler Infrastructure
4 //
5 // This file is distributed under the University of Illinois Open Source
6 // License. See LICENSE.TXT for details.
7 //
8 //===----------------------------------------------------------------------===//
9 //
10 // This file implements a class to represent arbitrary precision floating
11 // point values and provide a variety of arithmetic operations on them.
12 //
13 //===----------------------------------------------------------------------===//
14 
15 #include "llvm/ADT/APFloat.h"
16 #include "llvm/ADT/APSInt.h"
17 #include "llvm/ADT/StringRef.h"
18 #include "llvm/ADT/FoldingSet.h"
19 #include "llvm/Support/ErrorHandling.h"
20 #include "llvm/Support/MathExtras.h"
21 #include <limits.h>
22 #include <cstring>
23 
24 using namespace llvm;
25 
26 #define convolve(lhs, rhs) ((lhs) * 4 + (rhs))
27 
28 /* Assumed in hexadecimal significand parsing, and conversion to
29    hexadecimal strings.  */
30 #define COMPILE_TIME_ASSERT(cond) extern int CTAssert[(cond) ? 1 : -1]
31 COMPILE_TIME_ASSERT(integerPartWidth % 4 == 0);
32 
33 namespace llvm {
34 
35   /* Represents floating point arithmetic semantics.  */
36   struct fltSemantics {
37     /* The largest E such that 2^E is representable; this matches the
38        definition of IEEE 754.  */
39     exponent_t maxExponent;
40 
41     /* The smallest E such that 2^E is a normalized number; this
42        matches the definition of IEEE 754.  */
43     exponent_t minExponent;
44 
45     /* Number of bits in the significand.  This includes the integer
46        bit.  */
47     unsigned int precision;
48 
49     /* True if arithmetic is supported.  */
50     unsigned int arithmeticOK;
51   };
52 
53   const fltSemantics APFloat::IEEEhalf = { 15, -14, 11, true };
54   const fltSemantics APFloat::IEEEsingle = { 127, -126, 24, true };
55   const fltSemantics APFloat::IEEEdouble = { 1023, -1022, 53, true };
56   const fltSemantics APFloat::IEEEquad = { 16383, -16382, 113, true };
57   const fltSemantics APFloat::x87DoubleExtended = { 16383, -16382, 64, true };
58   const fltSemantics APFloat::Bogus = { 0, 0, 0, true };
59 
60   // The PowerPC format consists of two doubles.  It does not map cleanly
61   // onto the usual format above.  For now only storage of constants of
62   // this type is supported, no arithmetic.
63   const fltSemantics APFloat::PPCDoubleDouble = { 1023, -1022, 106, false };
64 
65   /* A tight upper bound on number of parts required to hold the value
66      pow(5, power) is
67 
68        power * 815 / (351 * integerPartWidth) + 1
69 
70      However, whilst the result may require only this many parts,
71      because we are multiplying two values to get it, the
72      multiplication may require an extra part with the excess part
73      being zero (consider the trivial case of 1 * 1, tcFullMultiply
74      requires two parts to hold the single-part result).  So we add an
75      extra one to guarantee enough space whilst multiplying.  */
76   const unsigned int maxExponent = 16383;
77   const unsigned int maxPrecision = 113;
78   const unsigned int maxPowerOfFiveExponent = maxExponent + maxPrecision - 1;
79   const unsigned int maxPowerOfFiveParts = 2 + ((maxPowerOfFiveExponent * 815)
80                                                 / (351 * integerPartWidth));
81 }
82 
83 /* A bunch of private, handy routines.  */
84 
85 static inline unsigned int
partCountForBits(unsigned int bits)86 partCountForBits(unsigned int bits)
87 {
88   return ((bits) + integerPartWidth - 1) / integerPartWidth;
89 }
90 
91 /* Returns 0U-9U.  Return values >= 10U are not digits.  */
92 static inline unsigned int
decDigitValue(unsigned int c)93 decDigitValue(unsigned int c)
94 {
95   return c - '0';
96 }
97 
98 static unsigned int
hexDigitValue(unsigned int c)99 hexDigitValue(unsigned int c)
100 {
101   unsigned int r;
102 
103   r = c - '0';
104   if (r <= 9)
105     return r;
106 
107   r = c - 'A';
108   if (r <= 5)
109     return r + 10;
110 
111   r = c - 'a';
112   if (r <= 5)
113     return r + 10;
114 
115   return -1U;
116 }
117 
118 static inline void
assertArithmeticOK(const llvm::fltSemantics & semantics)119 assertArithmeticOK(const llvm::fltSemantics &semantics) {
120   assert(semantics.arithmeticOK &&
121          "Compile-time arithmetic does not support these semantics");
122 }
123 
124 /* Return the value of a decimal exponent of the form
125    [+-]ddddddd.
126 
127    If the exponent overflows, returns a large exponent with the
128    appropriate sign.  */
129 static int
readExponent(StringRef::iterator begin,StringRef::iterator end)130 readExponent(StringRef::iterator begin, StringRef::iterator end)
131 {
132   bool isNegative;
133   unsigned int absExponent;
134   const unsigned int overlargeExponent = 24000;  /* FIXME.  */
135   StringRef::iterator p = begin;
136 
137   assert(p != end && "Exponent has no digits");
138 
139   isNegative = (*p == '-');
140   if (*p == '-' || *p == '+') {
141     p++;
142     assert(p != end && "Exponent has no digits");
143   }
144 
145   absExponent = decDigitValue(*p++);
146   assert(absExponent < 10U && "Invalid character in exponent");
147 
148   for (; p != end; ++p) {
149     unsigned int value;
150 
151     value = decDigitValue(*p);
152     assert(value < 10U && "Invalid character in exponent");
153 
154     value += absExponent * 10;
155     if (absExponent >= overlargeExponent) {
156       absExponent = overlargeExponent;
157       p = end;  /* outwit assert below */
158       break;
159     }
160     absExponent = value;
161   }
162 
163   assert(p == end && "Invalid exponent in exponent");
164 
165   if (isNegative)
166     return -(int) absExponent;
167   else
168     return (int) absExponent;
169 }
170 
171 /* This is ugly and needs cleaning up, but I don't immediately see
172    how whilst remaining safe.  */
173 static int
totalExponent(StringRef::iterator p,StringRef::iterator end,int exponentAdjustment)174 totalExponent(StringRef::iterator p, StringRef::iterator end,
175               int exponentAdjustment)
176 {
177   int unsignedExponent;
178   bool negative, overflow;
179   int exponent = 0;
180 
181   assert(p != end && "Exponent has no digits");
182 
183   negative = *p == '-';
184   if (*p == '-' || *p == '+') {
185     p++;
186     assert(p != end && "Exponent has no digits");
187   }
188 
189   unsignedExponent = 0;
190   overflow = false;
191   for (; p != end; ++p) {
192     unsigned int value;
193 
194     value = decDigitValue(*p);
195     assert(value < 10U && "Invalid character in exponent");
196 
197     unsignedExponent = unsignedExponent * 10 + value;
198     if (unsignedExponent > 32767)
199       overflow = true;
200   }
201 
202   if (exponentAdjustment > 32767 || exponentAdjustment < -32768)
203     overflow = true;
204 
205   if (!overflow) {
206     exponent = unsignedExponent;
207     if (negative)
208       exponent = -exponent;
209     exponent += exponentAdjustment;
210     if (exponent > 32767 || exponent < -32768)
211       overflow = true;
212   }
213 
214   if (overflow)
215     exponent = negative ? -32768: 32767;
216 
217   return exponent;
218 }
219 
220 static StringRef::iterator
skipLeadingZeroesAndAnyDot(StringRef::iterator begin,StringRef::iterator end,StringRef::iterator * dot)221 skipLeadingZeroesAndAnyDot(StringRef::iterator begin, StringRef::iterator end,
222                            StringRef::iterator *dot)
223 {
224   StringRef::iterator p = begin;
225   *dot = end;
226   while (*p == '0' && p != end)
227     p++;
228 
229   if (*p == '.') {
230     *dot = p++;
231 
232     assert(end - begin != 1 && "Significand has no digits");
233 
234     while (*p == '0' && p != end)
235       p++;
236   }
237 
238   return p;
239 }
240 
241 /* Given a normal decimal floating point number of the form
242 
243      dddd.dddd[eE][+-]ddd
244 
245    where the decimal point and exponent are optional, fill out the
246    structure D.  Exponent is appropriate if the significand is
247    treated as an integer, and normalizedExponent if the significand
248    is taken to have the decimal point after a single leading
249    non-zero digit.
250 
251    If the value is zero, V->firstSigDigit points to a non-digit, and
252    the return exponent is zero.
253 */
254 struct decimalInfo {
255   const char *firstSigDigit;
256   const char *lastSigDigit;
257   int exponent;
258   int normalizedExponent;
259 };
260 
261 static void
interpretDecimal(StringRef::iterator begin,StringRef::iterator end,decimalInfo * D)262 interpretDecimal(StringRef::iterator begin, StringRef::iterator end,
263                  decimalInfo *D)
264 {
265   StringRef::iterator dot = end;
266   StringRef::iterator p = skipLeadingZeroesAndAnyDot (begin, end, &dot);
267 
268   D->firstSigDigit = p;
269   D->exponent = 0;
270   D->normalizedExponent = 0;
271 
272   for (; p != end; ++p) {
273     if (*p == '.') {
274       assert(dot == end && "String contains multiple dots");
275       dot = p++;
276       if (p == end)
277         break;
278     }
279     if (decDigitValue(*p) >= 10U)
280       break;
281   }
282 
283   if (p != end) {
284     assert((*p == 'e' || *p == 'E') && "Invalid character in significand");
285     assert(p != begin && "Significand has no digits");
286     assert((dot == end || p - begin != 1) && "Significand has no digits");
287 
288     /* p points to the first non-digit in the string */
289     D->exponent = readExponent(p + 1, end);
290 
291     /* Implied decimal point?  */
292     if (dot == end)
293       dot = p;
294   }
295 
296   /* If number is all zeroes accept any exponent.  */
297   if (p != D->firstSigDigit) {
298     /* Drop insignificant trailing zeroes.  */
299     if (p != begin) {
300       do
301         do
302           p--;
303         while (p != begin && *p == '0');
304       while (p != begin && *p == '.');
305     }
306 
307     /* Adjust the exponents for any decimal point.  */
308     D->exponent += static_cast<exponent_t>((dot - p) - (dot > p));
309     D->normalizedExponent = (D->exponent +
310               static_cast<exponent_t>((p - D->firstSigDigit)
311                                       - (dot > D->firstSigDigit && dot < p)));
312   }
313 
314   D->lastSigDigit = p;
315 }
316 
317 /* Return the trailing fraction of a hexadecimal number.
318    DIGITVALUE is the first hex digit of the fraction, P points to
319    the next digit.  */
320 static lostFraction
trailingHexadecimalFraction(StringRef::iterator p,StringRef::iterator end,unsigned int digitValue)321 trailingHexadecimalFraction(StringRef::iterator p, StringRef::iterator end,
322                             unsigned int digitValue)
323 {
324   unsigned int hexDigit;
325 
326   /* If the first trailing digit isn't 0 or 8 we can work out the
327      fraction immediately.  */
328   if (digitValue > 8)
329     return lfMoreThanHalf;
330   else if (digitValue < 8 && digitValue > 0)
331     return lfLessThanHalf;
332 
333   /* Otherwise we need to find the first non-zero digit.  */
334   while (*p == '0')
335     p++;
336 
337   assert(p != end && "Invalid trailing hexadecimal fraction!");
338 
339   hexDigit = hexDigitValue(*p);
340 
341   /* If we ran off the end it is exactly zero or one-half, otherwise
342      a little more.  */
343   if (hexDigit == -1U)
344     return digitValue == 0 ? lfExactlyZero: lfExactlyHalf;
345   else
346     return digitValue == 0 ? lfLessThanHalf: lfMoreThanHalf;
347 }
348 
349 /* Return the fraction lost were a bignum truncated losing the least
350    significant BITS bits.  */
351 static lostFraction
lostFractionThroughTruncation(const integerPart * parts,unsigned int partCount,unsigned int bits)352 lostFractionThroughTruncation(const integerPart *parts,
353                               unsigned int partCount,
354                               unsigned int bits)
355 {
356   unsigned int lsb;
357 
358   lsb = APInt::tcLSB(parts, partCount);
359 
360   /* Note this is guaranteed true if bits == 0, or LSB == -1U.  */
361   if (bits <= lsb)
362     return lfExactlyZero;
363   if (bits == lsb + 1)
364     return lfExactlyHalf;
365   if (bits <= partCount * integerPartWidth &&
366       APInt::tcExtractBit(parts, bits - 1))
367     return lfMoreThanHalf;
368 
369   return lfLessThanHalf;
370 }
371 
372 /* Shift DST right BITS bits noting lost fraction.  */
373 static lostFraction
shiftRight(integerPart * dst,unsigned int parts,unsigned int bits)374 shiftRight(integerPart *dst, unsigned int parts, unsigned int bits)
375 {
376   lostFraction lost_fraction;
377 
378   lost_fraction = lostFractionThroughTruncation(dst, parts, bits);
379 
380   APInt::tcShiftRight(dst, parts, bits);
381 
382   return lost_fraction;
383 }
384 
385 /* Combine the effect of two lost fractions.  */
386 static lostFraction
combineLostFractions(lostFraction moreSignificant,lostFraction lessSignificant)387 combineLostFractions(lostFraction moreSignificant,
388                      lostFraction lessSignificant)
389 {
390   if (lessSignificant != lfExactlyZero) {
391     if (moreSignificant == lfExactlyZero)
392       moreSignificant = lfLessThanHalf;
393     else if (moreSignificant == lfExactlyHalf)
394       moreSignificant = lfMoreThanHalf;
395   }
396 
397   return moreSignificant;
398 }
399 
400 /* The error from the true value, in half-ulps, on multiplying two
401    floating point numbers, which differ from the value they
402    approximate by at most HUE1 and HUE2 half-ulps, is strictly less
403    than the returned value.
404 
405    See "How to Read Floating Point Numbers Accurately" by William D
406    Clinger.  */
407 static unsigned int
HUerrBound(bool inexactMultiply,unsigned int HUerr1,unsigned int HUerr2)408 HUerrBound(bool inexactMultiply, unsigned int HUerr1, unsigned int HUerr2)
409 {
410   assert(HUerr1 < 2 || HUerr2 < 2 || (HUerr1 + HUerr2 < 8));
411 
412   if (HUerr1 + HUerr2 == 0)
413     return inexactMultiply * 2;  /* <= inexactMultiply half-ulps.  */
414   else
415     return inexactMultiply + 2 * (HUerr1 + HUerr2);
416 }
417 
418 /* The number of ulps from the boundary (zero, or half if ISNEAREST)
419    when the least significant BITS are truncated.  BITS cannot be
420    zero.  */
421 static integerPart
ulpsFromBoundary(const integerPart * parts,unsigned int bits,bool isNearest)422 ulpsFromBoundary(const integerPart *parts, unsigned int bits, bool isNearest)
423 {
424   unsigned int count, partBits;
425   integerPart part, boundary;
426 
427   assert(bits != 0);
428 
429   bits--;
430   count = bits / integerPartWidth;
431   partBits = bits % integerPartWidth + 1;
432 
433   part = parts[count] & (~(integerPart) 0 >> (integerPartWidth - partBits));
434 
435   if (isNearest)
436     boundary = (integerPart) 1 << (partBits - 1);
437   else
438     boundary = 0;
439 
440   if (count == 0) {
441     if (part - boundary <= boundary - part)
442       return part - boundary;
443     else
444       return boundary - part;
445   }
446 
447   if (part == boundary) {
448     while (--count)
449       if (parts[count])
450         return ~(integerPart) 0; /* A lot.  */
451 
452     return parts[0];
453   } else if (part == boundary - 1) {
454     while (--count)
455       if (~parts[count])
456         return ~(integerPart) 0; /* A lot.  */
457 
458     return -parts[0];
459   }
460 
461   return ~(integerPart) 0; /* A lot.  */
462 }
463 
464 /* Place pow(5, power) in DST, and return the number of parts used.
465    DST must be at least one part larger than size of the answer.  */
466 static unsigned int
powerOf5(integerPart * dst,unsigned int power)467 powerOf5(integerPart *dst, unsigned int power)
468 {
469   static const integerPart firstEightPowers[] = { 1, 5, 25, 125, 625, 3125,
470                                                   15625, 78125 };
471   integerPart pow5s[maxPowerOfFiveParts * 2 + 5];
472   pow5s[0] = 78125 * 5;
473 
474   unsigned int partsCount[16] = { 1 };
475   integerPart scratch[maxPowerOfFiveParts], *p1, *p2, *pow5;
476   unsigned int result;
477   assert(power <= maxExponent);
478 
479   p1 = dst;
480   p2 = scratch;
481 
482   *p1 = firstEightPowers[power & 7];
483   power >>= 3;
484 
485   result = 1;
486   pow5 = pow5s;
487 
488   for (unsigned int n = 0; power; power >>= 1, n++) {
489     unsigned int pc;
490 
491     pc = partsCount[n];
492 
493     /* Calculate pow(5,pow(2,n+3)) if we haven't yet.  */
494     if (pc == 0) {
495       pc = partsCount[n - 1];
496       APInt::tcFullMultiply(pow5, pow5 - pc, pow5 - pc, pc, pc);
497       pc *= 2;
498       if (pow5[pc - 1] == 0)
499         pc--;
500       partsCount[n] = pc;
501     }
502 
503     if (power & 1) {
504       integerPart *tmp;
505 
506       APInt::tcFullMultiply(p2, p1, pow5, result, pc);
507       result += pc;
508       if (p2[result - 1] == 0)
509         result--;
510 
511       /* Now result is in p1 with partsCount parts and p2 is scratch
512          space.  */
513       tmp = p1, p1 = p2, p2 = tmp;
514     }
515 
516     pow5 += pc;
517   }
518 
519   if (p1 != dst)
520     APInt::tcAssign(dst, p1, result);
521 
522   return result;
523 }
524 
525 /* Zero at the end to avoid modular arithmetic when adding one; used
526    when rounding up during hexadecimal output.  */
527 static const char hexDigitsLower[] = "0123456789abcdef0";
528 static const char hexDigitsUpper[] = "0123456789ABCDEF0";
529 static const char infinityL[] = "infinity";
530 static const char infinityU[] = "INFINITY";
531 static const char NaNL[] = "nan";
532 static const char NaNU[] = "NAN";
533 
534 /* Write out an integerPart in hexadecimal, starting with the most
535    significant nibble.  Write out exactly COUNT hexdigits, return
536    COUNT.  */
537 static unsigned int
partAsHex(char * dst,integerPart part,unsigned int count,const char * hexDigitChars)538 partAsHex (char *dst, integerPart part, unsigned int count,
539            const char *hexDigitChars)
540 {
541   unsigned int result = count;
542 
543   assert(count != 0 && count <= integerPartWidth / 4);
544 
545   part >>= (integerPartWidth - 4 * count);
546   while (count--) {
547     dst[count] = hexDigitChars[part & 0xf];
548     part >>= 4;
549   }
550 
551   return result;
552 }
553 
554 /* Write out an unsigned decimal integer.  */
555 static char *
writeUnsignedDecimal(char * dst,unsigned int n)556 writeUnsignedDecimal (char *dst, unsigned int n)
557 {
558   char buff[40], *p;
559 
560   p = buff;
561   do
562     *p++ = '0' + n % 10;
563   while (n /= 10);
564 
565   do
566     *dst++ = *--p;
567   while (p != buff);
568 
569   return dst;
570 }
571 
572 /* Write out a signed decimal integer.  */
573 static char *
writeSignedDecimal(char * dst,int value)574 writeSignedDecimal (char *dst, int value)
575 {
576   if (value < 0) {
577     *dst++ = '-';
578     dst = writeUnsignedDecimal(dst, -(unsigned) value);
579   } else
580     dst = writeUnsignedDecimal(dst, value);
581 
582   return dst;
583 }
584 
585 /* Constructors.  */
586 void
initialize(const fltSemantics * ourSemantics)587 APFloat::initialize(const fltSemantics *ourSemantics)
588 {
589   unsigned int count;
590 
591   semantics = ourSemantics;
592   count = partCount();
593   if (count > 1)
594     significand.parts = new integerPart[count];
595 }
596 
597 void
freeSignificand()598 APFloat::freeSignificand()
599 {
600   if (partCount() > 1)
601     delete [] significand.parts;
602 }
603 
604 void
assign(const APFloat & rhs)605 APFloat::assign(const APFloat &rhs)
606 {
607   assert(semantics == rhs.semantics);
608 
609   sign = rhs.sign;
610   category = rhs.category;
611   exponent = rhs.exponent;
612   sign2 = rhs.sign2;
613   exponent2 = rhs.exponent2;
614   if (category == fcNormal || category == fcNaN)
615     copySignificand(rhs);
616 }
617 
618 void
copySignificand(const APFloat & rhs)619 APFloat::copySignificand(const APFloat &rhs)
620 {
621   assert(category == fcNormal || category == fcNaN);
622   assert(rhs.partCount() >= partCount());
623 
624   APInt::tcAssign(significandParts(), rhs.significandParts(),
625                   partCount());
626 }
627 
628 /* Make this number a NaN, with an arbitrary but deterministic value
629    for the significand.  If double or longer, this is a signalling NaN,
630    which may not be ideal.  If float, this is QNaN(0).  */
makeNaN(bool SNaN,bool Negative,const APInt * fill)631 void APFloat::makeNaN(bool SNaN, bool Negative, const APInt *fill)
632 {
633   category = fcNaN;
634   sign = Negative;
635 
636   integerPart *significand = significandParts();
637   unsigned numParts = partCount();
638 
639   // Set the significand bits to the fill.
640   if (!fill || fill->getNumWords() < numParts)
641     APInt::tcSet(significand, 0, numParts);
642   if (fill) {
643     APInt::tcAssign(significand, fill->getRawData(),
644                     std::min(fill->getNumWords(), numParts));
645 
646     // Zero out the excess bits of the significand.
647     unsigned bitsToPreserve = semantics->precision - 1;
648     unsigned part = bitsToPreserve / 64;
649     bitsToPreserve %= 64;
650     significand[part] &= ((1ULL << bitsToPreserve) - 1);
651     for (part++; part != numParts; ++part)
652       significand[part] = 0;
653   }
654 
655   unsigned QNaNBit = semantics->precision - 2;
656 
657   if (SNaN) {
658     // We always have to clear the QNaN bit to make it an SNaN.
659     APInt::tcClearBit(significand, QNaNBit);
660 
661     // If there are no bits set in the payload, we have to set
662     // *something* to make it a NaN instead of an infinity;
663     // conventionally, this is the next bit down from the QNaN bit.
664     if (APInt::tcIsZero(significand, numParts))
665       APInt::tcSetBit(significand, QNaNBit - 1);
666   } else {
667     // We always have to set the QNaN bit to make it a QNaN.
668     APInt::tcSetBit(significand, QNaNBit);
669   }
670 
671   // For x87 extended precision, we want to make a NaN, not a
672   // pseudo-NaN.  Maybe we should expose the ability to make
673   // pseudo-NaNs?
674   if (semantics == &APFloat::x87DoubleExtended)
675     APInt::tcSetBit(significand, QNaNBit + 1);
676 }
677 
makeNaN(const fltSemantics & Sem,bool SNaN,bool Negative,const APInt * fill)678 APFloat APFloat::makeNaN(const fltSemantics &Sem, bool SNaN, bool Negative,
679                          const APInt *fill) {
680   APFloat value(Sem, uninitialized);
681   value.makeNaN(SNaN, Negative, fill);
682   return value;
683 }
684 
685 APFloat &
operator =(const APFloat & rhs)686 APFloat::operator=(const APFloat &rhs)
687 {
688   if (this != &rhs) {
689     if (semantics != rhs.semantics) {
690       freeSignificand();
691       initialize(rhs.semantics);
692     }
693     assign(rhs);
694   }
695 
696   return *this;
697 }
698 
699 bool
bitwiseIsEqual(const APFloat & rhs) const700 APFloat::bitwiseIsEqual(const APFloat &rhs) const {
701   if (this == &rhs)
702     return true;
703   if (semantics != rhs.semantics ||
704       category != rhs.category ||
705       sign != rhs.sign)
706     return false;
707   if (semantics==(const llvm::fltSemantics*)&PPCDoubleDouble &&
708       sign2 != rhs.sign2)
709     return false;
710   if (category==fcZero || category==fcInfinity)
711     return true;
712   else if (category==fcNormal && exponent!=rhs.exponent)
713     return false;
714   else if (semantics==(const llvm::fltSemantics*)&PPCDoubleDouble &&
715            exponent2!=rhs.exponent2)
716     return false;
717   else {
718     int i= partCount();
719     const integerPart* p=significandParts();
720     const integerPart* q=rhs.significandParts();
721     for (; i>0; i--, p++, q++) {
722       if (*p != *q)
723         return false;
724     }
725     return true;
726   }
727 }
728 
APFloat(const fltSemantics & ourSemantics,integerPart value)729 APFloat::APFloat(const fltSemantics &ourSemantics, integerPart value)
730   : exponent2(0), sign2(0) {
731   assertArithmeticOK(ourSemantics);
732   initialize(&ourSemantics);
733   sign = 0;
734   zeroSignificand();
735   exponent = ourSemantics.precision - 1;
736   significandParts()[0] = value;
737   normalize(rmNearestTiesToEven, lfExactlyZero);
738 }
739 
APFloat(const fltSemantics & ourSemantics)740 APFloat::APFloat(const fltSemantics &ourSemantics) : exponent2(0), sign2(0) {
741   assertArithmeticOK(ourSemantics);
742   initialize(&ourSemantics);
743   category = fcZero;
744   sign = false;
745 }
746 
APFloat(const fltSemantics & ourSemantics,uninitializedTag tag)747 APFloat::APFloat(const fltSemantics &ourSemantics, uninitializedTag tag)
748   : exponent2(0), sign2(0) {
749   assertArithmeticOK(ourSemantics);
750   // Allocates storage if necessary but does not initialize it.
751   initialize(&ourSemantics);
752 }
753 
APFloat(const fltSemantics & ourSemantics,fltCategory ourCategory,bool negative)754 APFloat::APFloat(const fltSemantics &ourSemantics,
755                  fltCategory ourCategory, bool negative)
756   : exponent2(0), sign2(0) {
757   assertArithmeticOK(ourSemantics);
758   initialize(&ourSemantics);
759   category = ourCategory;
760   sign = negative;
761   if (category == fcNormal)
762     category = fcZero;
763   else if (ourCategory == fcNaN)
764     makeNaN();
765 }
766 
APFloat(const fltSemantics & ourSemantics,StringRef text)767 APFloat::APFloat(const fltSemantics &ourSemantics, StringRef text)
768   : exponent2(0), sign2(0) {
769   assertArithmeticOK(ourSemantics);
770   initialize(&ourSemantics);
771   convertFromString(text, rmNearestTiesToEven);
772 }
773 
APFloat(const APFloat & rhs)774 APFloat::APFloat(const APFloat &rhs) : exponent2(0), sign2(0) {
775   initialize(rhs.semantics);
776   assign(rhs);
777 }
778 
~APFloat()779 APFloat::~APFloat()
780 {
781   freeSignificand();
782 }
783 
784 // Profile - This method 'profiles' an APFloat for use with FoldingSet.
Profile(FoldingSetNodeID & ID) const785 void APFloat::Profile(FoldingSetNodeID& ID) const {
786   ID.Add(bitcastToAPInt());
787 }
788 
789 unsigned int
partCount() const790 APFloat::partCount() const
791 {
792   return partCountForBits(semantics->precision + 1);
793 }
794 
795 unsigned int
semanticsPrecision(const fltSemantics & semantics)796 APFloat::semanticsPrecision(const fltSemantics &semantics)
797 {
798   return semantics.precision;
799 }
800 
801 const integerPart *
significandParts() const802 APFloat::significandParts() const
803 {
804   return const_cast<APFloat *>(this)->significandParts();
805 }
806 
807 integerPart *
significandParts()808 APFloat::significandParts()
809 {
810   assert(category == fcNormal || category == fcNaN);
811 
812   if (partCount() > 1)
813     return significand.parts;
814   else
815     return &significand.part;
816 }
817 
818 void
zeroSignificand()819 APFloat::zeroSignificand()
820 {
821   category = fcNormal;
822   APInt::tcSet(significandParts(), 0, partCount());
823 }
824 
825 /* Increment an fcNormal floating point number's significand.  */
826 void
incrementSignificand()827 APFloat::incrementSignificand()
828 {
829   integerPart carry;
830 
831   carry = APInt::tcIncrement(significandParts(), partCount());
832 
833   /* Our callers should never cause us to overflow.  */
834   assert(carry == 0);
835 }
836 
837 /* Add the significand of the RHS.  Returns the carry flag.  */
838 integerPart
addSignificand(const APFloat & rhs)839 APFloat::addSignificand(const APFloat &rhs)
840 {
841   integerPart *parts;
842 
843   parts = significandParts();
844 
845   assert(semantics == rhs.semantics);
846   assert(exponent == rhs.exponent);
847 
848   return APInt::tcAdd(parts, rhs.significandParts(), 0, partCount());
849 }
850 
851 /* Subtract the significand of the RHS with a borrow flag.  Returns
852    the borrow flag.  */
853 integerPart
subtractSignificand(const APFloat & rhs,integerPart borrow)854 APFloat::subtractSignificand(const APFloat &rhs, integerPart borrow)
855 {
856   integerPart *parts;
857 
858   parts = significandParts();
859 
860   assert(semantics == rhs.semantics);
861   assert(exponent == rhs.exponent);
862 
863   return APInt::tcSubtract(parts, rhs.significandParts(), borrow,
864                            partCount());
865 }
866 
867 /* Multiply the significand of the RHS.  If ADDEND is non-NULL, add it
868    on to the full-precision result of the multiplication.  Returns the
869    lost fraction.  */
870 lostFraction
multiplySignificand(const APFloat & rhs,const APFloat * addend)871 APFloat::multiplySignificand(const APFloat &rhs, const APFloat *addend)
872 {
873   unsigned int omsb;        // One, not zero, based MSB.
874   unsigned int partsCount, newPartsCount, precision;
875   integerPart *lhsSignificand;
876   integerPart scratch[4];
877   integerPart *fullSignificand;
878   lostFraction lost_fraction;
879   bool ignored;
880 
881   assert(semantics == rhs.semantics);
882 
883   precision = semantics->precision;
884   newPartsCount = partCountForBits(precision * 2);
885 
886   if (newPartsCount > 4)
887     fullSignificand = new integerPart[newPartsCount];
888   else
889     fullSignificand = scratch;
890 
891   lhsSignificand = significandParts();
892   partsCount = partCount();
893 
894   APInt::tcFullMultiply(fullSignificand, lhsSignificand,
895                         rhs.significandParts(), partsCount, partsCount);
896 
897   lost_fraction = lfExactlyZero;
898   omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
899   exponent += rhs.exponent;
900 
901   if (addend) {
902     Significand savedSignificand = significand;
903     const fltSemantics *savedSemantics = semantics;
904     fltSemantics extendedSemantics;
905     opStatus status;
906     unsigned int extendedPrecision;
907 
908     /* Normalize our MSB.  */
909     extendedPrecision = precision + precision - 1;
910     if (omsb != extendedPrecision) {
911       APInt::tcShiftLeft(fullSignificand, newPartsCount,
912                          extendedPrecision - omsb);
913       exponent -= extendedPrecision - omsb;
914     }
915 
916     /* Create new semantics.  */
917     extendedSemantics = *semantics;
918     extendedSemantics.precision = extendedPrecision;
919 
920     if (newPartsCount == 1)
921       significand.part = fullSignificand[0];
922     else
923       significand.parts = fullSignificand;
924     semantics = &extendedSemantics;
925 
926     APFloat extendedAddend(*addend);
927     status = extendedAddend.convert(extendedSemantics, rmTowardZero, &ignored);
928     assert(status == opOK);
929     lost_fraction = addOrSubtractSignificand(extendedAddend, false);
930 
931     /* Restore our state.  */
932     if (newPartsCount == 1)
933       fullSignificand[0] = significand.part;
934     significand = savedSignificand;
935     semantics = savedSemantics;
936 
937     omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
938   }
939 
940   exponent -= (precision - 1);
941 
942   if (omsb > precision) {
943     unsigned int bits, significantParts;
944     lostFraction lf;
945 
946     bits = omsb - precision;
947     significantParts = partCountForBits(omsb);
948     lf = shiftRight(fullSignificand, significantParts, bits);
949     lost_fraction = combineLostFractions(lf, lost_fraction);
950     exponent += bits;
951   }
952 
953   APInt::tcAssign(lhsSignificand, fullSignificand, partsCount);
954 
955   if (newPartsCount > 4)
956     delete [] fullSignificand;
957 
958   return lost_fraction;
959 }
960 
961 /* Multiply the significands of LHS and RHS to DST.  */
962 lostFraction
divideSignificand(const APFloat & rhs)963 APFloat::divideSignificand(const APFloat &rhs)
964 {
965   unsigned int bit, i, partsCount;
966   const integerPart *rhsSignificand;
967   integerPart *lhsSignificand, *dividend, *divisor;
968   integerPart scratch[4];
969   lostFraction lost_fraction;
970 
971   assert(semantics == rhs.semantics);
972 
973   lhsSignificand = significandParts();
974   rhsSignificand = rhs.significandParts();
975   partsCount = partCount();
976 
977   if (partsCount > 2)
978     dividend = new integerPart[partsCount * 2];
979   else
980     dividend = scratch;
981 
982   divisor = dividend + partsCount;
983 
984   /* Copy the dividend and divisor as they will be modified in-place.  */
985   for (i = 0; i < partsCount; i++) {
986     dividend[i] = lhsSignificand[i];
987     divisor[i] = rhsSignificand[i];
988     lhsSignificand[i] = 0;
989   }
990 
991   exponent -= rhs.exponent;
992 
993   unsigned int precision = semantics->precision;
994 
995   /* Normalize the divisor.  */
996   bit = precision - APInt::tcMSB(divisor, partsCount) - 1;
997   if (bit) {
998     exponent += bit;
999     APInt::tcShiftLeft(divisor, partsCount, bit);
1000   }
1001 
1002   /* Normalize the dividend.  */
1003   bit = precision - APInt::tcMSB(dividend, partsCount) - 1;
1004   if (bit) {
1005     exponent -= bit;
1006     APInt::tcShiftLeft(dividend, partsCount, bit);
1007   }
1008 
1009   /* Ensure the dividend >= divisor initially for the loop below.
1010      Incidentally, this means that the division loop below is
1011      guaranteed to set the integer bit to one.  */
1012   if (APInt::tcCompare(dividend, divisor, partsCount) < 0) {
1013     exponent--;
1014     APInt::tcShiftLeft(dividend, partsCount, 1);
1015     assert(APInt::tcCompare(dividend, divisor, partsCount) >= 0);
1016   }
1017 
1018   /* Long division.  */
1019   for (bit = precision; bit; bit -= 1) {
1020     if (APInt::tcCompare(dividend, divisor, partsCount) >= 0) {
1021       APInt::tcSubtract(dividend, divisor, 0, partsCount);
1022       APInt::tcSetBit(lhsSignificand, bit - 1);
1023     }
1024 
1025     APInt::tcShiftLeft(dividend, partsCount, 1);
1026   }
1027 
1028   /* Figure out the lost fraction.  */
1029   int cmp = APInt::tcCompare(dividend, divisor, partsCount);
1030 
1031   if (cmp > 0)
1032     lost_fraction = lfMoreThanHalf;
1033   else if (cmp == 0)
1034     lost_fraction = lfExactlyHalf;
1035   else if (APInt::tcIsZero(dividend, partsCount))
1036     lost_fraction = lfExactlyZero;
1037   else
1038     lost_fraction = lfLessThanHalf;
1039 
1040   if (partsCount > 2)
1041     delete [] dividend;
1042 
1043   return lost_fraction;
1044 }
1045 
1046 unsigned int
significandMSB() const1047 APFloat::significandMSB() const
1048 {
1049   return APInt::tcMSB(significandParts(), partCount());
1050 }
1051 
1052 unsigned int
significandLSB() const1053 APFloat::significandLSB() const
1054 {
1055   return APInt::tcLSB(significandParts(), partCount());
1056 }
1057 
1058 /* Note that a zero result is NOT normalized to fcZero.  */
1059 lostFraction
shiftSignificandRight(unsigned int bits)1060 APFloat::shiftSignificandRight(unsigned int bits)
1061 {
1062   /* Our exponent should not overflow.  */
1063   assert((exponent_t) (exponent + bits) >= exponent);
1064 
1065   exponent += bits;
1066 
1067   return shiftRight(significandParts(), partCount(), bits);
1068 }
1069 
1070 /* Shift the significand left BITS bits, subtract BITS from its exponent.  */
1071 void
shiftSignificandLeft(unsigned int bits)1072 APFloat::shiftSignificandLeft(unsigned int bits)
1073 {
1074   assert(bits < semantics->precision);
1075 
1076   if (bits) {
1077     unsigned int partsCount = partCount();
1078 
1079     APInt::tcShiftLeft(significandParts(), partsCount, bits);
1080     exponent -= bits;
1081 
1082     assert(!APInt::tcIsZero(significandParts(), partsCount));
1083   }
1084 }
1085 
1086 APFloat::cmpResult
compareAbsoluteValue(const APFloat & rhs) const1087 APFloat::compareAbsoluteValue(const APFloat &rhs) const
1088 {
1089   int compare;
1090 
1091   assert(semantics == rhs.semantics);
1092   assert(category == fcNormal);
1093   assert(rhs.category == fcNormal);
1094 
1095   compare = exponent - rhs.exponent;
1096 
1097   /* If exponents are equal, do an unsigned bignum comparison of the
1098      significands.  */
1099   if (compare == 0)
1100     compare = APInt::tcCompare(significandParts(), rhs.significandParts(),
1101                                partCount());
1102 
1103   if (compare > 0)
1104     return cmpGreaterThan;
1105   else if (compare < 0)
1106     return cmpLessThan;
1107   else
1108     return cmpEqual;
1109 }
1110 
1111 /* Handle overflow.  Sign is preserved.  We either become infinity or
1112    the largest finite number.  */
1113 APFloat::opStatus
handleOverflow(roundingMode rounding_mode)1114 APFloat::handleOverflow(roundingMode rounding_mode)
1115 {
1116   /* Infinity?  */
1117   if (rounding_mode == rmNearestTiesToEven ||
1118       rounding_mode == rmNearestTiesToAway ||
1119       (rounding_mode == rmTowardPositive && !sign) ||
1120       (rounding_mode == rmTowardNegative && sign)) {
1121     category = fcInfinity;
1122     return (opStatus) (opOverflow | opInexact);
1123   }
1124 
1125   /* Otherwise we become the largest finite number.  */
1126   category = fcNormal;
1127   exponent = semantics->maxExponent;
1128   APInt::tcSetLeastSignificantBits(significandParts(), partCount(),
1129                                    semantics->precision);
1130 
1131   return opInexact;
1132 }
1133 
1134 /* Returns TRUE if, when truncating the current number, with BIT the
1135    new LSB, with the given lost fraction and rounding mode, the result
1136    would need to be rounded away from zero (i.e., by increasing the
1137    signficand).  This routine must work for fcZero of both signs, and
1138    fcNormal numbers.  */
1139 bool
roundAwayFromZero(roundingMode rounding_mode,lostFraction lost_fraction,unsigned int bit) const1140 APFloat::roundAwayFromZero(roundingMode rounding_mode,
1141                            lostFraction lost_fraction,
1142                            unsigned int bit) const
1143 {
1144   /* NaNs and infinities should not have lost fractions.  */
1145   assert(category == fcNormal || category == fcZero);
1146 
1147   /* Current callers never pass this so we don't handle it.  */
1148   assert(lost_fraction != lfExactlyZero);
1149 
1150   switch (rounding_mode) {
1151   default:
1152     llvm_unreachable(0);
1153 
1154   case rmNearestTiesToAway:
1155     return lost_fraction == lfExactlyHalf || lost_fraction == lfMoreThanHalf;
1156 
1157   case rmNearestTiesToEven:
1158     if (lost_fraction == lfMoreThanHalf)
1159       return true;
1160 
1161     /* Our zeroes don't have a significand to test.  */
1162     if (lost_fraction == lfExactlyHalf && category != fcZero)
1163       return APInt::tcExtractBit(significandParts(), bit);
1164 
1165     return false;
1166 
1167   case rmTowardZero:
1168     return false;
1169 
1170   case rmTowardPositive:
1171     return sign == false;
1172 
1173   case rmTowardNegative:
1174     return sign == true;
1175   }
1176 }
1177 
1178 APFloat::opStatus
normalize(roundingMode rounding_mode,lostFraction lost_fraction)1179 APFloat::normalize(roundingMode rounding_mode,
1180                    lostFraction lost_fraction)
1181 {
1182   unsigned int omsb;                /* One, not zero, based MSB.  */
1183   int exponentChange;
1184 
1185   if (category != fcNormal)
1186     return opOK;
1187 
1188   /* Before rounding normalize the exponent of fcNormal numbers.  */
1189   omsb = significandMSB() + 1;
1190 
1191   if (omsb) {
1192     /* OMSB is numbered from 1.  We want to place it in the integer
1193        bit numbered PRECISON if possible, with a compensating change in
1194        the exponent.  */
1195     exponentChange = omsb - semantics->precision;
1196 
1197     /* If the resulting exponent is too high, overflow according to
1198        the rounding mode.  */
1199     if (exponent + exponentChange > semantics->maxExponent)
1200       return handleOverflow(rounding_mode);
1201 
1202     /* Subnormal numbers have exponent minExponent, and their MSB
1203        is forced based on that.  */
1204     if (exponent + exponentChange < semantics->minExponent)
1205       exponentChange = semantics->minExponent - exponent;
1206 
1207     /* Shifting left is easy as we don't lose precision.  */
1208     if (exponentChange < 0) {
1209       assert(lost_fraction == lfExactlyZero);
1210 
1211       shiftSignificandLeft(-exponentChange);
1212 
1213       return opOK;
1214     }
1215 
1216     if (exponentChange > 0) {
1217       lostFraction lf;
1218 
1219       /* Shift right and capture any new lost fraction.  */
1220       lf = shiftSignificandRight(exponentChange);
1221 
1222       lost_fraction = combineLostFractions(lf, lost_fraction);
1223 
1224       /* Keep OMSB up-to-date.  */
1225       if (omsb > (unsigned) exponentChange)
1226         omsb -= exponentChange;
1227       else
1228         omsb = 0;
1229     }
1230   }
1231 
1232   /* Now round the number according to rounding_mode given the lost
1233      fraction.  */
1234 
1235   /* As specified in IEEE 754, since we do not trap we do not report
1236      underflow for exact results.  */
1237   if (lost_fraction == lfExactlyZero) {
1238     /* Canonicalize zeroes.  */
1239     if (omsb == 0)
1240       category = fcZero;
1241 
1242     return opOK;
1243   }
1244 
1245   /* Increment the significand if we're rounding away from zero.  */
1246   if (roundAwayFromZero(rounding_mode, lost_fraction, 0)) {
1247     if (omsb == 0)
1248       exponent = semantics->minExponent;
1249 
1250     incrementSignificand();
1251     omsb = significandMSB() + 1;
1252 
1253     /* Did the significand increment overflow?  */
1254     if (omsb == (unsigned) semantics->precision + 1) {
1255       /* Renormalize by incrementing the exponent and shifting our
1256          significand right one.  However if we already have the
1257          maximum exponent we overflow to infinity.  */
1258       if (exponent == semantics->maxExponent) {
1259         category = fcInfinity;
1260 
1261         return (opStatus) (opOverflow | opInexact);
1262       }
1263 
1264       shiftSignificandRight(1);
1265 
1266       return opInexact;
1267     }
1268   }
1269 
1270   /* The normal case - we were and are not denormal, and any
1271      significand increment above didn't overflow.  */
1272   if (omsb == semantics->precision)
1273     return opInexact;
1274 
1275   /* We have a non-zero denormal.  */
1276   assert(omsb < semantics->precision);
1277 
1278   /* Canonicalize zeroes.  */
1279   if (omsb == 0)
1280     category = fcZero;
1281 
1282   /* The fcZero case is a denormal that underflowed to zero.  */
1283   return (opStatus) (opUnderflow | opInexact);
1284 }
1285 
1286 APFloat::opStatus
addOrSubtractSpecials(const APFloat & rhs,bool subtract)1287 APFloat::addOrSubtractSpecials(const APFloat &rhs, bool subtract)
1288 {
1289   switch (convolve(category, rhs.category)) {
1290   default:
1291     llvm_unreachable(0);
1292 
1293   case convolve(fcNaN, fcZero):
1294   case convolve(fcNaN, fcNormal):
1295   case convolve(fcNaN, fcInfinity):
1296   case convolve(fcNaN, fcNaN):
1297   case convolve(fcNormal, fcZero):
1298   case convolve(fcInfinity, fcNormal):
1299   case convolve(fcInfinity, fcZero):
1300     return opOK;
1301 
1302   case convolve(fcZero, fcNaN):
1303   case convolve(fcNormal, fcNaN):
1304   case convolve(fcInfinity, fcNaN):
1305     category = fcNaN;
1306     copySignificand(rhs);
1307     return opOK;
1308 
1309   case convolve(fcNormal, fcInfinity):
1310   case convolve(fcZero, fcInfinity):
1311     category = fcInfinity;
1312     sign = rhs.sign ^ subtract;
1313     return opOK;
1314 
1315   case convolve(fcZero, fcNormal):
1316     assign(rhs);
1317     sign = rhs.sign ^ subtract;
1318     return opOK;
1319 
1320   case convolve(fcZero, fcZero):
1321     /* Sign depends on rounding mode; handled by caller.  */
1322     return opOK;
1323 
1324   case convolve(fcInfinity, fcInfinity):
1325     /* Differently signed infinities can only be validly
1326        subtracted.  */
1327     if (((sign ^ rhs.sign)!=0) != subtract) {
1328       makeNaN();
1329       return opInvalidOp;
1330     }
1331 
1332     return opOK;
1333 
1334   case convolve(fcNormal, fcNormal):
1335     return opDivByZero;
1336   }
1337 }
1338 
1339 /* Add or subtract two normal numbers.  */
1340 lostFraction
addOrSubtractSignificand(const APFloat & rhs,bool subtract)1341 APFloat::addOrSubtractSignificand(const APFloat &rhs, bool subtract)
1342 {
1343   integerPart carry;
1344   lostFraction lost_fraction;
1345   int bits;
1346 
1347   /* Determine if the operation on the absolute values is effectively
1348      an addition or subtraction.  */
1349   subtract ^= (sign ^ rhs.sign) ? true : false;
1350 
1351   /* Are we bigger exponent-wise than the RHS?  */
1352   bits = exponent - rhs.exponent;
1353 
1354   /* Subtraction is more subtle than one might naively expect.  */
1355   if (subtract) {
1356     APFloat temp_rhs(rhs);
1357     bool reverse;
1358 
1359     if (bits == 0) {
1360       reverse = compareAbsoluteValue(temp_rhs) == cmpLessThan;
1361       lost_fraction = lfExactlyZero;
1362     } else if (bits > 0) {
1363       lost_fraction = temp_rhs.shiftSignificandRight(bits - 1);
1364       shiftSignificandLeft(1);
1365       reverse = false;
1366     } else {
1367       lost_fraction = shiftSignificandRight(-bits - 1);
1368       temp_rhs.shiftSignificandLeft(1);
1369       reverse = true;
1370     }
1371 
1372     if (reverse) {
1373       carry = temp_rhs.subtractSignificand
1374         (*this, lost_fraction != lfExactlyZero);
1375       copySignificand(temp_rhs);
1376       sign = !sign;
1377     } else {
1378       carry = subtractSignificand
1379         (temp_rhs, lost_fraction != lfExactlyZero);
1380     }
1381 
1382     /* Invert the lost fraction - it was on the RHS and
1383        subtracted.  */
1384     if (lost_fraction == lfLessThanHalf)
1385       lost_fraction = lfMoreThanHalf;
1386     else if (lost_fraction == lfMoreThanHalf)
1387       lost_fraction = lfLessThanHalf;
1388 
1389     /* The code above is intended to ensure that no borrow is
1390        necessary.  */
1391     assert(!carry);
1392   } else {
1393     if (bits > 0) {
1394       APFloat temp_rhs(rhs);
1395 
1396       lost_fraction = temp_rhs.shiftSignificandRight(bits);
1397       carry = addSignificand(temp_rhs);
1398     } else {
1399       lost_fraction = shiftSignificandRight(-bits);
1400       carry = addSignificand(rhs);
1401     }
1402 
1403     /* We have a guard bit; generating a carry cannot happen.  */
1404     assert(!carry);
1405   }
1406 
1407   return lost_fraction;
1408 }
1409 
1410 APFloat::opStatus
multiplySpecials(const APFloat & rhs)1411 APFloat::multiplySpecials(const APFloat &rhs)
1412 {
1413   switch (convolve(category, rhs.category)) {
1414   default:
1415     llvm_unreachable(0);
1416 
1417   case convolve(fcNaN, fcZero):
1418   case convolve(fcNaN, fcNormal):
1419   case convolve(fcNaN, fcInfinity):
1420   case convolve(fcNaN, fcNaN):
1421     return opOK;
1422 
1423   case convolve(fcZero, fcNaN):
1424   case convolve(fcNormal, fcNaN):
1425   case convolve(fcInfinity, fcNaN):
1426     category = fcNaN;
1427     copySignificand(rhs);
1428     return opOK;
1429 
1430   case convolve(fcNormal, fcInfinity):
1431   case convolve(fcInfinity, fcNormal):
1432   case convolve(fcInfinity, fcInfinity):
1433     category = fcInfinity;
1434     return opOK;
1435 
1436   case convolve(fcZero, fcNormal):
1437   case convolve(fcNormal, fcZero):
1438   case convolve(fcZero, fcZero):
1439     category = fcZero;
1440     return opOK;
1441 
1442   case convolve(fcZero, fcInfinity):
1443   case convolve(fcInfinity, fcZero):
1444     makeNaN();
1445     return opInvalidOp;
1446 
1447   case convolve(fcNormal, fcNormal):
1448     return opOK;
1449   }
1450 }
1451 
1452 APFloat::opStatus
divideSpecials(const APFloat & rhs)1453 APFloat::divideSpecials(const APFloat &rhs)
1454 {
1455   switch (convolve(category, rhs.category)) {
1456   default:
1457     llvm_unreachable(0);
1458 
1459   case convolve(fcNaN, fcZero):
1460   case convolve(fcNaN, fcNormal):
1461   case convolve(fcNaN, fcInfinity):
1462   case convolve(fcNaN, fcNaN):
1463   case convolve(fcInfinity, fcZero):
1464   case convolve(fcInfinity, fcNormal):
1465   case convolve(fcZero, fcInfinity):
1466   case convolve(fcZero, fcNormal):
1467     return opOK;
1468 
1469   case convolve(fcZero, fcNaN):
1470   case convolve(fcNormal, fcNaN):
1471   case convolve(fcInfinity, fcNaN):
1472     category = fcNaN;
1473     copySignificand(rhs);
1474     return opOK;
1475 
1476   case convolve(fcNormal, fcInfinity):
1477     category = fcZero;
1478     return opOK;
1479 
1480   case convolve(fcNormal, fcZero):
1481     category = fcInfinity;
1482     return opDivByZero;
1483 
1484   case convolve(fcInfinity, fcInfinity):
1485   case convolve(fcZero, fcZero):
1486     makeNaN();
1487     return opInvalidOp;
1488 
1489   case convolve(fcNormal, fcNormal):
1490     return opOK;
1491   }
1492 }
1493 
1494 APFloat::opStatus
modSpecials(const APFloat & rhs)1495 APFloat::modSpecials(const APFloat &rhs)
1496 {
1497   switch (convolve(category, rhs.category)) {
1498   default:
1499     llvm_unreachable(0);
1500 
1501   case convolve(fcNaN, fcZero):
1502   case convolve(fcNaN, fcNormal):
1503   case convolve(fcNaN, fcInfinity):
1504   case convolve(fcNaN, fcNaN):
1505   case convolve(fcZero, fcInfinity):
1506   case convolve(fcZero, fcNormal):
1507   case convolve(fcNormal, fcInfinity):
1508     return opOK;
1509 
1510   case convolve(fcZero, fcNaN):
1511   case convolve(fcNormal, fcNaN):
1512   case convolve(fcInfinity, fcNaN):
1513     category = fcNaN;
1514     copySignificand(rhs);
1515     return opOK;
1516 
1517   case convolve(fcNormal, fcZero):
1518   case convolve(fcInfinity, fcZero):
1519   case convolve(fcInfinity, fcNormal):
1520   case convolve(fcInfinity, fcInfinity):
1521   case convolve(fcZero, fcZero):
1522     makeNaN();
1523     return opInvalidOp;
1524 
1525   case convolve(fcNormal, fcNormal):
1526     return opOK;
1527   }
1528 }
1529 
1530 /* Change sign.  */
1531 void
changeSign()1532 APFloat::changeSign()
1533 {
1534   /* Look mummy, this one's easy.  */
1535   sign = !sign;
1536 }
1537 
1538 void
clearSign()1539 APFloat::clearSign()
1540 {
1541   /* So is this one. */
1542   sign = 0;
1543 }
1544 
1545 void
copySign(const APFloat & rhs)1546 APFloat::copySign(const APFloat &rhs)
1547 {
1548   /* And this one. */
1549   sign = rhs.sign;
1550 }
1551 
1552 /* Normalized addition or subtraction.  */
1553 APFloat::opStatus
addOrSubtract(const APFloat & rhs,roundingMode rounding_mode,bool subtract)1554 APFloat::addOrSubtract(const APFloat &rhs, roundingMode rounding_mode,
1555                        bool subtract)
1556 {
1557   opStatus fs;
1558 
1559   assertArithmeticOK(*semantics);
1560 
1561   fs = addOrSubtractSpecials(rhs, subtract);
1562 
1563   /* This return code means it was not a simple case.  */
1564   if (fs == opDivByZero) {
1565     lostFraction lost_fraction;
1566 
1567     lost_fraction = addOrSubtractSignificand(rhs, subtract);
1568     fs = normalize(rounding_mode, lost_fraction);
1569 
1570     /* Can only be zero if we lost no fraction.  */
1571     assert(category != fcZero || lost_fraction == lfExactlyZero);
1572   }
1573 
1574   /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1575      positive zero unless rounding to minus infinity, except that
1576      adding two like-signed zeroes gives that zero.  */
1577   if (category == fcZero) {
1578     if (rhs.category != fcZero || (sign == rhs.sign) == subtract)
1579       sign = (rounding_mode == rmTowardNegative);
1580   }
1581 
1582   return fs;
1583 }
1584 
1585 /* Normalized addition.  */
1586 APFloat::opStatus
add(const APFloat & rhs,roundingMode rounding_mode)1587 APFloat::add(const APFloat &rhs, roundingMode rounding_mode)
1588 {
1589   return addOrSubtract(rhs, rounding_mode, false);
1590 }
1591 
1592 /* Normalized subtraction.  */
1593 APFloat::opStatus
subtract(const APFloat & rhs,roundingMode rounding_mode)1594 APFloat::subtract(const APFloat &rhs, roundingMode rounding_mode)
1595 {
1596   return addOrSubtract(rhs, rounding_mode, true);
1597 }
1598 
1599 /* Normalized multiply.  */
1600 APFloat::opStatus
multiply(const APFloat & rhs,roundingMode rounding_mode)1601 APFloat::multiply(const APFloat &rhs, roundingMode rounding_mode)
1602 {
1603   opStatus fs;
1604 
1605   assertArithmeticOK(*semantics);
1606   sign ^= rhs.sign;
1607   fs = multiplySpecials(rhs);
1608 
1609   if (category == fcNormal) {
1610     lostFraction lost_fraction = multiplySignificand(rhs, 0);
1611     fs = normalize(rounding_mode, lost_fraction);
1612     if (lost_fraction != lfExactlyZero)
1613       fs = (opStatus) (fs | opInexact);
1614   }
1615 
1616   return fs;
1617 }
1618 
1619 /* Normalized divide.  */
1620 APFloat::opStatus
divide(const APFloat & rhs,roundingMode rounding_mode)1621 APFloat::divide(const APFloat &rhs, roundingMode rounding_mode)
1622 {
1623   opStatus fs;
1624 
1625   assertArithmeticOK(*semantics);
1626   sign ^= rhs.sign;
1627   fs = divideSpecials(rhs);
1628 
1629   if (category == fcNormal) {
1630     lostFraction lost_fraction = divideSignificand(rhs);
1631     fs = normalize(rounding_mode, lost_fraction);
1632     if (lost_fraction != lfExactlyZero)
1633       fs = (opStatus) (fs | opInexact);
1634   }
1635 
1636   return fs;
1637 }
1638 
1639 /* Normalized remainder.  This is not currently correct in all cases.  */
1640 APFloat::opStatus
remainder(const APFloat & rhs)1641 APFloat::remainder(const APFloat &rhs)
1642 {
1643   opStatus fs;
1644   APFloat V = *this;
1645   unsigned int origSign = sign;
1646 
1647   assertArithmeticOK(*semantics);
1648   fs = V.divide(rhs, rmNearestTiesToEven);
1649   if (fs == opDivByZero)
1650     return fs;
1651 
1652   int parts = partCount();
1653   integerPart *x = new integerPart[parts];
1654   bool ignored;
1655   fs = V.convertToInteger(x, parts * integerPartWidth, true,
1656                           rmNearestTiesToEven, &ignored);
1657   if (fs==opInvalidOp)
1658     return fs;
1659 
1660   fs = V.convertFromZeroExtendedInteger(x, parts * integerPartWidth, true,
1661                                         rmNearestTiesToEven);
1662   assert(fs==opOK);   // should always work
1663 
1664   fs = V.multiply(rhs, rmNearestTiesToEven);
1665   assert(fs==opOK || fs==opInexact);   // should not overflow or underflow
1666 
1667   fs = subtract(V, rmNearestTiesToEven);
1668   assert(fs==opOK || fs==opInexact);   // likewise
1669 
1670   if (isZero())
1671     sign = origSign;    // IEEE754 requires this
1672   delete[] x;
1673   return fs;
1674 }
1675 
1676 /* Normalized llvm frem (C fmod).
1677    This is not currently correct in all cases.  */
1678 APFloat::opStatus
mod(const APFloat & rhs,roundingMode rounding_mode)1679 APFloat::mod(const APFloat &rhs, roundingMode rounding_mode)
1680 {
1681   opStatus fs;
1682   assertArithmeticOK(*semantics);
1683   fs = modSpecials(rhs);
1684 
1685   if (category == fcNormal && rhs.category == fcNormal) {
1686     APFloat V = *this;
1687     unsigned int origSign = sign;
1688 
1689     fs = V.divide(rhs, rmNearestTiesToEven);
1690     if (fs == opDivByZero)
1691       return fs;
1692 
1693     int parts = partCount();
1694     integerPart *x = new integerPart[parts];
1695     bool ignored;
1696     fs = V.convertToInteger(x, parts * integerPartWidth, true,
1697                             rmTowardZero, &ignored);
1698     if (fs==opInvalidOp)
1699       return fs;
1700 
1701     fs = V.convertFromZeroExtendedInteger(x, parts * integerPartWidth, true,
1702                                           rmNearestTiesToEven);
1703     assert(fs==opOK);   // should always work
1704 
1705     fs = V.multiply(rhs, rounding_mode);
1706     assert(fs==opOK || fs==opInexact);   // should not overflow or underflow
1707 
1708     fs = subtract(V, rounding_mode);
1709     assert(fs==opOK || fs==opInexact);   // likewise
1710 
1711     if (isZero())
1712       sign = origSign;    // IEEE754 requires this
1713     delete[] x;
1714   }
1715   return fs;
1716 }
1717 
1718 /* Normalized fused-multiply-add.  */
1719 APFloat::opStatus
fusedMultiplyAdd(const APFloat & multiplicand,const APFloat & addend,roundingMode rounding_mode)1720 APFloat::fusedMultiplyAdd(const APFloat &multiplicand,
1721                           const APFloat &addend,
1722                           roundingMode rounding_mode)
1723 {
1724   opStatus fs;
1725 
1726   assertArithmeticOK(*semantics);
1727 
1728   /* Post-multiplication sign, before addition.  */
1729   sign ^= multiplicand.sign;
1730 
1731   /* If and only if all arguments are normal do we need to do an
1732      extended-precision calculation.  */
1733   if (category == fcNormal &&
1734       multiplicand.category == fcNormal &&
1735       addend.category == fcNormal) {
1736     lostFraction lost_fraction;
1737 
1738     lost_fraction = multiplySignificand(multiplicand, &addend);
1739     fs = normalize(rounding_mode, lost_fraction);
1740     if (lost_fraction != lfExactlyZero)
1741       fs = (opStatus) (fs | opInexact);
1742 
1743     /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1744        positive zero unless rounding to minus infinity, except that
1745        adding two like-signed zeroes gives that zero.  */
1746     if (category == fcZero && sign != addend.sign)
1747       sign = (rounding_mode == rmTowardNegative);
1748   } else {
1749     fs = multiplySpecials(multiplicand);
1750 
1751     /* FS can only be opOK or opInvalidOp.  There is no more work
1752        to do in the latter case.  The IEEE-754R standard says it is
1753        implementation-defined in this case whether, if ADDEND is a
1754        quiet NaN, we raise invalid op; this implementation does so.
1755 
1756        If we need to do the addition we can do so with normal
1757        precision.  */
1758     if (fs == opOK)
1759       fs = addOrSubtract(addend, rounding_mode, false);
1760   }
1761 
1762   return fs;
1763 }
1764 
1765 /* Comparison requires normalized numbers.  */
1766 APFloat::cmpResult
compare(const APFloat & rhs) const1767 APFloat::compare(const APFloat &rhs) const
1768 {
1769   cmpResult result;
1770 
1771   assertArithmeticOK(*semantics);
1772   assert(semantics == rhs.semantics);
1773 
1774   switch (convolve(category, rhs.category)) {
1775   default:
1776     llvm_unreachable(0);
1777 
1778   case convolve(fcNaN, fcZero):
1779   case convolve(fcNaN, fcNormal):
1780   case convolve(fcNaN, fcInfinity):
1781   case convolve(fcNaN, fcNaN):
1782   case convolve(fcZero, fcNaN):
1783   case convolve(fcNormal, fcNaN):
1784   case convolve(fcInfinity, fcNaN):
1785     return cmpUnordered;
1786 
1787   case convolve(fcInfinity, fcNormal):
1788   case convolve(fcInfinity, fcZero):
1789   case convolve(fcNormal, fcZero):
1790     if (sign)
1791       return cmpLessThan;
1792     else
1793       return cmpGreaterThan;
1794 
1795   case convolve(fcNormal, fcInfinity):
1796   case convolve(fcZero, fcInfinity):
1797   case convolve(fcZero, fcNormal):
1798     if (rhs.sign)
1799       return cmpGreaterThan;
1800     else
1801       return cmpLessThan;
1802 
1803   case convolve(fcInfinity, fcInfinity):
1804     if (sign == rhs.sign)
1805       return cmpEqual;
1806     else if (sign)
1807       return cmpLessThan;
1808     else
1809       return cmpGreaterThan;
1810 
1811   case convolve(fcZero, fcZero):
1812     return cmpEqual;
1813 
1814   case convolve(fcNormal, fcNormal):
1815     break;
1816   }
1817 
1818   /* Two normal numbers.  Do they have the same sign?  */
1819   if (sign != rhs.sign) {
1820     if (sign)
1821       result = cmpLessThan;
1822     else
1823       result = cmpGreaterThan;
1824   } else {
1825     /* Compare absolute values; invert result if negative.  */
1826     result = compareAbsoluteValue(rhs);
1827 
1828     if (sign) {
1829       if (result == cmpLessThan)
1830         result = cmpGreaterThan;
1831       else if (result == cmpGreaterThan)
1832         result = cmpLessThan;
1833     }
1834   }
1835 
1836   return result;
1837 }
1838 
1839 /// APFloat::convert - convert a value of one floating point type to another.
1840 /// The return value corresponds to the IEEE754 exceptions.  *losesInfo
1841 /// records whether the transformation lost information, i.e. whether
1842 /// converting the result back to the original type will produce the
1843 /// original value (this is almost the same as return value==fsOK, but there
1844 /// are edge cases where this is not so).
1845 
1846 APFloat::opStatus
convert(const fltSemantics & toSemantics,roundingMode rounding_mode,bool * losesInfo)1847 APFloat::convert(const fltSemantics &toSemantics,
1848                  roundingMode rounding_mode, bool *losesInfo)
1849 {
1850   lostFraction lostFraction;
1851   unsigned int newPartCount, oldPartCount;
1852   opStatus fs;
1853 
1854   assertArithmeticOK(*semantics);
1855   assertArithmeticOK(toSemantics);
1856   lostFraction = lfExactlyZero;
1857   newPartCount = partCountForBits(toSemantics.precision + 1);
1858   oldPartCount = partCount();
1859 
1860   /* Handle storage complications.  If our new form is wider,
1861      re-allocate our bit pattern into wider storage.  If it is
1862      narrower, we ignore the excess parts, but if narrowing to a
1863      single part we need to free the old storage.
1864      Be careful not to reference significandParts for zeroes
1865      and infinities, since it aborts.  */
1866   if (newPartCount > oldPartCount) {
1867     integerPart *newParts;
1868     newParts = new integerPart[newPartCount];
1869     APInt::tcSet(newParts, 0, newPartCount);
1870     if (category==fcNormal || category==fcNaN)
1871       APInt::tcAssign(newParts, significandParts(), oldPartCount);
1872     freeSignificand();
1873     significand.parts = newParts;
1874   } else if (newPartCount < oldPartCount) {
1875     /* Capture any lost fraction through truncation of parts so we get
1876        correct rounding whilst normalizing.  */
1877     if (category==fcNormal)
1878       lostFraction = lostFractionThroughTruncation
1879         (significandParts(), oldPartCount, toSemantics.precision);
1880     if (newPartCount == 1) {
1881         integerPart newPart = 0;
1882         if (category==fcNormal || category==fcNaN)
1883           newPart = significandParts()[0];
1884         freeSignificand();
1885         significand.part = newPart;
1886     }
1887   }
1888 
1889   if (category == fcNormal) {
1890     /* Re-interpret our bit-pattern.  */
1891     exponent += toSemantics.precision - semantics->precision;
1892     semantics = &toSemantics;
1893     fs = normalize(rounding_mode, lostFraction);
1894     *losesInfo = (fs != opOK);
1895   } else if (category == fcNaN) {
1896     int shift = toSemantics.precision - semantics->precision;
1897     // Do this now so significandParts gets the right answer
1898     const fltSemantics *oldSemantics = semantics;
1899     semantics = &toSemantics;
1900     *losesInfo = false;
1901     // No normalization here, just truncate
1902     if (shift>0)
1903       APInt::tcShiftLeft(significandParts(), newPartCount, shift);
1904     else if (shift < 0) {
1905       unsigned ushift = -shift;
1906       // Figure out if we are losing information.  This happens
1907       // if are shifting out something other than 0s, or if the x87 long
1908       // double input did not have its integer bit set (pseudo-NaN), or if the
1909       // x87 long double input did not have its QNan bit set (because the x87
1910       // hardware sets this bit when converting a lower-precision NaN to
1911       // x87 long double).
1912       if (APInt::tcLSB(significandParts(), newPartCount) < ushift)
1913         *losesInfo = true;
1914       if (oldSemantics == &APFloat::x87DoubleExtended &&
1915           (!(*significandParts() & 0x8000000000000000ULL) ||
1916            !(*significandParts() & 0x4000000000000000ULL)))
1917         *losesInfo = true;
1918       APInt::tcShiftRight(significandParts(), newPartCount, ushift);
1919     }
1920     // gcc forces the Quiet bit on, which means (float)(double)(float_sNan)
1921     // does not give you back the same bits.  This is dubious, and we
1922     // don't currently do it.  You're really supposed to get
1923     // an invalid operation signal at runtime, but nobody does that.
1924     fs = opOK;
1925   } else {
1926     semantics = &toSemantics;
1927     fs = opOK;
1928     *losesInfo = false;
1929   }
1930 
1931   return fs;
1932 }
1933 
1934 /* Convert a floating point number to an integer according to the
1935    rounding mode.  If the rounded integer value is out of range this
1936    returns an invalid operation exception and the contents of the
1937    destination parts are unspecified.  If the rounded value is in
1938    range but the floating point number is not the exact integer, the C
1939    standard doesn't require an inexact exception to be raised.  IEEE
1940    854 does require it so we do that.
1941 
1942    Note that for conversions to integer type the C standard requires
1943    round-to-zero to always be used.  */
1944 APFloat::opStatus
convertToSignExtendedInteger(integerPart * parts,unsigned int width,bool isSigned,roundingMode rounding_mode,bool * isExact) const1945 APFloat::convertToSignExtendedInteger(integerPart *parts, unsigned int width,
1946                                       bool isSigned,
1947                                       roundingMode rounding_mode,
1948                                       bool *isExact) const
1949 {
1950   lostFraction lost_fraction;
1951   const integerPart *src;
1952   unsigned int dstPartsCount, truncatedBits;
1953 
1954   assertArithmeticOK(*semantics);
1955 
1956   *isExact = false;
1957 
1958   /* Handle the three special cases first.  */
1959   if (category == fcInfinity || category == fcNaN)
1960     return opInvalidOp;
1961 
1962   dstPartsCount = partCountForBits(width);
1963 
1964   if (category == fcZero) {
1965     APInt::tcSet(parts, 0, dstPartsCount);
1966     // Negative zero can't be represented as an int.
1967     *isExact = !sign;
1968     return opOK;
1969   }
1970 
1971   src = significandParts();
1972 
1973   /* Step 1: place our absolute value, with any fraction truncated, in
1974      the destination.  */
1975   if (exponent < 0) {
1976     /* Our absolute value is less than one; truncate everything.  */
1977     APInt::tcSet(parts, 0, dstPartsCount);
1978     /* For exponent -1 the integer bit represents .5, look at that.
1979        For smaller exponents leftmost truncated bit is 0. */
1980     truncatedBits = semantics->precision -1U - exponent;
1981   } else {
1982     /* We want the most significant (exponent + 1) bits; the rest are
1983        truncated.  */
1984     unsigned int bits = exponent + 1U;
1985 
1986     /* Hopelessly large in magnitude?  */
1987     if (bits > width)
1988       return opInvalidOp;
1989 
1990     if (bits < semantics->precision) {
1991       /* We truncate (semantics->precision - bits) bits.  */
1992       truncatedBits = semantics->precision - bits;
1993       APInt::tcExtract(parts, dstPartsCount, src, bits, truncatedBits);
1994     } else {
1995       /* We want at least as many bits as are available.  */
1996       APInt::tcExtract(parts, dstPartsCount, src, semantics->precision, 0);
1997       APInt::tcShiftLeft(parts, dstPartsCount, bits - semantics->precision);
1998       truncatedBits = 0;
1999     }
2000   }
2001 
2002   /* Step 2: work out any lost fraction, and increment the absolute
2003      value if we would round away from zero.  */
2004   if (truncatedBits) {
2005     lost_fraction = lostFractionThroughTruncation(src, partCount(),
2006                                                   truncatedBits);
2007     if (lost_fraction != lfExactlyZero &&
2008         roundAwayFromZero(rounding_mode, lost_fraction, truncatedBits)) {
2009       if (APInt::tcIncrement(parts, dstPartsCount))
2010         return opInvalidOp;     /* Overflow.  */
2011     }
2012   } else {
2013     lost_fraction = lfExactlyZero;
2014   }
2015 
2016   /* Step 3: check if we fit in the destination.  */
2017   unsigned int omsb = APInt::tcMSB(parts, dstPartsCount) + 1;
2018 
2019   if (sign) {
2020     if (!isSigned) {
2021       /* Negative numbers cannot be represented as unsigned.  */
2022       if (omsb != 0)
2023         return opInvalidOp;
2024     } else {
2025       /* It takes omsb bits to represent the unsigned integer value.
2026          We lose a bit for the sign, but care is needed as the
2027          maximally negative integer is a special case.  */
2028       if (omsb == width && APInt::tcLSB(parts, dstPartsCount) + 1 != omsb)
2029         return opInvalidOp;
2030 
2031       /* This case can happen because of rounding.  */
2032       if (omsb > width)
2033         return opInvalidOp;
2034     }
2035 
2036     APInt::tcNegate (parts, dstPartsCount);
2037   } else {
2038     if (omsb >= width + !isSigned)
2039       return opInvalidOp;
2040   }
2041 
2042   if (lost_fraction == lfExactlyZero) {
2043     *isExact = true;
2044     return opOK;
2045   } else
2046     return opInexact;
2047 }
2048 
2049 /* Same as convertToSignExtendedInteger, except we provide
2050    deterministic values in case of an invalid operation exception,
2051    namely zero for NaNs and the minimal or maximal value respectively
2052    for underflow or overflow.
2053    The *isExact output tells whether the result is exact, in the sense
2054    that converting it back to the original floating point type produces
2055    the original value.  This is almost equivalent to result==opOK,
2056    except for negative zeroes.
2057 */
2058 APFloat::opStatus
convertToInteger(integerPart * parts,unsigned int width,bool isSigned,roundingMode rounding_mode,bool * isExact) const2059 APFloat::convertToInteger(integerPart *parts, unsigned int width,
2060                           bool isSigned,
2061                           roundingMode rounding_mode, bool *isExact) const
2062 {
2063   opStatus fs;
2064 
2065   fs = convertToSignExtendedInteger(parts, width, isSigned, rounding_mode,
2066                                     isExact);
2067 
2068   if (fs == opInvalidOp) {
2069     unsigned int bits, dstPartsCount;
2070 
2071     dstPartsCount = partCountForBits(width);
2072 
2073     if (category == fcNaN)
2074       bits = 0;
2075     else if (sign)
2076       bits = isSigned;
2077     else
2078       bits = width - isSigned;
2079 
2080     APInt::tcSetLeastSignificantBits(parts, dstPartsCount, bits);
2081     if (sign && isSigned)
2082       APInt::tcShiftLeft(parts, dstPartsCount, width - 1);
2083   }
2084 
2085   return fs;
2086 }
2087 
2088 /* Same as convertToInteger(integerPart*, ...), except the result is returned in
2089    an APSInt, whose initial bit-width and signed-ness are used to determine the
2090    precision of the conversion.
2091  */
2092 APFloat::opStatus
convertToInteger(APSInt & result,roundingMode rounding_mode,bool * isExact) const2093 APFloat::convertToInteger(APSInt &result,
2094                           roundingMode rounding_mode, bool *isExact) const
2095 {
2096   unsigned bitWidth = result.getBitWidth();
2097   SmallVector<uint64_t, 4> parts(result.getNumWords());
2098   opStatus status = convertToInteger(
2099     parts.data(), bitWidth, result.isSigned(), rounding_mode, isExact);
2100   // Keeps the original signed-ness.
2101   result = APInt(bitWidth, parts);
2102   return status;
2103 }
2104 
2105 /* Convert an unsigned integer SRC to a floating point number,
2106    rounding according to ROUNDING_MODE.  The sign of the floating
2107    point number is not modified.  */
2108 APFloat::opStatus
convertFromUnsignedParts(const integerPart * src,unsigned int srcCount,roundingMode rounding_mode)2109 APFloat::convertFromUnsignedParts(const integerPart *src,
2110                                   unsigned int srcCount,
2111                                   roundingMode rounding_mode)
2112 {
2113   unsigned int omsb, precision, dstCount;
2114   integerPart *dst;
2115   lostFraction lost_fraction;
2116 
2117   assertArithmeticOK(*semantics);
2118   category = fcNormal;
2119   omsb = APInt::tcMSB(src, srcCount) + 1;
2120   dst = significandParts();
2121   dstCount = partCount();
2122   precision = semantics->precision;
2123 
2124   /* We want the most significant PRECISON bits of SRC.  There may not
2125      be that many; extract what we can.  */
2126   if (precision <= omsb) {
2127     exponent = omsb - 1;
2128     lost_fraction = lostFractionThroughTruncation(src, srcCount,
2129                                                   omsb - precision);
2130     APInt::tcExtract(dst, dstCount, src, precision, omsb - precision);
2131   } else {
2132     exponent = precision - 1;
2133     lost_fraction = lfExactlyZero;
2134     APInt::tcExtract(dst, dstCount, src, omsb, 0);
2135   }
2136 
2137   return normalize(rounding_mode, lost_fraction);
2138 }
2139 
2140 APFloat::opStatus
convertFromAPInt(const APInt & Val,bool isSigned,roundingMode rounding_mode)2141 APFloat::convertFromAPInt(const APInt &Val,
2142                           bool isSigned,
2143                           roundingMode rounding_mode)
2144 {
2145   unsigned int partCount = Val.getNumWords();
2146   APInt api = Val;
2147 
2148   sign = false;
2149   if (isSigned && api.isNegative()) {
2150     sign = true;
2151     api = -api;
2152   }
2153 
2154   return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode);
2155 }
2156 
2157 /* Convert a two's complement integer SRC to a floating point number,
2158    rounding according to ROUNDING_MODE.  ISSIGNED is true if the
2159    integer is signed, in which case it must be sign-extended.  */
2160 APFloat::opStatus
convertFromSignExtendedInteger(const integerPart * src,unsigned int srcCount,bool isSigned,roundingMode rounding_mode)2161 APFloat::convertFromSignExtendedInteger(const integerPart *src,
2162                                         unsigned int srcCount,
2163                                         bool isSigned,
2164                                         roundingMode rounding_mode)
2165 {
2166   opStatus status;
2167 
2168   assertArithmeticOK(*semantics);
2169   if (isSigned &&
2170       APInt::tcExtractBit(src, srcCount * integerPartWidth - 1)) {
2171     integerPart *copy;
2172 
2173     /* If we're signed and negative negate a copy.  */
2174     sign = true;
2175     copy = new integerPart[srcCount];
2176     APInt::tcAssign(copy, src, srcCount);
2177     APInt::tcNegate(copy, srcCount);
2178     status = convertFromUnsignedParts(copy, srcCount, rounding_mode);
2179     delete [] copy;
2180   } else {
2181     sign = false;
2182     status = convertFromUnsignedParts(src, srcCount, rounding_mode);
2183   }
2184 
2185   return status;
2186 }
2187 
2188 /* FIXME: should this just take a const APInt reference?  */
2189 APFloat::opStatus
convertFromZeroExtendedInteger(const integerPart * parts,unsigned int width,bool isSigned,roundingMode rounding_mode)2190 APFloat::convertFromZeroExtendedInteger(const integerPart *parts,
2191                                         unsigned int width, bool isSigned,
2192                                         roundingMode rounding_mode)
2193 {
2194   unsigned int partCount = partCountForBits(width);
2195   APInt api = APInt(width, makeArrayRef(parts, partCount));
2196 
2197   sign = false;
2198   if (isSigned && APInt::tcExtractBit(parts, width - 1)) {
2199     sign = true;
2200     api = -api;
2201   }
2202 
2203   return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode);
2204 }
2205 
2206 APFloat::opStatus
convertFromHexadecimalString(StringRef s,roundingMode rounding_mode)2207 APFloat::convertFromHexadecimalString(StringRef s, roundingMode rounding_mode)
2208 {
2209   lostFraction lost_fraction = lfExactlyZero;
2210   integerPart *significand;
2211   unsigned int bitPos, partsCount;
2212   StringRef::iterator dot, firstSignificantDigit;
2213 
2214   zeroSignificand();
2215   exponent = 0;
2216   category = fcNormal;
2217 
2218   significand = significandParts();
2219   partsCount = partCount();
2220   bitPos = partsCount * integerPartWidth;
2221 
2222   /* Skip leading zeroes and any (hexa)decimal point.  */
2223   StringRef::iterator begin = s.begin();
2224   StringRef::iterator end = s.end();
2225   StringRef::iterator p = skipLeadingZeroesAndAnyDot(begin, end, &dot);
2226   firstSignificantDigit = p;
2227 
2228   for (; p != end;) {
2229     integerPart hex_value;
2230 
2231     if (*p == '.') {
2232       assert(dot == end && "String contains multiple dots");
2233       dot = p++;
2234       if (p == end) {
2235         break;
2236       }
2237     }
2238 
2239     hex_value = hexDigitValue(*p);
2240     if (hex_value == -1U) {
2241       break;
2242     }
2243 
2244     p++;
2245 
2246     if (p == end) {
2247       break;
2248     } else {
2249       /* Store the number whilst 4-bit nibbles remain.  */
2250       if (bitPos) {
2251         bitPos -= 4;
2252         hex_value <<= bitPos % integerPartWidth;
2253         significand[bitPos / integerPartWidth] |= hex_value;
2254       } else {
2255         lost_fraction = trailingHexadecimalFraction(p, end, hex_value);
2256         while (p != end && hexDigitValue(*p) != -1U)
2257           p++;
2258         break;
2259       }
2260     }
2261   }
2262 
2263   /* Hex floats require an exponent but not a hexadecimal point.  */
2264   assert(p != end && "Hex strings require an exponent");
2265   assert((*p == 'p' || *p == 'P') && "Invalid character in significand");
2266   assert(p != begin && "Significand has no digits");
2267   assert((dot == end || p - begin != 1) && "Significand has no digits");
2268 
2269   /* Ignore the exponent if we are zero.  */
2270   if (p != firstSignificantDigit) {
2271     int expAdjustment;
2272 
2273     /* Implicit hexadecimal point?  */
2274     if (dot == end)
2275       dot = p;
2276 
2277     /* Calculate the exponent adjustment implicit in the number of
2278        significant digits.  */
2279     expAdjustment = static_cast<int>(dot - firstSignificantDigit);
2280     if (expAdjustment < 0)
2281       expAdjustment++;
2282     expAdjustment = expAdjustment * 4 - 1;
2283 
2284     /* Adjust for writing the significand starting at the most
2285        significant nibble.  */
2286     expAdjustment += semantics->precision;
2287     expAdjustment -= partsCount * integerPartWidth;
2288 
2289     /* Adjust for the given exponent.  */
2290     exponent = totalExponent(p + 1, end, expAdjustment);
2291   }
2292 
2293   return normalize(rounding_mode, lost_fraction);
2294 }
2295 
2296 APFloat::opStatus
roundSignificandWithExponent(const integerPart * decSigParts,unsigned sigPartCount,int exp,roundingMode rounding_mode)2297 APFloat::roundSignificandWithExponent(const integerPart *decSigParts,
2298                                       unsigned sigPartCount, int exp,
2299                                       roundingMode rounding_mode)
2300 {
2301   unsigned int parts, pow5PartCount;
2302   fltSemantics calcSemantics = { 32767, -32767, 0, true };
2303   integerPart pow5Parts[maxPowerOfFiveParts];
2304   bool isNearest;
2305 
2306   isNearest = (rounding_mode == rmNearestTiesToEven ||
2307                rounding_mode == rmNearestTiesToAway);
2308 
2309   parts = partCountForBits(semantics->precision + 11);
2310 
2311   /* Calculate pow(5, abs(exp)).  */
2312   pow5PartCount = powerOf5(pow5Parts, exp >= 0 ? exp: -exp);
2313 
2314   for (;; parts *= 2) {
2315     opStatus sigStatus, powStatus;
2316     unsigned int excessPrecision, truncatedBits;
2317 
2318     calcSemantics.precision = parts * integerPartWidth - 1;
2319     excessPrecision = calcSemantics.precision - semantics->precision;
2320     truncatedBits = excessPrecision;
2321 
2322     APFloat decSig(calcSemantics, fcZero, sign);
2323     APFloat pow5(calcSemantics, fcZero, false);
2324 
2325     sigStatus = decSig.convertFromUnsignedParts(decSigParts, sigPartCount,
2326                                                 rmNearestTiesToEven);
2327     powStatus = pow5.convertFromUnsignedParts(pow5Parts, pow5PartCount,
2328                                               rmNearestTiesToEven);
2329     /* Add exp, as 10^n = 5^n * 2^n.  */
2330     decSig.exponent += exp;
2331 
2332     lostFraction calcLostFraction;
2333     integerPart HUerr, HUdistance;
2334     unsigned int powHUerr;
2335 
2336     if (exp >= 0) {
2337       /* multiplySignificand leaves the precision-th bit set to 1.  */
2338       calcLostFraction = decSig.multiplySignificand(pow5, NULL);
2339       powHUerr = powStatus != opOK;
2340     } else {
2341       calcLostFraction = decSig.divideSignificand(pow5);
2342       /* Denormal numbers have less precision.  */
2343       if (decSig.exponent < semantics->minExponent) {
2344         excessPrecision += (semantics->minExponent - decSig.exponent);
2345         truncatedBits = excessPrecision;
2346         if (excessPrecision > calcSemantics.precision)
2347           excessPrecision = calcSemantics.precision;
2348       }
2349       /* Extra half-ulp lost in reciprocal of exponent.  */
2350       powHUerr = (powStatus == opOK && calcLostFraction == lfExactlyZero) ? 0:2;
2351     }
2352 
2353     /* Both multiplySignificand and divideSignificand return the
2354        result with the integer bit set.  */
2355     assert(APInt::tcExtractBit
2356            (decSig.significandParts(), calcSemantics.precision - 1) == 1);
2357 
2358     HUerr = HUerrBound(calcLostFraction != lfExactlyZero, sigStatus != opOK,
2359                        powHUerr);
2360     HUdistance = 2 * ulpsFromBoundary(decSig.significandParts(),
2361                                       excessPrecision, isNearest);
2362 
2363     /* Are we guaranteed to round correctly if we truncate?  */
2364     if (HUdistance >= HUerr) {
2365       APInt::tcExtract(significandParts(), partCount(), decSig.significandParts(),
2366                        calcSemantics.precision - excessPrecision,
2367                        excessPrecision);
2368       /* Take the exponent of decSig.  If we tcExtract-ed less bits
2369          above we must adjust our exponent to compensate for the
2370          implicit right shift.  */
2371       exponent = (decSig.exponent + semantics->precision
2372                   - (calcSemantics.precision - excessPrecision));
2373       calcLostFraction = lostFractionThroughTruncation(decSig.significandParts(),
2374                                                        decSig.partCount(),
2375                                                        truncatedBits);
2376       return normalize(rounding_mode, calcLostFraction);
2377     }
2378   }
2379 }
2380 
2381 APFloat::opStatus
convertFromDecimalString(StringRef str,roundingMode rounding_mode)2382 APFloat::convertFromDecimalString(StringRef str, roundingMode rounding_mode)
2383 {
2384   decimalInfo D;
2385   opStatus fs;
2386 
2387   /* Scan the text.  */
2388   StringRef::iterator p = str.begin();
2389   interpretDecimal(p, str.end(), &D);
2390 
2391   /* Handle the quick cases.  First the case of no significant digits,
2392      i.e. zero, and then exponents that are obviously too large or too
2393      small.  Writing L for log 10 / log 2, a number d.ddddd*10^exp
2394      definitely overflows if
2395 
2396            (exp - 1) * L >= maxExponent
2397 
2398      and definitely underflows to zero where
2399 
2400            (exp + 1) * L <= minExponent - precision
2401 
2402      With integer arithmetic the tightest bounds for L are
2403 
2404            93/28 < L < 196/59            [ numerator <= 256 ]
2405            42039/12655 < L < 28738/8651  [ numerator <= 65536 ]
2406   */
2407 
2408   if (decDigitValue(*D.firstSigDigit) >= 10U) {
2409     category = fcZero;
2410     fs = opOK;
2411 
2412   /* Check whether the normalized exponent is high enough to overflow
2413      max during the log-rebasing in the max-exponent check below. */
2414   } else if (D.normalizedExponent - 1 > INT_MAX / 42039) {
2415     fs = handleOverflow(rounding_mode);
2416 
2417   /* If it wasn't, then it also wasn't high enough to overflow max
2418      during the log-rebasing in the min-exponent check.  Check that it
2419      won't overflow min in either check, then perform the min-exponent
2420      check. */
2421   } else if (D.normalizedExponent - 1 < INT_MIN / 42039 ||
2422              (D.normalizedExponent + 1) * 28738 <=
2423                8651 * (semantics->minExponent - (int) semantics->precision)) {
2424     /* Underflow to zero and round.  */
2425     zeroSignificand();
2426     fs = normalize(rounding_mode, lfLessThanHalf);
2427 
2428   /* We can finally safely perform the max-exponent check. */
2429   } else if ((D.normalizedExponent - 1) * 42039
2430              >= 12655 * semantics->maxExponent) {
2431     /* Overflow and round.  */
2432     fs = handleOverflow(rounding_mode);
2433   } else {
2434     integerPart *decSignificand;
2435     unsigned int partCount;
2436 
2437     /* A tight upper bound on number of bits required to hold an
2438        N-digit decimal integer is N * 196 / 59.  Allocate enough space
2439        to hold the full significand, and an extra part required by
2440        tcMultiplyPart.  */
2441     partCount = static_cast<unsigned int>(D.lastSigDigit - D.firstSigDigit) + 1;
2442     partCount = partCountForBits(1 + 196 * partCount / 59);
2443     decSignificand = new integerPart[partCount + 1];
2444     partCount = 0;
2445 
2446     /* Convert to binary efficiently - we do almost all multiplication
2447        in an integerPart.  When this would overflow do we do a single
2448        bignum multiplication, and then revert again to multiplication
2449        in an integerPart.  */
2450     do {
2451       integerPart decValue, val, multiplier;
2452 
2453       val = 0;
2454       multiplier = 1;
2455 
2456       do {
2457         if (*p == '.') {
2458           p++;
2459           if (p == str.end()) {
2460             break;
2461           }
2462         }
2463         decValue = decDigitValue(*p++);
2464         assert(decValue < 10U && "Invalid character in significand");
2465         multiplier *= 10;
2466         val = val * 10 + decValue;
2467         /* The maximum number that can be multiplied by ten with any
2468            digit added without overflowing an integerPart.  */
2469       } while (p <= D.lastSigDigit && multiplier <= (~ (integerPart) 0 - 9) / 10);
2470 
2471       /* Multiply out the current part.  */
2472       APInt::tcMultiplyPart(decSignificand, decSignificand, multiplier, val,
2473                             partCount, partCount + 1, false);
2474 
2475       /* If we used another part (likely but not guaranteed), increase
2476          the count.  */
2477       if (decSignificand[partCount])
2478         partCount++;
2479     } while (p <= D.lastSigDigit);
2480 
2481     category = fcNormal;
2482     fs = roundSignificandWithExponent(decSignificand, partCount,
2483                                       D.exponent, rounding_mode);
2484 
2485     delete [] decSignificand;
2486   }
2487 
2488   return fs;
2489 }
2490 
2491 APFloat::opStatus
convertFromString(StringRef str,roundingMode rounding_mode)2492 APFloat::convertFromString(StringRef str, roundingMode rounding_mode)
2493 {
2494   assertArithmeticOK(*semantics);
2495   assert(!str.empty() && "Invalid string length");
2496 
2497   /* Handle a leading minus sign.  */
2498   StringRef::iterator p = str.begin();
2499   size_t slen = str.size();
2500   sign = *p == '-' ? 1 : 0;
2501   if (*p == '-' || *p == '+') {
2502     p++;
2503     slen--;
2504     assert(slen && "String has no digits");
2505   }
2506 
2507   if (slen >= 2 && p[0] == '0' && (p[1] == 'x' || p[1] == 'X')) {
2508     assert(slen - 2 && "Invalid string");
2509     return convertFromHexadecimalString(StringRef(p + 2, slen - 2),
2510                                         rounding_mode);
2511   }
2512 
2513   return convertFromDecimalString(StringRef(p, slen), rounding_mode);
2514 }
2515 
2516 /* Write out a hexadecimal representation of the floating point value
2517    to DST, which must be of sufficient size, in the C99 form
2518    [-]0xh.hhhhp[+-]d.  Return the number of characters written,
2519    excluding the terminating NUL.
2520 
2521    If UPPERCASE, the output is in upper case, otherwise in lower case.
2522 
2523    HEXDIGITS digits appear altogether, rounding the value if
2524    necessary.  If HEXDIGITS is 0, the minimal precision to display the
2525    number precisely is used instead.  If nothing would appear after
2526    the decimal point it is suppressed.
2527 
2528    The decimal exponent is always printed and has at least one digit.
2529    Zero values display an exponent of zero.  Infinities and NaNs
2530    appear as "infinity" or "nan" respectively.
2531 
2532    The above rules are as specified by C99.  There is ambiguity about
2533    what the leading hexadecimal digit should be.  This implementation
2534    uses whatever is necessary so that the exponent is displayed as
2535    stored.  This implies the exponent will fall within the IEEE format
2536    range, and the leading hexadecimal digit will be 0 (for denormals),
2537    1 (normal numbers) or 2 (normal numbers rounded-away-from-zero with
2538    any other digits zero).
2539 */
2540 unsigned int
convertToHexString(char * dst,unsigned int hexDigits,bool upperCase,roundingMode rounding_mode) const2541 APFloat::convertToHexString(char *dst, unsigned int hexDigits,
2542                             bool upperCase, roundingMode rounding_mode) const
2543 {
2544   char *p;
2545 
2546   assertArithmeticOK(*semantics);
2547 
2548   p = dst;
2549   if (sign)
2550     *dst++ = '-';
2551 
2552   switch (category) {
2553   case fcInfinity:
2554     memcpy (dst, upperCase ? infinityU: infinityL, sizeof infinityU - 1);
2555     dst += sizeof infinityL - 1;
2556     break;
2557 
2558   case fcNaN:
2559     memcpy (dst, upperCase ? NaNU: NaNL, sizeof NaNU - 1);
2560     dst += sizeof NaNU - 1;
2561     break;
2562 
2563   case fcZero:
2564     *dst++ = '0';
2565     *dst++ = upperCase ? 'X': 'x';
2566     *dst++ = '0';
2567     if (hexDigits > 1) {
2568       *dst++ = '.';
2569       memset (dst, '0', hexDigits - 1);
2570       dst += hexDigits - 1;
2571     }
2572     *dst++ = upperCase ? 'P': 'p';
2573     *dst++ = '0';
2574     break;
2575 
2576   case fcNormal:
2577     dst = convertNormalToHexString (dst, hexDigits, upperCase, rounding_mode);
2578     break;
2579   }
2580 
2581   *dst = 0;
2582 
2583   return static_cast<unsigned int>(dst - p);
2584 }
2585 
2586 /* Does the hard work of outputting the correctly rounded hexadecimal
2587    form of a normal floating point number with the specified number of
2588    hexadecimal digits.  If HEXDIGITS is zero the minimum number of
2589    digits necessary to print the value precisely is output.  */
2590 char *
convertNormalToHexString(char * dst,unsigned int hexDigits,bool upperCase,roundingMode rounding_mode) const2591 APFloat::convertNormalToHexString(char *dst, unsigned int hexDigits,
2592                                   bool upperCase,
2593                                   roundingMode rounding_mode) const
2594 {
2595   unsigned int count, valueBits, shift, partsCount, outputDigits;
2596   const char *hexDigitChars;
2597   const integerPart *significand;
2598   char *p;
2599   bool roundUp;
2600 
2601   *dst++ = '0';
2602   *dst++ = upperCase ? 'X': 'x';
2603 
2604   roundUp = false;
2605   hexDigitChars = upperCase ? hexDigitsUpper: hexDigitsLower;
2606 
2607   significand = significandParts();
2608   partsCount = partCount();
2609 
2610   /* +3 because the first digit only uses the single integer bit, so
2611      we have 3 virtual zero most-significant-bits.  */
2612   valueBits = semantics->precision + 3;
2613   shift = integerPartWidth - valueBits % integerPartWidth;
2614 
2615   /* The natural number of digits required ignoring trailing
2616      insignificant zeroes.  */
2617   outputDigits = (valueBits - significandLSB () + 3) / 4;
2618 
2619   /* hexDigits of zero means use the required number for the
2620      precision.  Otherwise, see if we are truncating.  If we are,
2621      find out if we need to round away from zero.  */
2622   if (hexDigits) {
2623     if (hexDigits < outputDigits) {
2624       /* We are dropping non-zero bits, so need to check how to round.
2625          "bits" is the number of dropped bits.  */
2626       unsigned int bits;
2627       lostFraction fraction;
2628 
2629       bits = valueBits - hexDigits * 4;
2630       fraction = lostFractionThroughTruncation (significand, partsCount, bits);
2631       roundUp = roundAwayFromZero(rounding_mode, fraction, bits);
2632     }
2633     outputDigits = hexDigits;
2634   }
2635 
2636   /* Write the digits consecutively, and start writing in the location
2637      of the hexadecimal point.  We move the most significant digit
2638      left and add the hexadecimal point later.  */
2639   p = ++dst;
2640 
2641   count = (valueBits + integerPartWidth - 1) / integerPartWidth;
2642 
2643   while (outputDigits && count) {
2644     integerPart part;
2645 
2646     /* Put the most significant integerPartWidth bits in "part".  */
2647     if (--count == partsCount)
2648       part = 0;  /* An imaginary higher zero part.  */
2649     else
2650       part = significand[count] << shift;
2651 
2652     if (count && shift)
2653       part |= significand[count - 1] >> (integerPartWidth - shift);
2654 
2655     /* Convert as much of "part" to hexdigits as we can.  */
2656     unsigned int curDigits = integerPartWidth / 4;
2657 
2658     if (curDigits > outputDigits)
2659       curDigits = outputDigits;
2660     dst += partAsHex (dst, part, curDigits, hexDigitChars);
2661     outputDigits -= curDigits;
2662   }
2663 
2664   if (roundUp) {
2665     char *q = dst;
2666 
2667     /* Note that hexDigitChars has a trailing '0'.  */
2668     do {
2669       q--;
2670       *q = hexDigitChars[hexDigitValue (*q) + 1];
2671     } while (*q == '0');
2672     assert(q >= p);
2673   } else {
2674     /* Add trailing zeroes.  */
2675     memset (dst, '0', outputDigits);
2676     dst += outputDigits;
2677   }
2678 
2679   /* Move the most significant digit to before the point, and if there
2680      is something after the decimal point add it.  This must come
2681      after rounding above.  */
2682   p[-1] = p[0];
2683   if (dst -1 == p)
2684     dst--;
2685   else
2686     p[0] = '.';
2687 
2688   /* Finally output the exponent.  */
2689   *dst++ = upperCase ? 'P': 'p';
2690 
2691   return writeSignedDecimal (dst, exponent);
2692 }
2693 
2694 // For good performance it is desirable for different APFloats
2695 // to produce different integers.
2696 uint32_t
getHashValue() const2697 APFloat::getHashValue() const
2698 {
2699   if (category==fcZero) return sign<<8 | semantics->precision ;
2700   else if (category==fcInfinity) return sign<<9 | semantics->precision;
2701   else if (category==fcNaN) return 1<<10 | semantics->precision;
2702   else {
2703     uint32_t hash = sign<<11 | semantics->precision | exponent<<12;
2704     const integerPart* p = significandParts();
2705     for (int i=partCount(); i>0; i--, p++)
2706       hash ^= ((uint32_t)*p) ^ (uint32_t)((*p)>>32);
2707     return hash;
2708   }
2709 }
2710 
2711 // Conversion from APFloat to/from host float/double.  It may eventually be
2712 // possible to eliminate these and have everybody deal with APFloats, but that
2713 // will take a while.  This approach will not easily extend to long double.
2714 // Current implementation requires integerPartWidth==64, which is correct at
2715 // the moment but could be made more general.
2716 
2717 // Denormals have exponent minExponent in APFloat, but minExponent-1 in
2718 // the actual IEEE respresentations.  We compensate for that here.
2719 
2720 APInt
convertF80LongDoubleAPFloatToAPInt() const2721 APFloat::convertF80LongDoubleAPFloatToAPInt() const
2722 {
2723   assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended);
2724   assert(partCount()==2);
2725 
2726   uint64_t myexponent, mysignificand;
2727 
2728   if (category==fcNormal) {
2729     myexponent = exponent+16383; //bias
2730     mysignificand = significandParts()[0];
2731     if (myexponent==1 && !(mysignificand & 0x8000000000000000ULL))
2732       myexponent = 0;   // denormal
2733   } else if (category==fcZero) {
2734     myexponent = 0;
2735     mysignificand = 0;
2736   } else if (category==fcInfinity) {
2737     myexponent = 0x7fff;
2738     mysignificand = 0x8000000000000000ULL;
2739   } else {
2740     assert(category == fcNaN && "Unknown category");
2741     myexponent = 0x7fff;
2742     mysignificand = significandParts()[0];
2743   }
2744 
2745   uint64_t words[2];
2746   words[0] = mysignificand;
2747   words[1] =  ((uint64_t)(sign & 1) << 15) |
2748               (myexponent & 0x7fffLL);
2749   return APInt(80, words);
2750 }
2751 
2752 APInt
convertPPCDoubleDoubleAPFloatToAPInt() const2753 APFloat::convertPPCDoubleDoubleAPFloatToAPInt() const
2754 {
2755   assert(semantics == (const llvm::fltSemantics*)&PPCDoubleDouble);
2756   assert(partCount()==2);
2757 
2758   uint64_t myexponent, mysignificand, myexponent2, mysignificand2;
2759 
2760   if (category==fcNormal) {
2761     myexponent = exponent + 1023; //bias
2762     myexponent2 = exponent2 + 1023;
2763     mysignificand = significandParts()[0];
2764     mysignificand2 = significandParts()[1];
2765     if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
2766       myexponent = 0;   // denormal
2767     if (myexponent2==1 && !(mysignificand2 & 0x10000000000000LL))
2768       myexponent2 = 0;   // denormal
2769   } else if (category==fcZero) {
2770     myexponent = 0;
2771     mysignificand = 0;
2772     myexponent2 = 0;
2773     mysignificand2 = 0;
2774   } else if (category==fcInfinity) {
2775     myexponent = 0x7ff;
2776     myexponent2 = 0;
2777     mysignificand = 0;
2778     mysignificand2 = 0;
2779   } else {
2780     assert(category == fcNaN && "Unknown category");
2781     myexponent = 0x7ff;
2782     mysignificand = significandParts()[0];
2783     myexponent2 = exponent2;
2784     mysignificand2 = significandParts()[1];
2785   }
2786 
2787   uint64_t words[2];
2788   words[0] =  ((uint64_t)(sign & 1) << 63) |
2789               ((myexponent & 0x7ff) <<  52) |
2790               (mysignificand & 0xfffffffffffffLL);
2791   words[1] =  ((uint64_t)(sign2 & 1) << 63) |
2792               ((myexponent2 & 0x7ff) <<  52) |
2793               (mysignificand2 & 0xfffffffffffffLL);
2794   return APInt(128, words);
2795 }
2796 
2797 APInt
convertQuadrupleAPFloatToAPInt() const2798 APFloat::convertQuadrupleAPFloatToAPInt() const
2799 {
2800   assert(semantics == (const llvm::fltSemantics*)&IEEEquad);
2801   assert(partCount()==2);
2802 
2803   uint64_t myexponent, mysignificand, mysignificand2;
2804 
2805   if (category==fcNormal) {
2806     myexponent = exponent+16383; //bias
2807     mysignificand = significandParts()[0];
2808     mysignificand2 = significandParts()[1];
2809     if (myexponent==1 && !(mysignificand2 & 0x1000000000000LL))
2810       myexponent = 0;   // denormal
2811   } else if (category==fcZero) {
2812     myexponent = 0;
2813     mysignificand = mysignificand2 = 0;
2814   } else if (category==fcInfinity) {
2815     myexponent = 0x7fff;
2816     mysignificand = mysignificand2 = 0;
2817   } else {
2818     assert(category == fcNaN && "Unknown category!");
2819     myexponent = 0x7fff;
2820     mysignificand = significandParts()[0];
2821     mysignificand2 = significandParts()[1];
2822   }
2823 
2824   uint64_t words[2];
2825   words[0] = mysignificand;
2826   words[1] = ((uint64_t)(sign & 1) << 63) |
2827              ((myexponent & 0x7fff) << 48) |
2828              (mysignificand2 & 0xffffffffffffLL);
2829 
2830   return APInt(128, words);
2831 }
2832 
2833 APInt
convertDoubleAPFloatToAPInt() const2834 APFloat::convertDoubleAPFloatToAPInt() const
2835 {
2836   assert(semantics == (const llvm::fltSemantics*)&IEEEdouble);
2837   assert(partCount()==1);
2838 
2839   uint64_t myexponent, mysignificand;
2840 
2841   if (category==fcNormal) {
2842     myexponent = exponent+1023; //bias
2843     mysignificand = *significandParts();
2844     if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
2845       myexponent = 0;   // denormal
2846   } else if (category==fcZero) {
2847     myexponent = 0;
2848     mysignificand = 0;
2849   } else if (category==fcInfinity) {
2850     myexponent = 0x7ff;
2851     mysignificand = 0;
2852   } else {
2853     assert(category == fcNaN && "Unknown category!");
2854     myexponent = 0x7ff;
2855     mysignificand = *significandParts();
2856   }
2857 
2858   return APInt(64, ((((uint64_t)(sign & 1) << 63) |
2859                      ((myexponent & 0x7ff) <<  52) |
2860                      (mysignificand & 0xfffffffffffffLL))));
2861 }
2862 
2863 APInt
convertFloatAPFloatToAPInt() const2864 APFloat::convertFloatAPFloatToAPInt() const
2865 {
2866   assert(semantics == (const llvm::fltSemantics*)&IEEEsingle);
2867   assert(partCount()==1);
2868 
2869   uint32_t myexponent, mysignificand;
2870 
2871   if (category==fcNormal) {
2872     myexponent = exponent+127; //bias
2873     mysignificand = (uint32_t)*significandParts();
2874     if (myexponent == 1 && !(mysignificand & 0x800000))
2875       myexponent = 0;   // denormal
2876   } else if (category==fcZero) {
2877     myexponent = 0;
2878     mysignificand = 0;
2879   } else if (category==fcInfinity) {
2880     myexponent = 0xff;
2881     mysignificand = 0;
2882   } else {
2883     assert(category == fcNaN && "Unknown category!");
2884     myexponent = 0xff;
2885     mysignificand = (uint32_t)*significandParts();
2886   }
2887 
2888   return APInt(32, (((sign&1) << 31) | ((myexponent&0xff) << 23) |
2889                     (mysignificand & 0x7fffff)));
2890 }
2891 
2892 APInt
convertHalfAPFloatToAPInt() const2893 APFloat::convertHalfAPFloatToAPInt() const
2894 {
2895   assert(semantics == (const llvm::fltSemantics*)&IEEEhalf);
2896   assert(partCount()==1);
2897 
2898   uint32_t myexponent, mysignificand;
2899 
2900   if (category==fcNormal) {
2901     myexponent = exponent+15; //bias
2902     mysignificand = (uint32_t)*significandParts();
2903     if (myexponent == 1 && !(mysignificand & 0x400))
2904       myexponent = 0;   // denormal
2905   } else if (category==fcZero) {
2906     myexponent = 0;
2907     mysignificand = 0;
2908   } else if (category==fcInfinity) {
2909     myexponent = 0x1f;
2910     mysignificand = 0;
2911   } else {
2912     assert(category == fcNaN && "Unknown category!");
2913     myexponent = 0x1f;
2914     mysignificand = (uint32_t)*significandParts();
2915   }
2916 
2917   return APInt(16, (((sign&1) << 15) | ((myexponent&0x1f) << 10) |
2918                     (mysignificand & 0x3ff)));
2919 }
2920 
2921 // This function creates an APInt that is just a bit map of the floating
2922 // point constant as it would appear in memory.  It is not a conversion,
2923 // and treating the result as a normal integer is unlikely to be useful.
2924 
2925 APInt
bitcastToAPInt() const2926 APFloat::bitcastToAPInt() const
2927 {
2928   if (semantics == (const llvm::fltSemantics*)&IEEEhalf)
2929     return convertHalfAPFloatToAPInt();
2930 
2931   if (semantics == (const llvm::fltSemantics*)&IEEEsingle)
2932     return convertFloatAPFloatToAPInt();
2933 
2934   if (semantics == (const llvm::fltSemantics*)&IEEEdouble)
2935     return convertDoubleAPFloatToAPInt();
2936 
2937   if (semantics == (const llvm::fltSemantics*)&IEEEquad)
2938     return convertQuadrupleAPFloatToAPInt();
2939 
2940   if (semantics == (const llvm::fltSemantics*)&PPCDoubleDouble)
2941     return convertPPCDoubleDoubleAPFloatToAPInt();
2942 
2943   assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended &&
2944          "unknown format!");
2945   return convertF80LongDoubleAPFloatToAPInt();
2946 }
2947 
2948 float
convertToFloat() const2949 APFloat::convertToFloat() const
2950 {
2951   assert(semantics == (const llvm::fltSemantics*)&IEEEsingle &&
2952          "Float semantics are not IEEEsingle");
2953   APInt api = bitcastToAPInt();
2954   return api.bitsToFloat();
2955 }
2956 
2957 double
convertToDouble() const2958 APFloat::convertToDouble() const
2959 {
2960   assert(semantics == (const llvm::fltSemantics*)&IEEEdouble &&
2961          "Float semantics are not IEEEdouble");
2962   APInt api = bitcastToAPInt();
2963   return api.bitsToDouble();
2964 }
2965 
2966 /// Integer bit is explicit in this format.  Intel hardware (387 and later)
2967 /// does not support these bit patterns:
2968 ///  exponent = all 1's, integer bit 0, significand 0 ("pseudoinfinity")
2969 ///  exponent = all 1's, integer bit 0, significand nonzero ("pseudoNaN")
2970 ///  exponent = 0, integer bit 1 ("pseudodenormal")
2971 ///  exponent!=0 nor all 1's, integer bit 0 ("unnormal")
2972 /// At the moment, the first two are treated as NaNs, the second two as Normal.
2973 void
initFromF80LongDoubleAPInt(const APInt & api)2974 APFloat::initFromF80LongDoubleAPInt(const APInt &api)
2975 {
2976   assert(api.getBitWidth()==80);
2977   uint64_t i1 = api.getRawData()[0];
2978   uint64_t i2 = api.getRawData()[1];
2979   uint64_t myexponent = (i2 & 0x7fff);
2980   uint64_t mysignificand = i1;
2981 
2982   initialize(&APFloat::x87DoubleExtended);
2983   assert(partCount()==2);
2984 
2985   sign = static_cast<unsigned int>(i2>>15);
2986   if (myexponent==0 && mysignificand==0) {
2987     // exponent, significand meaningless
2988     category = fcZero;
2989   } else if (myexponent==0x7fff && mysignificand==0x8000000000000000ULL) {
2990     // exponent, significand meaningless
2991     category = fcInfinity;
2992   } else if (myexponent==0x7fff && mysignificand!=0x8000000000000000ULL) {
2993     // exponent meaningless
2994     category = fcNaN;
2995     significandParts()[0] = mysignificand;
2996     significandParts()[1] = 0;
2997   } else {
2998     category = fcNormal;
2999     exponent = myexponent - 16383;
3000     significandParts()[0] = mysignificand;
3001     significandParts()[1] = 0;
3002     if (myexponent==0)          // denormal
3003       exponent = -16382;
3004   }
3005 }
3006 
3007 void
initFromPPCDoubleDoubleAPInt(const APInt & api)3008 APFloat::initFromPPCDoubleDoubleAPInt(const APInt &api)
3009 {
3010   assert(api.getBitWidth()==128);
3011   uint64_t i1 = api.getRawData()[0];
3012   uint64_t i2 = api.getRawData()[1];
3013   uint64_t myexponent = (i1 >> 52) & 0x7ff;
3014   uint64_t mysignificand = i1 & 0xfffffffffffffLL;
3015   uint64_t myexponent2 = (i2 >> 52) & 0x7ff;
3016   uint64_t mysignificand2 = i2 & 0xfffffffffffffLL;
3017 
3018   initialize(&APFloat::PPCDoubleDouble);
3019   assert(partCount()==2);
3020 
3021   sign = static_cast<unsigned int>(i1>>63);
3022   sign2 = static_cast<unsigned int>(i2>>63);
3023   if (myexponent==0 && mysignificand==0) {
3024     // exponent, significand meaningless
3025     // exponent2 and significand2 are required to be 0; we don't check
3026     category = fcZero;
3027   } else if (myexponent==0x7ff && mysignificand==0) {
3028     // exponent, significand meaningless
3029     // exponent2 and significand2 are required to be 0; we don't check
3030     category = fcInfinity;
3031   } else if (myexponent==0x7ff && mysignificand!=0) {
3032     // exponent meaningless.  So is the whole second word, but keep it
3033     // for determinism.
3034     category = fcNaN;
3035     exponent2 = myexponent2;
3036     significandParts()[0] = mysignificand;
3037     significandParts()[1] = mysignificand2;
3038   } else {
3039     category = fcNormal;
3040     // Note there is no category2; the second word is treated as if it is
3041     // fcNormal, although it might be something else considered by itself.
3042     exponent = myexponent - 1023;
3043     exponent2 = myexponent2 - 1023;
3044     significandParts()[0] = mysignificand;
3045     significandParts()[1] = mysignificand2;
3046     if (myexponent==0)          // denormal
3047       exponent = -1022;
3048     else
3049       significandParts()[0] |= 0x10000000000000LL;  // integer bit
3050     if (myexponent2==0)
3051       exponent2 = -1022;
3052     else
3053       significandParts()[1] |= 0x10000000000000LL;  // integer bit
3054   }
3055 }
3056 
3057 void
initFromQuadrupleAPInt(const APInt & api)3058 APFloat::initFromQuadrupleAPInt(const APInt &api)
3059 {
3060   assert(api.getBitWidth()==128);
3061   uint64_t i1 = api.getRawData()[0];
3062   uint64_t i2 = api.getRawData()[1];
3063   uint64_t myexponent = (i2 >> 48) & 0x7fff;
3064   uint64_t mysignificand  = i1;
3065   uint64_t mysignificand2 = i2 & 0xffffffffffffLL;
3066 
3067   initialize(&APFloat::IEEEquad);
3068   assert(partCount()==2);
3069 
3070   sign = static_cast<unsigned int>(i2>>63);
3071   if (myexponent==0 &&
3072       (mysignificand==0 && mysignificand2==0)) {
3073     // exponent, significand meaningless
3074     category = fcZero;
3075   } else if (myexponent==0x7fff &&
3076              (mysignificand==0 && mysignificand2==0)) {
3077     // exponent, significand meaningless
3078     category = fcInfinity;
3079   } else if (myexponent==0x7fff &&
3080              (mysignificand!=0 || mysignificand2 !=0)) {
3081     // exponent meaningless
3082     category = fcNaN;
3083     significandParts()[0] = mysignificand;
3084     significandParts()[1] = mysignificand2;
3085   } else {
3086     category = fcNormal;
3087     exponent = myexponent - 16383;
3088     significandParts()[0] = mysignificand;
3089     significandParts()[1] = mysignificand2;
3090     if (myexponent==0)          // denormal
3091       exponent = -16382;
3092     else
3093       significandParts()[1] |= 0x1000000000000LL;  // integer bit
3094   }
3095 }
3096 
3097 void
initFromDoubleAPInt(const APInt & api)3098 APFloat::initFromDoubleAPInt(const APInt &api)
3099 {
3100   assert(api.getBitWidth()==64);
3101   uint64_t i = *api.getRawData();
3102   uint64_t myexponent = (i >> 52) & 0x7ff;
3103   uint64_t mysignificand = i & 0xfffffffffffffLL;
3104 
3105   initialize(&APFloat::IEEEdouble);
3106   assert(partCount()==1);
3107 
3108   sign = static_cast<unsigned int>(i>>63);
3109   if (myexponent==0 && mysignificand==0) {
3110     // exponent, significand meaningless
3111     category = fcZero;
3112   } else if (myexponent==0x7ff && mysignificand==0) {
3113     // exponent, significand meaningless
3114     category = fcInfinity;
3115   } else if (myexponent==0x7ff && mysignificand!=0) {
3116     // exponent meaningless
3117     category = fcNaN;
3118     *significandParts() = mysignificand;
3119   } else {
3120     category = fcNormal;
3121     exponent = myexponent - 1023;
3122     *significandParts() = mysignificand;
3123     if (myexponent==0)          // denormal
3124       exponent = -1022;
3125     else
3126       *significandParts() |= 0x10000000000000LL;  // integer bit
3127   }
3128 }
3129 
3130 void
initFromFloatAPInt(const APInt & api)3131 APFloat::initFromFloatAPInt(const APInt & api)
3132 {
3133   assert(api.getBitWidth()==32);
3134   uint32_t i = (uint32_t)*api.getRawData();
3135   uint32_t myexponent = (i >> 23) & 0xff;
3136   uint32_t mysignificand = i & 0x7fffff;
3137 
3138   initialize(&APFloat::IEEEsingle);
3139   assert(partCount()==1);
3140 
3141   sign = i >> 31;
3142   if (myexponent==0 && mysignificand==0) {
3143     // exponent, significand meaningless
3144     category = fcZero;
3145   } else if (myexponent==0xff && mysignificand==0) {
3146     // exponent, significand meaningless
3147     category = fcInfinity;
3148   } else if (myexponent==0xff && mysignificand!=0) {
3149     // sign, exponent, significand meaningless
3150     category = fcNaN;
3151     *significandParts() = mysignificand;
3152   } else {
3153     category = fcNormal;
3154     exponent = myexponent - 127;  //bias
3155     *significandParts() = mysignificand;
3156     if (myexponent==0)    // denormal
3157       exponent = -126;
3158     else
3159       *significandParts() |= 0x800000; // integer bit
3160   }
3161 }
3162 
3163 void
initFromHalfAPInt(const APInt & api)3164 APFloat::initFromHalfAPInt(const APInt & api)
3165 {
3166   assert(api.getBitWidth()==16);
3167   uint32_t i = (uint32_t)*api.getRawData();
3168   uint32_t myexponent = (i >> 10) & 0x1f;
3169   uint32_t mysignificand = i & 0x3ff;
3170 
3171   initialize(&APFloat::IEEEhalf);
3172   assert(partCount()==1);
3173 
3174   sign = i >> 15;
3175   if (myexponent==0 && mysignificand==0) {
3176     // exponent, significand meaningless
3177     category = fcZero;
3178   } else if (myexponent==0x1f && mysignificand==0) {
3179     // exponent, significand meaningless
3180     category = fcInfinity;
3181   } else if (myexponent==0x1f && mysignificand!=0) {
3182     // sign, exponent, significand meaningless
3183     category = fcNaN;
3184     *significandParts() = mysignificand;
3185   } else {
3186     category = fcNormal;
3187     exponent = myexponent - 15;  //bias
3188     *significandParts() = mysignificand;
3189     if (myexponent==0)    // denormal
3190       exponent = -14;
3191     else
3192       *significandParts() |= 0x400; // integer bit
3193   }
3194 }
3195 
3196 /// Treat api as containing the bits of a floating point number.  Currently
3197 /// we infer the floating point type from the size of the APInt.  The
3198 /// isIEEE argument distinguishes between PPC128 and IEEE128 (not meaningful
3199 /// when the size is anything else).
3200 void
initFromAPInt(const APInt & api,bool isIEEE)3201 APFloat::initFromAPInt(const APInt& api, bool isIEEE)
3202 {
3203   if (api.getBitWidth() == 16)
3204     return initFromHalfAPInt(api);
3205   else if (api.getBitWidth() == 32)
3206     return initFromFloatAPInt(api);
3207   else if (api.getBitWidth()==64)
3208     return initFromDoubleAPInt(api);
3209   else if (api.getBitWidth()==80)
3210     return initFromF80LongDoubleAPInt(api);
3211   else if (api.getBitWidth()==128)
3212     return (isIEEE ?
3213             initFromQuadrupleAPInt(api) : initFromPPCDoubleDoubleAPInt(api));
3214   else
3215     llvm_unreachable(0);
3216 }
3217 
3218 APFloat
getAllOnesValue(unsigned BitWidth,bool isIEEE)3219 APFloat::getAllOnesValue(unsigned BitWidth, bool isIEEE)
3220 {
3221   return APFloat(APInt::getAllOnesValue(BitWidth), isIEEE);
3222 }
3223 
getLargest(const fltSemantics & Sem,bool Negative)3224 APFloat APFloat::getLargest(const fltSemantics &Sem, bool Negative) {
3225   APFloat Val(Sem, fcNormal, Negative);
3226 
3227   // We want (in interchange format):
3228   //   sign = {Negative}
3229   //   exponent = 1..10
3230   //   significand = 1..1
3231 
3232   Val.exponent = Sem.maxExponent; // unbiased
3233 
3234   // 1-initialize all bits....
3235   Val.zeroSignificand();
3236   integerPart *significand = Val.significandParts();
3237   unsigned N = partCountForBits(Sem.precision);
3238   for (unsigned i = 0; i != N; ++i)
3239     significand[i] = ~((integerPart) 0);
3240 
3241   // ...and then clear the top bits for internal consistency.
3242   significand[N-1] &=
3243     (((integerPart) 1) << ((Sem.precision % integerPartWidth) - 1)) - 1;
3244 
3245   return Val;
3246 }
3247 
getSmallest(const fltSemantics & Sem,bool Negative)3248 APFloat APFloat::getSmallest(const fltSemantics &Sem, bool Negative) {
3249   APFloat Val(Sem, fcNormal, Negative);
3250 
3251   // We want (in interchange format):
3252   //   sign = {Negative}
3253   //   exponent = 0..0
3254   //   significand = 0..01
3255 
3256   Val.exponent = Sem.minExponent; // unbiased
3257   Val.zeroSignificand();
3258   Val.significandParts()[0] = 1;
3259   return Val;
3260 }
3261 
getSmallestNormalized(const fltSemantics & Sem,bool Negative)3262 APFloat APFloat::getSmallestNormalized(const fltSemantics &Sem, bool Negative) {
3263   APFloat Val(Sem, fcNormal, Negative);
3264 
3265   // We want (in interchange format):
3266   //   sign = {Negative}
3267   //   exponent = 0..0
3268   //   significand = 10..0
3269 
3270   Val.exponent = Sem.minExponent;
3271   Val.zeroSignificand();
3272   Val.significandParts()[partCountForBits(Sem.precision)-1] |=
3273     (((integerPart) 1) << ((Sem.precision % integerPartWidth) - 1));
3274 
3275   return Val;
3276 }
3277 
APFloat(const APInt & api,bool isIEEE)3278 APFloat::APFloat(const APInt& api, bool isIEEE) : exponent2(0), sign2(0) {
3279   initFromAPInt(api, isIEEE);
3280 }
3281 
APFloat(float f)3282 APFloat::APFloat(float f) : exponent2(0), sign2(0) {
3283   initFromAPInt(APInt::floatToBits(f));
3284 }
3285 
APFloat(double d)3286 APFloat::APFloat(double d) : exponent2(0), sign2(0) {
3287   initFromAPInt(APInt::doubleToBits(d));
3288 }
3289 
3290 namespace {
append(SmallVectorImpl<char> & Buffer,unsigned N,const char * Str)3291   static void append(SmallVectorImpl<char> &Buffer,
3292                      unsigned N, const char *Str) {
3293     unsigned Start = Buffer.size();
3294     Buffer.set_size(Start + N);
3295     memcpy(&Buffer[Start], Str, N);
3296   }
3297 
3298   template <unsigned N>
append(SmallVectorImpl<char> & Buffer,const char (& Str)[N])3299   void append(SmallVectorImpl<char> &Buffer, const char (&Str)[N]) {
3300     append(Buffer, N, Str);
3301   }
3302 
3303   /// Removes data from the given significand until it is no more
3304   /// precise than is required for the desired precision.
AdjustToPrecision(APInt & significand,int & exp,unsigned FormatPrecision)3305   void AdjustToPrecision(APInt &significand,
3306                          int &exp, unsigned FormatPrecision) {
3307     unsigned bits = significand.getActiveBits();
3308 
3309     // 196/59 is a very slight overestimate of lg_2(10).
3310     unsigned bitsRequired = (FormatPrecision * 196 + 58) / 59;
3311 
3312     if (bits <= bitsRequired) return;
3313 
3314     unsigned tensRemovable = (bits - bitsRequired) * 59 / 196;
3315     if (!tensRemovable) return;
3316 
3317     exp += tensRemovable;
3318 
3319     APInt divisor(significand.getBitWidth(), 1);
3320     APInt powten(significand.getBitWidth(), 10);
3321     while (true) {
3322       if (tensRemovable & 1)
3323         divisor *= powten;
3324       tensRemovable >>= 1;
3325       if (!tensRemovable) break;
3326       powten *= powten;
3327     }
3328 
3329     significand = significand.udiv(divisor);
3330 
3331     // Truncate the significand down to its active bit count, but
3332     // don't try to drop below 32.
3333     unsigned newPrecision = std::max(32U, significand.getActiveBits());
3334     significand = significand.trunc(newPrecision);
3335   }
3336 
3337 
AdjustToPrecision(SmallVectorImpl<char> & buffer,int & exp,unsigned FormatPrecision)3338   void AdjustToPrecision(SmallVectorImpl<char> &buffer,
3339                          int &exp, unsigned FormatPrecision) {
3340     unsigned N = buffer.size();
3341     if (N <= FormatPrecision) return;
3342 
3343     // The most significant figures are the last ones in the buffer.
3344     unsigned FirstSignificant = N - FormatPrecision;
3345 
3346     // Round.
3347     // FIXME: this probably shouldn't use 'round half up'.
3348 
3349     // Rounding down is just a truncation, except we also want to drop
3350     // trailing zeros from the new result.
3351     if (buffer[FirstSignificant - 1] < '5') {
3352       while (buffer[FirstSignificant] == '0')
3353         FirstSignificant++;
3354 
3355       exp += FirstSignificant;
3356       buffer.erase(&buffer[0], &buffer[FirstSignificant]);
3357       return;
3358     }
3359 
3360     // Rounding up requires a decimal add-with-carry.  If we continue
3361     // the carry, the newly-introduced zeros will just be truncated.
3362     for (unsigned I = FirstSignificant; I != N; ++I) {
3363       if (buffer[I] == '9') {
3364         FirstSignificant++;
3365       } else {
3366         buffer[I]++;
3367         break;
3368       }
3369     }
3370 
3371     // If we carried through, we have exactly one digit of precision.
3372     if (FirstSignificant == N) {
3373       exp += FirstSignificant;
3374       buffer.clear();
3375       buffer.push_back('1');
3376       return;
3377     }
3378 
3379     exp += FirstSignificant;
3380     buffer.erase(&buffer[0], &buffer[FirstSignificant]);
3381   }
3382 }
3383 
toString(SmallVectorImpl<char> & Str,unsigned FormatPrecision,unsigned FormatMaxPadding) const3384 void APFloat::toString(SmallVectorImpl<char> &Str,
3385                        unsigned FormatPrecision,
3386                        unsigned FormatMaxPadding) const {
3387   switch (category) {
3388   case fcInfinity:
3389     if (isNegative())
3390       return append(Str, "-Inf");
3391     else
3392       return append(Str, "+Inf");
3393 
3394   case fcNaN: return append(Str, "NaN");
3395 
3396   case fcZero:
3397     if (isNegative())
3398       Str.push_back('-');
3399 
3400     if (!FormatMaxPadding)
3401       append(Str, "0.0E+0");
3402     else
3403       Str.push_back('0');
3404     return;
3405 
3406   case fcNormal:
3407     break;
3408   }
3409 
3410   if (isNegative())
3411     Str.push_back('-');
3412 
3413   // Decompose the number into an APInt and an exponent.
3414   int exp = exponent - ((int) semantics->precision - 1);
3415   APInt significand(semantics->precision,
3416                     makeArrayRef(significandParts(),
3417                                  partCountForBits(semantics->precision)));
3418 
3419   // Set FormatPrecision if zero.  We want to do this before we
3420   // truncate trailing zeros, as those are part of the precision.
3421   if (!FormatPrecision) {
3422     // It's an interesting question whether to use the nominal
3423     // precision or the active precision here for denormals.
3424 
3425     // FormatPrecision = ceil(significandBits / lg_2(10))
3426     FormatPrecision = (semantics->precision * 59 + 195) / 196;
3427   }
3428 
3429   // Ignore trailing binary zeros.
3430   int trailingZeros = significand.countTrailingZeros();
3431   exp += trailingZeros;
3432   significand = significand.lshr(trailingZeros);
3433 
3434   // Change the exponent from 2^e to 10^e.
3435   if (exp == 0) {
3436     // Nothing to do.
3437   } else if (exp > 0) {
3438     // Just shift left.
3439     significand = significand.zext(semantics->precision + exp);
3440     significand <<= exp;
3441     exp = 0;
3442   } else { /* exp < 0 */
3443     int texp = -exp;
3444 
3445     // We transform this using the identity:
3446     //   (N)(2^-e) == (N)(5^e)(10^-e)
3447     // This means we have to multiply N (the significand) by 5^e.
3448     // To avoid overflow, we have to operate on numbers large
3449     // enough to store N * 5^e:
3450     //   log2(N * 5^e) == log2(N) + e * log2(5)
3451     //                 <= semantics->precision + e * 137 / 59
3452     //   (log_2(5) ~ 2.321928 < 2.322034 ~ 137/59)
3453 
3454     unsigned precision = semantics->precision + 137 * texp / 59;
3455 
3456     // Multiply significand by 5^e.
3457     //   N * 5^0101 == N * 5^(1*1) * 5^(0*2) * 5^(1*4) * 5^(0*8)
3458     significand = significand.zext(precision);
3459     APInt five_to_the_i(precision, 5);
3460     while (true) {
3461       if (texp & 1) significand *= five_to_the_i;
3462 
3463       texp >>= 1;
3464       if (!texp) break;
3465       five_to_the_i *= five_to_the_i;
3466     }
3467   }
3468 
3469   AdjustToPrecision(significand, exp, FormatPrecision);
3470 
3471   llvm::SmallVector<char, 256> buffer;
3472 
3473   // Fill the buffer.
3474   unsigned precision = significand.getBitWidth();
3475   APInt ten(precision, 10);
3476   APInt digit(precision, 0);
3477 
3478   bool inTrail = true;
3479   while (significand != 0) {
3480     // digit <- significand % 10
3481     // significand <- significand / 10
3482     APInt::udivrem(significand, ten, significand, digit);
3483 
3484     unsigned d = digit.getZExtValue();
3485 
3486     // Drop trailing zeros.
3487     if (inTrail && !d) exp++;
3488     else {
3489       buffer.push_back((char) ('0' + d));
3490       inTrail = false;
3491     }
3492   }
3493 
3494   assert(!buffer.empty() && "no characters in buffer!");
3495 
3496   // Drop down to FormatPrecision.
3497   // TODO: don't do more precise calculations above than are required.
3498   AdjustToPrecision(buffer, exp, FormatPrecision);
3499 
3500   unsigned NDigits = buffer.size();
3501 
3502   // Check whether we should use scientific notation.
3503   bool FormatScientific;
3504   if (!FormatMaxPadding)
3505     FormatScientific = true;
3506   else {
3507     if (exp >= 0) {
3508       // 765e3 --> 765000
3509       //              ^^^
3510       // But we shouldn't make the number look more precise than it is.
3511       FormatScientific = ((unsigned) exp > FormatMaxPadding ||
3512                           NDigits + (unsigned) exp > FormatPrecision);
3513     } else {
3514       // Power of the most significant digit.
3515       int MSD = exp + (int) (NDigits - 1);
3516       if (MSD >= 0) {
3517         // 765e-2 == 7.65
3518         FormatScientific = false;
3519       } else {
3520         // 765e-5 == 0.00765
3521         //           ^ ^^
3522         FormatScientific = ((unsigned) -MSD) > FormatMaxPadding;
3523       }
3524     }
3525   }
3526 
3527   // Scientific formatting is pretty straightforward.
3528   if (FormatScientific) {
3529     exp += (NDigits - 1);
3530 
3531     Str.push_back(buffer[NDigits-1]);
3532     Str.push_back('.');
3533     if (NDigits == 1)
3534       Str.push_back('0');
3535     else
3536       for (unsigned I = 1; I != NDigits; ++I)
3537         Str.push_back(buffer[NDigits-1-I]);
3538     Str.push_back('E');
3539 
3540     Str.push_back(exp >= 0 ? '+' : '-');
3541     if (exp < 0) exp = -exp;
3542     SmallVector<char, 6> expbuf;
3543     do {
3544       expbuf.push_back((char) ('0' + (exp % 10)));
3545       exp /= 10;
3546     } while (exp);
3547     for (unsigned I = 0, E = expbuf.size(); I != E; ++I)
3548       Str.push_back(expbuf[E-1-I]);
3549     return;
3550   }
3551 
3552   // Non-scientific, positive exponents.
3553   if (exp >= 0) {
3554     for (unsigned I = 0; I != NDigits; ++I)
3555       Str.push_back(buffer[NDigits-1-I]);
3556     for (unsigned I = 0; I != (unsigned) exp; ++I)
3557       Str.push_back('0');
3558     return;
3559   }
3560 
3561   // Non-scientific, negative exponents.
3562 
3563   // The number of digits to the left of the decimal point.
3564   int NWholeDigits = exp + (int) NDigits;
3565 
3566   unsigned I = 0;
3567   if (NWholeDigits > 0) {
3568     for (; I != (unsigned) NWholeDigits; ++I)
3569       Str.push_back(buffer[NDigits-I-1]);
3570     Str.push_back('.');
3571   } else {
3572     unsigned NZeros = 1 + (unsigned) -NWholeDigits;
3573 
3574     Str.push_back('0');
3575     Str.push_back('.');
3576     for (unsigned Z = 1; Z != NZeros; ++Z)
3577       Str.push_back('0');
3578   }
3579 
3580   for (; I != NDigits; ++I)
3581     Str.push_back(buffer[NDigits-I-1]);
3582 }
3583 
getExactInverse(APFloat * inv) const3584 bool APFloat::getExactInverse(APFloat *inv) const {
3585   // We can only guarantee the existence of an exact inverse for IEEE floats.
3586   if (semantics != &IEEEhalf && semantics != &IEEEsingle &&
3587       semantics != &IEEEdouble && semantics != &IEEEquad)
3588     return false;
3589 
3590   // Special floats and denormals have no exact inverse.
3591   if (category != fcNormal)
3592     return false;
3593 
3594   // Check that the number is a power of two by making sure that only the
3595   // integer bit is set in the significand.
3596   if (significandLSB() != semantics->precision - 1)
3597     return false;
3598 
3599   // Get the inverse.
3600   APFloat reciprocal(*semantics, 1ULL);
3601   if (reciprocal.divide(*this, rmNearestTiesToEven) != opOK)
3602     return false;
3603 
3604   // Avoid multiplication with a denormal, it is not safe on all platforms and
3605   // may be slower than a normal division.
3606   if (reciprocal.significandMSB() + 1 < reciprocal.semantics->precision)
3607     return false;
3608 
3609   assert(reciprocal.category == fcNormal &&
3610          reciprocal.significandLSB() == reciprocal.semantics->precision - 1);
3611 
3612   if (inv)
3613     *inv = reciprocal;
3614 
3615   return true;
3616 }
3617