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