Fix another roundToIntegral bug where very large values could become infinity. Probl...
[oota-llvm.git] / lib / Support / APFloat.cpp
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
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
94 decDigitValue(unsigned int c)
95 {
96   return c - '0';
97 }
98
99 static unsigned int
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
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
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
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
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
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
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
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
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
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
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
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
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
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 *
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 *
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
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
599 APFloat::freeSignificand()
600 {
601   if (partCount() > 1)
602     delete [] significand.parts;
603 }
604
605 void
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
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).  */
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
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 &
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
701 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
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
741 APFloat::APFloat(const fltSemantics &ourSemantics) : exponent2(0), sign2(0) {
742   assertArithmeticOK(ourSemantics);
743   initialize(&ourSemantics);
744   category = fcZero;
745   sign = false;
746 }
747
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
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
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
775 APFloat::APFloat(const APFloat &rhs) : exponent2(0), sign2(0) {
776   initialize(rhs.semantics);
777   assign(rhs);
778 }
779
780 APFloat::~APFloat()
781 {
782   freeSignificand();
783 }
784
785 // Profile - This method 'profiles' an APFloat for use with FoldingSet.
786 void APFloat::Profile(FoldingSetNodeID& ID) const {
787   ID.Add(bitcastToAPInt());
788 }
789
790 unsigned int
791 APFloat::partCount() const
792 {
793   return partCountForBits(semantics->precision + 1);
794 }
795
796 unsigned int
797 APFloat::semanticsPrecision(const fltSemantics &semantics)
798 {
799   return semantics.precision;
800 }
801
802 const integerPart *
803 APFloat::significandParts() const
804 {
805   return const_cast<APFloat *>(this)->significandParts();
806 }
807
808 integerPart *
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
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
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
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
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
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
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
1050 APFloat::significandMSB() const
1051 {
1052   return APInt::tcMSB(significandParts(), partCount());
1053 }
1054
1055 unsigned int
1056 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
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
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
1090 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
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
1143 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
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
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
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
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
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
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
1535 APFloat::changeSign()
1536 {
1537   /* Look mummy, this one's easy.  */
1538   sign = !sign;
1539 }
1540
1541 void
1542 APFloat::clearSign()
1543 {
1544   /* So is this one. */
1545   sign = 0;
1546 }
1547
1548 void
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
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
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
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
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
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
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
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
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 /* Rounding-mode corrrect round to integral value.  */
1769 APFloat::opStatus APFloat::roundToIntegral(roundingMode rounding_mode) {
1770   opStatus fs;
1771   assertArithmeticOK(*semantics);
1772
1773   // If the exponent is large enough, we know that this value is already
1774   // integral, and the arithmetic below would potentially cause it to saturate
1775   // to +/-Inf.  Bail out early instead.
1776   if (exponent+1 >= (int)semanticsPrecision(*semantics))
1777     return opOK;
1778
1779   // The algorithm here is quite simple: we add 2^(p-1), where p is the
1780   // precision of our format, and then subtract it back off again.  The choice
1781   // of rounding modes for the addition/subtraction determines the rounding mode
1782   // for our integral rounding as well.
1783   // NOTE: When the input value is negative, we do subtraction followed by
1784   // addition instead.
1785   APInt IntegerConstant(NextPowerOf2(semanticsPrecision(*semantics)), 1);
1786   IntegerConstant <<= semanticsPrecision(*semantics)-1;
1787   APFloat MagicConstant(*semantics);
1788   fs = MagicConstant.convertFromAPInt(IntegerConstant, false,
1789                                       rmNearestTiesToEven);
1790   MagicConstant.copySign(*this);
1791
1792   if (fs != opOK)
1793     return fs;
1794
1795   // Preserve the input sign so that we can handle 0.0/-0.0 cases correctly.
1796   bool inputSign = isNegative();
1797
1798   fs = add(MagicConstant, rounding_mode);
1799   if (fs != opOK && fs != opInexact)
1800     return fs;
1801
1802   fs = subtract(MagicConstant, rounding_mode);
1803
1804   // Restore the input sign.
1805   if (inputSign != isNegative())
1806     changeSign();
1807
1808   return fs;
1809 }
1810
1811
1812 /* Comparison requires normalized numbers.  */
1813 APFloat::cmpResult
1814 APFloat::compare(const APFloat &rhs) const
1815 {
1816   cmpResult result;
1817
1818   assertArithmeticOK(*semantics);
1819   assert(semantics == rhs.semantics);
1820
1821   switch (convolve(category, rhs.category)) {
1822   default:
1823     llvm_unreachable(0);
1824
1825   case convolve(fcNaN, fcZero):
1826   case convolve(fcNaN, fcNormal):
1827   case convolve(fcNaN, fcInfinity):
1828   case convolve(fcNaN, fcNaN):
1829   case convolve(fcZero, fcNaN):
1830   case convolve(fcNormal, fcNaN):
1831   case convolve(fcInfinity, fcNaN):
1832     return cmpUnordered;
1833
1834   case convolve(fcInfinity, fcNormal):
1835   case convolve(fcInfinity, fcZero):
1836   case convolve(fcNormal, fcZero):
1837     if (sign)
1838       return cmpLessThan;
1839     else
1840       return cmpGreaterThan;
1841
1842   case convolve(fcNormal, fcInfinity):
1843   case convolve(fcZero, fcInfinity):
1844   case convolve(fcZero, fcNormal):
1845     if (rhs.sign)
1846       return cmpGreaterThan;
1847     else
1848       return cmpLessThan;
1849
1850   case convolve(fcInfinity, fcInfinity):
1851     if (sign == rhs.sign)
1852       return cmpEqual;
1853     else if (sign)
1854       return cmpLessThan;
1855     else
1856       return cmpGreaterThan;
1857
1858   case convolve(fcZero, fcZero):
1859     return cmpEqual;
1860
1861   case convolve(fcNormal, fcNormal):
1862     break;
1863   }
1864
1865   /* Two normal numbers.  Do they have the same sign?  */
1866   if (sign != rhs.sign) {
1867     if (sign)
1868       result = cmpLessThan;
1869     else
1870       result = cmpGreaterThan;
1871   } else {
1872     /* Compare absolute values; invert result if negative.  */
1873     result = compareAbsoluteValue(rhs);
1874
1875     if (sign) {
1876       if (result == cmpLessThan)
1877         result = cmpGreaterThan;
1878       else if (result == cmpGreaterThan)
1879         result = cmpLessThan;
1880     }
1881   }
1882
1883   return result;
1884 }
1885
1886 /// APFloat::convert - convert a value of one floating point type to another.
1887 /// The return value corresponds to the IEEE754 exceptions.  *losesInfo
1888 /// records whether the transformation lost information, i.e. whether
1889 /// converting the result back to the original type will produce the
1890 /// original value (this is almost the same as return value==fsOK, but there
1891 /// are edge cases where this is not so).
1892
1893 APFloat::opStatus
1894 APFloat::convert(const fltSemantics &toSemantics,
1895                  roundingMode rounding_mode, bool *losesInfo)
1896 {
1897   lostFraction lostFraction;
1898   unsigned int newPartCount, oldPartCount;
1899   opStatus fs;
1900   int shift;
1901   const fltSemantics &fromSemantics = *semantics;
1902
1903   assertArithmeticOK(fromSemantics);
1904   assertArithmeticOK(toSemantics);
1905   lostFraction = lfExactlyZero;
1906   newPartCount = partCountForBits(toSemantics.precision + 1);
1907   oldPartCount = partCount();
1908   shift = toSemantics.precision - fromSemantics.precision;
1909
1910   bool X86SpecialNan = false;
1911   if (&fromSemantics == &APFloat::x87DoubleExtended &&
1912       &toSemantics != &APFloat::x87DoubleExtended && category == fcNaN &&
1913       (!(*significandParts() & 0x8000000000000000ULL) ||
1914        !(*significandParts() & 0x4000000000000000ULL))) {
1915     // x86 has some unusual NaNs which cannot be represented in any other
1916     // format; note them here.
1917     X86SpecialNan = true;
1918   }
1919
1920   // If this is a truncation, perform the shift before we narrow the storage.
1921   if (shift < 0 && (category==fcNormal || category==fcNaN))
1922     lostFraction = shiftRight(significandParts(), oldPartCount, -shift);
1923
1924   // Fix the storage so it can hold to new value.
1925   if (newPartCount > oldPartCount) {
1926     // The new type requires more storage; make it available.
1927     integerPart *newParts;
1928     newParts = new integerPart[newPartCount];
1929     APInt::tcSet(newParts, 0, newPartCount);
1930     if (category==fcNormal || category==fcNaN)
1931       APInt::tcAssign(newParts, significandParts(), oldPartCount);
1932     freeSignificand();
1933     significand.parts = newParts;
1934   } else if (newPartCount == 1 && oldPartCount != 1) {
1935     // Switch to built-in storage for a single part.
1936     integerPart newPart = 0;
1937     if (category==fcNormal || category==fcNaN)
1938       newPart = significandParts()[0];
1939     freeSignificand();
1940     significand.part = newPart;
1941   }
1942
1943   // Now that we have the right storage, switch the semantics.
1944   semantics = &toSemantics;
1945
1946   // If this is an extension, perform the shift now that the storage is
1947   // available.
1948   if (shift > 0 && (category==fcNormal || category==fcNaN))
1949     APInt::tcShiftLeft(significandParts(), newPartCount, shift);
1950
1951   if (category == fcNormal) {
1952     fs = normalize(rounding_mode, lostFraction);
1953     *losesInfo = (fs != opOK);
1954   } else if (category == fcNaN) {
1955     *losesInfo = lostFraction != lfExactlyZero || X86SpecialNan;
1956     // gcc forces the Quiet bit on, which means (float)(double)(float_sNan)
1957     // does not give you back the same bits.  This is dubious, and we
1958     // don't currently do it.  You're really supposed to get
1959     // an invalid operation signal at runtime, but nobody does that.
1960     fs = opOK;
1961   } else {
1962     *losesInfo = false;
1963     fs = opOK;
1964   }
1965
1966   return fs;
1967 }
1968
1969 /* Convert a floating point number to an integer according to the
1970    rounding mode.  If the rounded integer value is out of range this
1971    returns an invalid operation exception and the contents of the
1972    destination parts are unspecified.  If the rounded value is in
1973    range but the floating point number is not the exact integer, the C
1974    standard doesn't require an inexact exception to be raised.  IEEE
1975    854 does require it so we do that.
1976
1977    Note that for conversions to integer type the C standard requires
1978    round-to-zero to always be used.  */
1979 APFloat::opStatus
1980 APFloat::convertToSignExtendedInteger(integerPart *parts, unsigned int width,
1981                                       bool isSigned,
1982                                       roundingMode rounding_mode,
1983                                       bool *isExact) const
1984 {
1985   lostFraction lost_fraction;
1986   const integerPart *src;
1987   unsigned int dstPartsCount, truncatedBits;
1988
1989   assertArithmeticOK(*semantics);
1990
1991   *isExact = false;
1992
1993   /* Handle the three special cases first.  */
1994   if (category == fcInfinity || category == fcNaN)
1995     return opInvalidOp;
1996
1997   dstPartsCount = partCountForBits(width);
1998
1999   if (category == fcZero) {
2000     APInt::tcSet(parts, 0, dstPartsCount);
2001     // Negative zero can't be represented as an int.
2002     *isExact = !sign;
2003     return opOK;
2004   }
2005
2006   src = significandParts();
2007
2008   /* Step 1: place our absolute value, with any fraction truncated, in
2009      the destination.  */
2010   if (exponent < 0) {
2011     /* Our absolute value is less than one; truncate everything.  */
2012     APInt::tcSet(parts, 0, dstPartsCount);
2013     /* For exponent -1 the integer bit represents .5, look at that.
2014        For smaller exponents leftmost truncated bit is 0. */
2015     truncatedBits = semantics->precision -1U - exponent;
2016   } else {
2017     /* We want the most significant (exponent + 1) bits; the rest are
2018        truncated.  */
2019     unsigned int bits = exponent + 1U;
2020
2021     /* Hopelessly large in magnitude?  */
2022     if (bits > width)
2023       return opInvalidOp;
2024
2025     if (bits < semantics->precision) {
2026       /* We truncate (semantics->precision - bits) bits.  */
2027       truncatedBits = semantics->precision - bits;
2028       APInt::tcExtract(parts, dstPartsCount, src, bits, truncatedBits);
2029     } else {
2030       /* We want at least as many bits as are available.  */
2031       APInt::tcExtract(parts, dstPartsCount, src, semantics->precision, 0);
2032       APInt::tcShiftLeft(parts, dstPartsCount, bits - semantics->precision);
2033       truncatedBits = 0;
2034     }
2035   }
2036
2037   /* Step 2: work out any lost fraction, and increment the absolute
2038      value if we would round away from zero.  */
2039   if (truncatedBits) {
2040     lost_fraction = lostFractionThroughTruncation(src, partCount(),
2041                                                   truncatedBits);
2042     if (lost_fraction != lfExactlyZero &&
2043         roundAwayFromZero(rounding_mode, lost_fraction, truncatedBits)) {
2044       if (APInt::tcIncrement(parts, dstPartsCount))
2045         return opInvalidOp;     /* Overflow.  */
2046     }
2047   } else {
2048     lost_fraction = lfExactlyZero;
2049   }
2050
2051   /* Step 3: check if we fit in the destination.  */
2052   unsigned int omsb = APInt::tcMSB(parts, dstPartsCount) + 1;
2053
2054   if (sign) {
2055     if (!isSigned) {
2056       /* Negative numbers cannot be represented as unsigned.  */
2057       if (omsb != 0)
2058         return opInvalidOp;
2059     } else {
2060       /* It takes omsb bits to represent the unsigned integer value.
2061          We lose a bit for the sign, but care is needed as the
2062          maximally negative integer is a special case.  */
2063       if (omsb == width && APInt::tcLSB(parts, dstPartsCount) + 1 != omsb)
2064         return opInvalidOp;
2065
2066       /* This case can happen because of rounding.  */
2067       if (omsb > width)
2068         return opInvalidOp;
2069     }
2070
2071     APInt::tcNegate (parts, dstPartsCount);
2072   } else {
2073     if (omsb >= width + !isSigned)
2074       return opInvalidOp;
2075   }
2076
2077   if (lost_fraction == lfExactlyZero) {
2078     *isExact = true;
2079     return opOK;
2080   } else
2081     return opInexact;
2082 }
2083
2084 /* Same as convertToSignExtendedInteger, except we provide
2085    deterministic values in case of an invalid operation exception,
2086    namely zero for NaNs and the minimal or maximal value respectively
2087    for underflow or overflow.
2088    The *isExact output tells whether the result is exact, in the sense
2089    that converting it back to the original floating point type produces
2090    the original value.  This is almost equivalent to result==opOK,
2091    except for negative zeroes.
2092 */
2093 APFloat::opStatus
2094 APFloat::convertToInteger(integerPart *parts, unsigned int width,
2095                           bool isSigned,
2096                           roundingMode rounding_mode, bool *isExact) const
2097 {
2098   opStatus fs;
2099
2100   fs = convertToSignExtendedInteger(parts, width, isSigned, rounding_mode,
2101                                     isExact);
2102
2103   if (fs == opInvalidOp) {
2104     unsigned int bits, dstPartsCount;
2105
2106     dstPartsCount = partCountForBits(width);
2107
2108     if (category == fcNaN)
2109       bits = 0;
2110     else if (sign)
2111       bits = isSigned;
2112     else
2113       bits = width - isSigned;
2114
2115     APInt::tcSetLeastSignificantBits(parts, dstPartsCount, bits);
2116     if (sign && isSigned)
2117       APInt::tcShiftLeft(parts, dstPartsCount, width - 1);
2118   }
2119
2120   return fs;
2121 }
2122
2123 /* Same as convertToInteger(integerPart*, ...), except the result is returned in
2124    an APSInt, whose initial bit-width and signed-ness are used to determine the
2125    precision of the conversion.
2126  */
2127 APFloat::opStatus
2128 APFloat::convertToInteger(APSInt &result,
2129                           roundingMode rounding_mode, bool *isExact) const
2130 {
2131   unsigned bitWidth = result.getBitWidth();
2132   SmallVector<uint64_t, 4> parts(result.getNumWords());
2133   opStatus status = convertToInteger(
2134     parts.data(), bitWidth, result.isSigned(), rounding_mode, isExact);
2135   // Keeps the original signed-ness.
2136   result = APInt(bitWidth, parts);
2137   return status;
2138 }
2139
2140 /* Convert an unsigned integer SRC to a floating point number,
2141    rounding according to ROUNDING_MODE.  The sign of the floating
2142    point number is not modified.  */
2143 APFloat::opStatus
2144 APFloat::convertFromUnsignedParts(const integerPart *src,
2145                                   unsigned int srcCount,
2146                                   roundingMode rounding_mode)
2147 {
2148   unsigned int omsb, precision, dstCount;
2149   integerPart *dst;
2150   lostFraction lost_fraction;
2151
2152   assertArithmeticOK(*semantics);
2153   category = fcNormal;
2154   omsb = APInt::tcMSB(src, srcCount) + 1;
2155   dst = significandParts();
2156   dstCount = partCount();
2157   precision = semantics->precision;
2158
2159   /* We want the most significant PRECISION bits of SRC.  There may not
2160      be that many; extract what we can.  */
2161   if (precision <= omsb) {
2162     exponent = omsb - 1;
2163     lost_fraction = lostFractionThroughTruncation(src, srcCount,
2164                                                   omsb - precision);
2165     APInt::tcExtract(dst, dstCount, src, precision, omsb - precision);
2166   } else {
2167     exponent = precision - 1;
2168     lost_fraction = lfExactlyZero;
2169     APInt::tcExtract(dst, dstCount, src, omsb, 0);
2170   }
2171
2172   return normalize(rounding_mode, lost_fraction);
2173 }
2174
2175 APFloat::opStatus
2176 APFloat::convertFromAPInt(const APInt &Val,
2177                           bool isSigned,
2178                           roundingMode rounding_mode)
2179 {
2180   unsigned int partCount = Val.getNumWords();
2181   APInt api = Val;
2182
2183   sign = false;
2184   if (isSigned && api.isNegative()) {
2185     sign = true;
2186     api = -api;
2187   }
2188
2189   return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode);
2190 }
2191
2192 /* Convert a two's complement integer SRC to a floating point number,
2193    rounding according to ROUNDING_MODE.  ISSIGNED is true if the
2194    integer is signed, in which case it must be sign-extended.  */
2195 APFloat::opStatus
2196 APFloat::convertFromSignExtendedInteger(const integerPart *src,
2197                                         unsigned int srcCount,
2198                                         bool isSigned,
2199                                         roundingMode rounding_mode)
2200 {
2201   opStatus status;
2202
2203   assertArithmeticOK(*semantics);
2204   if (isSigned &&
2205       APInt::tcExtractBit(src, srcCount * integerPartWidth - 1)) {
2206     integerPart *copy;
2207
2208     /* If we're signed and negative negate a copy.  */
2209     sign = true;
2210     copy = new integerPart[srcCount];
2211     APInt::tcAssign(copy, src, srcCount);
2212     APInt::tcNegate(copy, srcCount);
2213     status = convertFromUnsignedParts(copy, srcCount, rounding_mode);
2214     delete [] copy;
2215   } else {
2216     sign = false;
2217     status = convertFromUnsignedParts(src, srcCount, rounding_mode);
2218   }
2219
2220   return status;
2221 }
2222
2223 /* FIXME: should this just take a const APInt reference?  */
2224 APFloat::opStatus
2225 APFloat::convertFromZeroExtendedInteger(const integerPart *parts,
2226                                         unsigned int width, bool isSigned,
2227                                         roundingMode rounding_mode)
2228 {
2229   unsigned int partCount = partCountForBits(width);
2230   APInt api = APInt(width, makeArrayRef(parts, partCount));
2231
2232   sign = false;
2233   if (isSigned && APInt::tcExtractBit(parts, width - 1)) {
2234     sign = true;
2235     api = -api;
2236   }
2237
2238   return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode);
2239 }
2240
2241 APFloat::opStatus
2242 APFloat::convertFromHexadecimalString(StringRef s, roundingMode rounding_mode)
2243 {
2244   lostFraction lost_fraction = lfExactlyZero;
2245   integerPart *significand;
2246   unsigned int bitPos, partsCount;
2247   StringRef::iterator dot, firstSignificantDigit;
2248
2249   zeroSignificand();
2250   exponent = 0;
2251   category = fcNormal;
2252
2253   significand = significandParts();
2254   partsCount = partCount();
2255   bitPos = partsCount * integerPartWidth;
2256
2257   /* Skip leading zeroes and any (hexa)decimal point.  */
2258   StringRef::iterator begin = s.begin();
2259   StringRef::iterator end = s.end();
2260   StringRef::iterator p = skipLeadingZeroesAndAnyDot(begin, end, &dot);
2261   firstSignificantDigit = p;
2262
2263   for (; p != end;) {
2264     integerPart hex_value;
2265
2266     if (*p == '.') {
2267       assert(dot == end && "String contains multiple dots");
2268       dot = p++;
2269       if (p == end) {
2270         break;
2271       }
2272     }
2273
2274     hex_value = hexDigitValue(*p);
2275     if (hex_value == -1U) {
2276       break;
2277     }
2278
2279     p++;
2280
2281     if (p == end) {
2282       break;
2283     } else {
2284       /* Store the number whilst 4-bit nibbles remain.  */
2285       if (bitPos) {
2286         bitPos -= 4;
2287         hex_value <<= bitPos % integerPartWidth;
2288         significand[bitPos / integerPartWidth] |= hex_value;
2289       } else {
2290         lost_fraction = trailingHexadecimalFraction(p, end, hex_value);
2291         while (p != end && hexDigitValue(*p) != -1U)
2292           p++;
2293         break;
2294       }
2295     }
2296   }
2297
2298   /* Hex floats require an exponent but not a hexadecimal point.  */
2299   assert(p != end && "Hex strings require an exponent");
2300   assert((*p == 'p' || *p == 'P') && "Invalid character in significand");
2301   assert(p != begin && "Significand has no digits");
2302   assert((dot == end || p - begin != 1) && "Significand has no digits");
2303
2304   /* Ignore the exponent if we are zero.  */
2305   if (p != firstSignificantDigit) {
2306     int expAdjustment;
2307
2308     /* Implicit hexadecimal point?  */
2309     if (dot == end)
2310       dot = p;
2311
2312     /* Calculate the exponent adjustment implicit in the number of
2313        significant digits.  */
2314     expAdjustment = static_cast<int>(dot - firstSignificantDigit);
2315     if (expAdjustment < 0)
2316       expAdjustment++;
2317     expAdjustment = expAdjustment * 4 - 1;
2318
2319     /* Adjust for writing the significand starting at the most
2320        significant nibble.  */
2321     expAdjustment += semantics->precision;
2322     expAdjustment -= partsCount * integerPartWidth;
2323
2324     /* Adjust for the given exponent.  */
2325     exponent = totalExponent(p + 1, end, expAdjustment);
2326   }
2327
2328   return normalize(rounding_mode, lost_fraction);
2329 }
2330
2331 APFloat::opStatus
2332 APFloat::roundSignificandWithExponent(const integerPart *decSigParts,
2333                                       unsigned sigPartCount, int exp,
2334                                       roundingMode rounding_mode)
2335 {
2336   unsigned int parts, pow5PartCount;
2337   fltSemantics calcSemantics = { 32767, -32767, 0, true };
2338   integerPart pow5Parts[maxPowerOfFiveParts];
2339   bool isNearest;
2340
2341   isNearest = (rounding_mode == rmNearestTiesToEven ||
2342                rounding_mode == rmNearestTiesToAway);
2343
2344   parts = partCountForBits(semantics->precision + 11);
2345
2346   /* Calculate pow(5, abs(exp)).  */
2347   pow5PartCount = powerOf5(pow5Parts, exp >= 0 ? exp: -exp);
2348
2349   for (;; parts *= 2) {
2350     opStatus sigStatus, powStatus;
2351     unsigned int excessPrecision, truncatedBits;
2352
2353     calcSemantics.precision = parts * integerPartWidth - 1;
2354     excessPrecision = calcSemantics.precision - semantics->precision;
2355     truncatedBits = excessPrecision;
2356
2357     APFloat decSig(calcSemantics, fcZero, sign);
2358     APFloat pow5(calcSemantics, fcZero, false);
2359
2360     sigStatus = decSig.convertFromUnsignedParts(decSigParts, sigPartCount,
2361                                                 rmNearestTiesToEven);
2362     powStatus = pow5.convertFromUnsignedParts(pow5Parts, pow5PartCount,
2363                                               rmNearestTiesToEven);
2364     /* Add exp, as 10^n = 5^n * 2^n.  */
2365     decSig.exponent += exp;
2366
2367     lostFraction calcLostFraction;
2368     integerPart HUerr, HUdistance;
2369     unsigned int powHUerr;
2370
2371     if (exp >= 0) {
2372       /* multiplySignificand leaves the precision-th bit set to 1.  */
2373       calcLostFraction = decSig.multiplySignificand(pow5, NULL);
2374       powHUerr = powStatus != opOK;
2375     } else {
2376       calcLostFraction = decSig.divideSignificand(pow5);
2377       /* Denormal numbers have less precision.  */
2378       if (decSig.exponent < semantics->minExponent) {
2379         excessPrecision += (semantics->minExponent - decSig.exponent);
2380         truncatedBits = excessPrecision;
2381         if (excessPrecision > calcSemantics.precision)
2382           excessPrecision = calcSemantics.precision;
2383       }
2384       /* Extra half-ulp lost in reciprocal of exponent.  */
2385       powHUerr = (powStatus == opOK && calcLostFraction == lfExactlyZero) ? 0:2;
2386     }
2387
2388     /* Both multiplySignificand and divideSignificand return the
2389        result with the integer bit set.  */
2390     assert(APInt::tcExtractBit
2391            (decSig.significandParts(), calcSemantics.precision - 1) == 1);
2392
2393     HUerr = HUerrBound(calcLostFraction != lfExactlyZero, sigStatus != opOK,
2394                        powHUerr);
2395     HUdistance = 2 * ulpsFromBoundary(decSig.significandParts(),
2396                                       excessPrecision, isNearest);
2397
2398     /* Are we guaranteed to round correctly if we truncate?  */
2399     if (HUdistance >= HUerr) {
2400       APInt::tcExtract(significandParts(), partCount(), decSig.significandParts(),
2401                        calcSemantics.precision - excessPrecision,
2402                        excessPrecision);
2403       /* Take the exponent of decSig.  If we tcExtract-ed less bits
2404          above we must adjust our exponent to compensate for the
2405          implicit right shift.  */
2406       exponent = (decSig.exponent + semantics->precision
2407                   - (calcSemantics.precision - excessPrecision));
2408       calcLostFraction = lostFractionThroughTruncation(decSig.significandParts(),
2409                                                        decSig.partCount(),
2410                                                        truncatedBits);
2411       return normalize(rounding_mode, calcLostFraction);
2412     }
2413   }
2414 }
2415
2416 APFloat::opStatus
2417 APFloat::convertFromDecimalString(StringRef str, roundingMode rounding_mode)
2418 {
2419   decimalInfo D;
2420   opStatus fs;
2421
2422   /* Scan the text.  */
2423   StringRef::iterator p = str.begin();
2424   interpretDecimal(p, str.end(), &D);
2425
2426   /* Handle the quick cases.  First the case of no significant digits,
2427      i.e. zero, and then exponents that are obviously too large or too
2428      small.  Writing L for log 10 / log 2, a number d.ddddd*10^exp
2429      definitely overflows if
2430
2431            (exp - 1) * L >= maxExponent
2432
2433      and definitely underflows to zero where
2434
2435            (exp + 1) * L <= minExponent - precision
2436
2437      With integer arithmetic the tightest bounds for L are
2438
2439            93/28 < L < 196/59            [ numerator <= 256 ]
2440            42039/12655 < L < 28738/8651  [ numerator <= 65536 ]
2441   */
2442
2443   if (decDigitValue(*D.firstSigDigit) >= 10U) {
2444     category = fcZero;
2445     fs = opOK;
2446
2447   /* Check whether the normalized exponent is high enough to overflow
2448      max during the log-rebasing in the max-exponent check below. */
2449   } else if (D.normalizedExponent - 1 > INT_MAX / 42039) {
2450     fs = handleOverflow(rounding_mode);
2451
2452   /* If it wasn't, then it also wasn't high enough to overflow max
2453      during the log-rebasing in the min-exponent check.  Check that it
2454      won't overflow min in either check, then perform the min-exponent
2455      check. */
2456   } else if (D.normalizedExponent - 1 < INT_MIN / 42039 ||
2457              (D.normalizedExponent + 1) * 28738 <=
2458                8651 * (semantics->minExponent - (int) semantics->precision)) {
2459     /* Underflow to zero and round.  */
2460     zeroSignificand();
2461     fs = normalize(rounding_mode, lfLessThanHalf);
2462
2463   /* We can finally safely perform the max-exponent check. */
2464   } else if ((D.normalizedExponent - 1) * 42039
2465              >= 12655 * semantics->maxExponent) {
2466     /* Overflow and round.  */
2467     fs = handleOverflow(rounding_mode);
2468   } else {
2469     integerPart *decSignificand;
2470     unsigned int partCount;
2471
2472     /* A tight upper bound on number of bits required to hold an
2473        N-digit decimal integer is N * 196 / 59.  Allocate enough space
2474        to hold the full significand, and an extra part required by
2475        tcMultiplyPart.  */
2476     partCount = static_cast<unsigned int>(D.lastSigDigit - D.firstSigDigit) + 1;
2477     partCount = partCountForBits(1 + 196 * partCount / 59);
2478     decSignificand = new integerPart[partCount + 1];
2479     partCount = 0;
2480
2481     /* Convert to binary efficiently - we do almost all multiplication
2482        in an integerPart.  When this would overflow do we do a single
2483        bignum multiplication, and then revert again to multiplication
2484        in an integerPart.  */
2485     do {
2486       integerPart decValue, val, multiplier;
2487
2488       val = 0;
2489       multiplier = 1;
2490
2491       do {
2492         if (*p == '.') {
2493           p++;
2494           if (p == str.end()) {
2495             break;
2496           }
2497         }
2498         decValue = decDigitValue(*p++);
2499         assert(decValue < 10U && "Invalid character in significand");
2500         multiplier *= 10;
2501         val = val * 10 + decValue;
2502         /* The maximum number that can be multiplied by ten with any
2503            digit added without overflowing an integerPart.  */
2504       } while (p <= D.lastSigDigit && multiplier <= (~ (integerPart) 0 - 9) / 10);
2505
2506       /* Multiply out the current part.  */
2507       APInt::tcMultiplyPart(decSignificand, decSignificand, multiplier, val,
2508                             partCount, partCount + 1, false);
2509
2510       /* If we used another part (likely but not guaranteed), increase
2511          the count.  */
2512       if (decSignificand[partCount])
2513         partCount++;
2514     } while (p <= D.lastSigDigit);
2515
2516     category = fcNormal;
2517     fs = roundSignificandWithExponent(decSignificand, partCount,
2518                                       D.exponent, rounding_mode);
2519
2520     delete [] decSignificand;
2521   }
2522
2523   return fs;
2524 }
2525
2526 APFloat::opStatus
2527 APFloat::convertFromString(StringRef str, roundingMode rounding_mode)
2528 {
2529   assertArithmeticOK(*semantics);
2530   assert(!str.empty() && "Invalid string length");
2531
2532   /* Handle a leading minus sign.  */
2533   StringRef::iterator p = str.begin();
2534   size_t slen = str.size();
2535   sign = *p == '-' ? 1 : 0;
2536   if (*p == '-' || *p == '+') {
2537     p++;
2538     slen--;
2539     assert(slen && "String has no digits");
2540   }
2541
2542   if (slen >= 2 && p[0] == '0' && (p[1] == 'x' || p[1] == 'X')) {
2543     assert(slen - 2 && "Invalid string");
2544     return convertFromHexadecimalString(StringRef(p + 2, slen - 2),
2545                                         rounding_mode);
2546   }
2547
2548   return convertFromDecimalString(StringRef(p, slen), rounding_mode);
2549 }
2550
2551 /* Write out a hexadecimal representation of the floating point value
2552    to DST, which must be of sufficient size, in the C99 form
2553    [-]0xh.hhhhp[+-]d.  Return the number of characters written,
2554    excluding the terminating NUL.
2555
2556    If UPPERCASE, the output is in upper case, otherwise in lower case.
2557
2558    HEXDIGITS digits appear altogether, rounding the value if
2559    necessary.  If HEXDIGITS is 0, the minimal precision to display the
2560    number precisely is used instead.  If nothing would appear after
2561    the decimal point it is suppressed.
2562
2563    The decimal exponent is always printed and has at least one digit.
2564    Zero values display an exponent of zero.  Infinities and NaNs
2565    appear as "infinity" or "nan" respectively.
2566
2567    The above rules are as specified by C99.  There is ambiguity about
2568    what the leading hexadecimal digit should be.  This implementation
2569    uses whatever is necessary so that the exponent is displayed as
2570    stored.  This implies the exponent will fall within the IEEE format
2571    range, and the leading hexadecimal digit will be 0 (for denormals),
2572    1 (normal numbers) or 2 (normal numbers rounded-away-from-zero with
2573    any other digits zero).
2574 */
2575 unsigned int
2576 APFloat::convertToHexString(char *dst, unsigned int hexDigits,
2577                             bool upperCase, roundingMode rounding_mode) const
2578 {
2579   char *p;
2580
2581   assertArithmeticOK(*semantics);
2582
2583   p = dst;
2584   if (sign)
2585     *dst++ = '-';
2586
2587   switch (category) {
2588   case fcInfinity:
2589     memcpy (dst, upperCase ? infinityU: infinityL, sizeof infinityU - 1);
2590     dst += sizeof infinityL - 1;
2591     break;
2592
2593   case fcNaN:
2594     memcpy (dst, upperCase ? NaNU: NaNL, sizeof NaNU - 1);
2595     dst += sizeof NaNU - 1;
2596     break;
2597
2598   case fcZero:
2599     *dst++ = '0';
2600     *dst++ = upperCase ? 'X': 'x';
2601     *dst++ = '0';
2602     if (hexDigits > 1) {
2603       *dst++ = '.';
2604       memset (dst, '0', hexDigits - 1);
2605       dst += hexDigits - 1;
2606     }
2607     *dst++ = upperCase ? 'P': 'p';
2608     *dst++ = '0';
2609     break;
2610
2611   case fcNormal:
2612     dst = convertNormalToHexString (dst, hexDigits, upperCase, rounding_mode);
2613     break;
2614   }
2615
2616   *dst = 0;
2617
2618   return static_cast<unsigned int>(dst - p);
2619 }
2620
2621 /* Does the hard work of outputting the correctly rounded hexadecimal
2622    form of a normal floating point number with the specified number of
2623    hexadecimal digits.  If HEXDIGITS is zero the minimum number of
2624    digits necessary to print the value precisely is output.  */
2625 char *
2626 APFloat::convertNormalToHexString(char *dst, unsigned int hexDigits,
2627                                   bool upperCase,
2628                                   roundingMode rounding_mode) const
2629 {
2630   unsigned int count, valueBits, shift, partsCount, outputDigits;
2631   const char *hexDigitChars;
2632   const integerPart *significand;
2633   char *p;
2634   bool roundUp;
2635
2636   *dst++ = '0';
2637   *dst++ = upperCase ? 'X': 'x';
2638
2639   roundUp = false;
2640   hexDigitChars = upperCase ? hexDigitsUpper: hexDigitsLower;
2641
2642   significand = significandParts();
2643   partsCount = partCount();
2644
2645   /* +3 because the first digit only uses the single integer bit, so
2646      we have 3 virtual zero most-significant-bits.  */
2647   valueBits = semantics->precision + 3;
2648   shift = integerPartWidth - valueBits % integerPartWidth;
2649
2650   /* The natural number of digits required ignoring trailing
2651      insignificant zeroes.  */
2652   outputDigits = (valueBits - significandLSB () + 3) / 4;
2653
2654   /* hexDigits of zero means use the required number for the
2655      precision.  Otherwise, see if we are truncating.  If we are,
2656      find out if we need to round away from zero.  */
2657   if (hexDigits) {
2658     if (hexDigits < outputDigits) {
2659       /* We are dropping non-zero bits, so need to check how to round.
2660          "bits" is the number of dropped bits.  */
2661       unsigned int bits;
2662       lostFraction fraction;
2663
2664       bits = valueBits - hexDigits * 4;
2665       fraction = lostFractionThroughTruncation (significand, partsCount, bits);
2666       roundUp = roundAwayFromZero(rounding_mode, fraction, bits);
2667     }
2668     outputDigits = hexDigits;
2669   }
2670
2671   /* Write the digits consecutively, and start writing in the location
2672      of the hexadecimal point.  We move the most significant digit
2673      left and add the hexadecimal point later.  */
2674   p = ++dst;
2675
2676   count = (valueBits + integerPartWidth - 1) / integerPartWidth;
2677
2678   while (outputDigits && count) {
2679     integerPart part;
2680
2681     /* Put the most significant integerPartWidth bits in "part".  */
2682     if (--count == partsCount)
2683       part = 0;  /* An imaginary higher zero part.  */
2684     else
2685       part = significand[count] << shift;
2686
2687     if (count && shift)
2688       part |= significand[count - 1] >> (integerPartWidth - shift);
2689
2690     /* Convert as much of "part" to hexdigits as we can.  */
2691     unsigned int curDigits = integerPartWidth / 4;
2692
2693     if (curDigits > outputDigits)
2694       curDigits = outputDigits;
2695     dst += partAsHex (dst, part, curDigits, hexDigitChars);
2696     outputDigits -= curDigits;
2697   }
2698
2699   if (roundUp) {
2700     char *q = dst;
2701
2702     /* Note that hexDigitChars has a trailing '0'.  */
2703     do {
2704       q--;
2705       *q = hexDigitChars[hexDigitValue (*q) + 1];
2706     } while (*q == '0');
2707     assert(q >= p);
2708   } else {
2709     /* Add trailing zeroes.  */
2710     memset (dst, '0', outputDigits);
2711     dst += outputDigits;
2712   }
2713
2714   /* Move the most significant digit to before the point, and if there
2715      is something after the decimal point add it.  This must come
2716      after rounding above.  */
2717   p[-1] = p[0];
2718   if (dst -1 == p)
2719     dst--;
2720   else
2721     p[0] = '.';
2722
2723   /* Finally output the exponent.  */
2724   *dst++ = upperCase ? 'P': 'p';
2725
2726   return writeSignedDecimal (dst, exponent);
2727 }
2728
2729 hash_code llvm::hash_value(const APFloat &Arg) {
2730   if (Arg.category != APFloat::fcNormal)
2731     return hash_combine((uint8_t)Arg.category,
2732                         // NaN has no sign, fix it at zero.
2733                         Arg.isNaN() ? (uint8_t)0 : (uint8_t)Arg.sign,
2734                         Arg.semantics->precision);
2735
2736   // Normal floats need their exponent and significand hashed.
2737   return hash_combine((uint8_t)Arg.category, (uint8_t)Arg.sign,
2738                       Arg.semantics->precision, Arg.exponent,
2739                       hash_combine_range(
2740                         Arg.significandParts(),
2741                         Arg.significandParts() + Arg.partCount()));
2742 }
2743
2744 // Conversion from APFloat to/from host float/double.  It may eventually be
2745 // possible to eliminate these and have everybody deal with APFloats, but that
2746 // will take a while.  This approach will not easily extend to long double.
2747 // Current implementation requires integerPartWidth==64, which is correct at
2748 // the moment but could be made more general.
2749
2750 // Denormals have exponent minExponent in APFloat, but minExponent-1 in
2751 // the actual IEEE respresentations.  We compensate for that here.
2752
2753 APInt
2754 APFloat::convertF80LongDoubleAPFloatToAPInt() const
2755 {
2756   assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended);
2757   assert(partCount()==2);
2758
2759   uint64_t myexponent, mysignificand;
2760
2761   if (category==fcNormal) {
2762     myexponent = exponent+16383; //bias
2763     mysignificand = significandParts()[0];
2764     if (myexponent==1 && !(mysignificand & 0x8000000000000000ULL))
2765       myexponent = 0;   // denormal
2766   } else if (category==fcZero) {
2767     myexponent = 0;
2768     mysignificand = 0;
2769   } else if (category==fcInfinity) {
2770     myexponent = 0x7fff;
2771     mysignificand = 0x8000000000000000ULL;
2772   } else {
2773     assert(category == fcNaN && "Unknown category");
2774     myexponent = 0x7fff;
2775     mysignificand = significandParts()[0];
2776   }
2777
2778   uint64_t words[2];
2779   words[0] = mysignificand;
2780   words[1] =  ((uint64_t)(sign & 1) << 15) |
2781               (myexponent & 0x7fffLL);
2782   return APInt(80, words);
2783 }
2784
2785 APInt
2786 APFloat::convertPPCDoubleDoubleAPFloatToAPInt() const
2787 {
2788   assert(semantics == (const llvm::fltSemantics*)&PPCDoubleDouble);
2789   assert(partCount()==2);
2790
2791   uint64_t myexponent, mysignificand, myexponent2, mysignificand2;
2792
2793   if (category==fcNormal) {
2794     myexponent = exponent + 1023; //bias
2795     myexponent2 = exponent2 + 1023;
2796     mysignificand = significandParts()[0];
2797     mysignificand2 = significandParts()[1];
2798     if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
2799       myexponent = 0;   // denormal
2800     if (myexponent2==1 && !(mysignificand2 & 0x10000000000000LL))
2801       myexponent2 = 0;   // denormal
2802   } else if (category==fcZero) {
2803     myexponent = 0;
2804     mysignificand = 0;
2805     myexponent2 = 0;
2806     mysignificand2 = 0;
2807   } else if (category==fcInfinity) {
2808     myexponent = 0x7ff;
2809     myexponent2 = 0;
2810     mysignificand = 0;
2811     mysignificand2 = 0;
2812   } else {
2813     assert(category == fcNaN && "Unknown category");
2814     myexponent = 0x7ff;
2815     mysignificand = significandParts()[0];
2816     myexponent2 = exponent2;
2817     mysignificand2 = significandParts()[1];
2818   }
2819
2820   uint64_t words[2];
2821   words[0] =  ((uint64_t)(sign & 1) << 63) |
2822               ((myexponent & 0x7ff) <<  52) |
2823               (mysignificand & 0xfffffffffffffLL);
2824   words[1] =  ((uint64_t)(sign2 & 1) << 63) |
2825               ((myexponent2 & 0x7ff) <<  52) |
2826               (mysignificand2 & 0xfffffffffffffLL);
2827   return APInt(128, words);
2828 }
2829
2830 APInt
2831 APFloat::convertQuadrupleAPFloatToAPInt() const
2832 {
2833   assert(semantics == (const llvm::fltSemantics*)&IEEEquad);
2834   assert(partCount()==2);
2835
2836   uint64_t myexponent, mysignificand, mysignificand2;
2837
2838   if (category==fcNormal) {
2839     myexponent = exponent+16383; //bias
2840     mysignificand = significandParts()[0];
2841     mysignificand2 = significandParts()[1];
2842     if (myexponent==1 && !(mysignificand2 & 0x1000000000000LL))
2843       myexponent = 0;   // denormal
2844   } else if (category==fcZero) {
2845     myexponent = 0;
2846     mysignificand = mysignificand2 = 0;
2847   } else if (category==fcInfinity) {
2848     myexponent = 0x7fff;
2849     mysignificand = mysignificand2 = 0;
2850   } else {
2851     assert(category == fcNaN && "Unknown category!");
2852     myexponent = 0x7fff;
2853     mysignificand = significandParts()[0];
2854     mysignificand2 = significandParts()[1];
2855   }
2856
2857   uint64_t words[2];
2858   words[0] = mysignificand;
2859   words[1] = ((uint64_t)(sign & 1) << 63) |
2860              ((myexponent & 0x7fff) << 48) |
2861              (mysignificand2 & 0xffffffffffffLL);
2862
2863   return APInt(128, words);
2864 }
2865
2866 APInt
2867 APFloat::convertDoubleAPFloatToAPInt() const
2868 {
2869   assert(semantics == (const llvm::fltSemantics*)&IEEEdouble);
2870   assert(partCount()==1);
2871
2872   uint64_t myexponent, mysignificand;
2873
2874   if (category==fcNormal) {
2875     myexponent = exponent+1023; //bias
2876     mysignificand = *significandParts();
2877     if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
2878       myexponent = 0;   // denormal
2879   } else if (category==fcZero) {
2880     myexponent = 0;
2881     mysignificand = 0;
2882   } else if (category==fcInfinity) {
2883     myexponent = 0x7ff;
2884     mysignificand = 0;
2885   } else {
2886     assert(category == fcNaN && "Unknown category!");
2887     myexponent = 0x7ff;
2888     mysignificand = *significandParts();
2889   }
2890
2891   return APInt(64, ((((uint64_t)(sign & 1) << 63) |
2892                      ((myexponent & 0x7ff) <<  52) |
2893                      (mysignificand & 0xfffffffffffffLL))));
2894 }
2895
2896 APInt
2897 APFloat::convertFloatAPFloatToAPInt() const
2898 {
2899   assert(semantics == (const llvm::fltSemantics*)&IEEEsingle);
2900   assert(partCount()==1);
2901
2902   uint32_t myexponent, mysignificand;
2903
2904   if (category==fcNormal) {
2905     myexponent = exponent+127; //bias
2906     mysignificand = (uint32_t)*significandParts();
2907     if (myexponent == 1 && !(mysignificand & 0x800000))
2908       myexponent = 0;   // denormal
2909   } else if (category==fcZero) {
2910     myexponent = 0;
2911     mysignificand = 0;
2912   } else if (category==fcInfinity) {
2913     myexponent = 0xff;
2914     mysignificand = 0;
2915   } else {
2916     assert(category == fcNaN && "Unknown category!");
2917     myexponent = 0xff;
2918     mysignificand = (uint32_t)*significandParts();
2919   }
2920
2921   return APInt(32, (((sign&1) << 31) | ((myexponent&0xff) << 23) |
2922                     (mysignificand & 0x7fffff)));
2923 }
2924
2925 APInt
2926 APFloat::convertHalfAPFloatToAPInt() const
2927 {
2928   assert(semantics == (const llvm::fltSemantics*)&IEEEhalf);
2929   assert(partCount()==1);
2930
2931   uint32_t myexponent, mysignificand;
2932
2933   if (category==fcNormal) {
2934     myexponent = exponent+15; //bias
2935     mysignificand = (uint32_t)*significandParts();
2936     if (myexponent == 1 && !(mysignificand & 0x400))
2937       myexponent = 0;   // denormal
2938   } else if (category==fcZero) {
2939     myexponent = 0;
2940     mysignificand = 0;
2941   } else if (category==fcInfinity) {
2942     myexponent = 0x1f;
2943     mysignificand = 0;
2944   } else {
2945     assert(category == fcNaN && "Unknown category!");
2946     myexponent = 0x1f;
2947     mysignificand = (uint32_t)*significandParts();
2948   }
2949
2950   return APInt(16, (((sign&1) << 15) | ((myexponent&0x1f) << 10) |
2951                     (mysignificand & 0x3ff)));
2952 }
2953
2954 // This function creates an APInt that is just a bit map of the floating
2955 // point constant as it would appear in memory.  It is not a conversion,
2956 // and treating the result as a normal integer is unlikely to be useful.
2957
2958 APInt
2959 APFloat::bitcastToAPInt() const
2960 {
2961   if (semantics == (const llvm::fltSemantics*)&IEEEhalf)
2962     return convertHalfAPFloatToAPInt();
2963
2964   if (semantics == (const llvm::fltSemantics*)&IEEEsingle)
2965     return convertFloatAPFloatToAPInt();
2966
2967   if (semantics == (const llvm::fltSemantics*)&IEEEdouble)
2968     return convertDoubleAPFloatToAPInt();
2969
2970   if (semantics == (const llvm::fltSemantics*)&IEEEquad)
2971     return convertQuadrupleAPFloatToAPInt();
2972
2973   if (semantics == (const llvm::fltSemantics*)&PPCDoubleDouble)
2974     return convertPPCDoubleDoubleAPFloatToAPInt();
2975
2976   assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended &&
2977          "unknown format!");
2978   return convertF80LongDoubleAPFloatToAPInt();
2979 }
2980
2981 float
2982 APFloat::convertToFloat() const
2983 {
2984   assert(semantics == (const llvm::fltSemantics*)&IEEEsingle &&
2985          "Float semantics are not IEEEsingle");
2986   APInt api = bitcastToAPInt();
2987   return api.bitsToFloat();
2988 }
2989
2990 double
2991 APFloat::convertToDouble() const
2992 {
2993   assert(semantics == (const llvm::fltSemantics*)&IEEEdouble &&
2994          "Float semantics are not IEEEdouble");
2995   APInt api = bitcastToAPInt();
2996   return api.bitsToDouble();
2997 }
2998
2999 /// Integer bit is explicit in this format.  Intel hardware (387 and later)
3000 /// does not support these bit patterns:
3001 ///  exponent = all 1's, integer bit 0, significand 0 ("pseudoinfinity")
3002 ///  exponent = all 1's, integer bit 0, significand nonzero ("pseudoNaN")
3003 ///  exponent = 0, integer bit 1 ("pseudodenormal")
3004 ///  exponent!=0 nor all 1's, integer bit 0 ("unnormal")
3005 /// At the moment, the first two are treated as NaNs, the second two as Normal.
3006 void
3007 APFloat::initFromF80LongDoubleAPInt(const APInt &api)
3008 {
3009   assert(api.getBitWidth()==80);
3010   uint64_t i1 = api.getRawData()[0];
3011   uint64_t i2 = api.getRawData()[1];
3012   uint64_t myexponent = (i2 & 0x7fff);
3013   uint64_t mysignificand = i1;
3014
3015   initialize(&APFloat::x87DoubleExtended);
3016   assert(partCount()==2);
3017
3018   sign = static_cast<unsigned int>(i2>>15);
3019   if (myexponent==0 && mysignificand==0) {
3020     // exponent, significand meaningless
3021     category = fcZero;
3022   } else if (myexponent==0x7fff && mysignificand==0x8000000000000000ULL) {
3023     // exponent, significand meaningless
3024     category = fcInfinity;
3025   } else if (myexponent==0x7fff && mysignificand!=0x8000000000000000ULL) {
3026     // exponent meaningless
3027     category = fcNaN;
3028     significandParts()[0] = mysignificand;
3029     significandParts()[1] = 0;
3030   } else {
3031     category = fcNormal;
3032     exponent = myexponent - 16383;
3033     significandParts()[0] = mysignificand;
3034     significandParts()[1] = 0;
3035     if (myexponent==0)          // denormal
3036       exponent = -16382;
3037   }
3038 }
3039
3040 void
3041 APFloat::initFromPPCDoubleDoubleAPInt(const APInt &api)
3042 {
3043   assert(api.getBitWidth()==128);
3044   uint64_t i1 = api.getRawData()[0];
3045   uint64_t i2 = api.getRawData()[1];
3046   uint64_t myexponent = (i1 >> 52) & 0x7ff;
3047   uint64_t mysignificand = i1 & 0xfffffffffffffLL;
3048   uint64_t myexponent2 = (i2 >> 52) & 0x7ff;
3049   uint64_t mysignificand2 = i2 & 0xfffffffffffffLL;
3050
3051   initialize(&APFloat::PPCDoubleDouble);
3052   assert(partCount()==2);
3053
3054   sign = static_cast<unsigned int>(i1>>63);
3055   sign2 = static_cast<unsigned int>(i2>>63);
3056   if (myexponent==0 && mysignificand==0) {
3057     // exponent, significand meaningless
3058     // exponent2 and significand2 are required to be 0; we don't check
3059     category = fcZero;
3060   } else if (myexponent==0x7ff && mysignificand==0) {
3061     // exponent, significand meaningless
3062     // exponent2 and significand2 are required to be 0; we don't check
3063     category = fcInfinity;
3064   } else if (myexponent==0x7ff && mysignificand!=0) {
3065     // exponent meaningless.  So is the whole second word, but keep it
3066     // for determinism.
3067     category = fcNaN;
3068     exponent2 = myexponent2;
3069     significandParts()[0] = mysignificand;
3070     significandParts()[1] = mysignificand2;
3071   } else {
3072     category = fcNormal;
3073     // Note there is no category2; the second word is treated as if it is
3074     // fcNormal, although it might be something else considered by itself.
3075     exponent = myexponent - 1023;
3076     exponent2 = myexponent2 - 1023;
3077     significandParts()[0] = mysignificand;
3078     significandParts()[1] = mysignificand2;
3079     if (myexponent==0)          // denormal
3080       exponent = -1022;
3081     else
3082       significandParts()[0] |= 0x10000000000000LL;  // integer bit
3083     if (myexponent2==0)
3084       exponent2 = -1022;
3085     else
3086       significandParts()[1] |= 0x10000000000000LL;  // integer bit
3087   }
3088 }
3089
3090 void
3091 APFloat::initFromQuadrupleAPInt(const APInt &api)
3092 {
3093   assert(api.getBitWidth()==128);
3094   uint64_t i1 = api.getRawData()[0];
3095   uint64_t i2 = api.getRawData()[1];
3096   uint64_t myexponent = (i2 >> 48) & 0x7fff;
3097   uint64_t mysignificand  = i1;
3098   uint64_t mysignificand2 = i2 & 0xffffffffffffLL;
3099
3100   initialize(&APFloat::IEEEquad);
3101   assert(partCount()==2);
3102
3103   sign = static_cast<unsigned int>(i2>>63);
3104   if (myexponent==0 &&
3105       (mysignificand==0 && mysignificand2==0)) {
3106     // exponent, significand meaningless
3107     category = fcZero;
3108   } else if (myexponent==0x7fff &&
3109              (mysignificand==0 && mysignificand2==0)) {
3110     // exponent, significand meaningless
3111     category = fcInfinity;
3112   } else if (myexponent==0x7fff &&
3113              (mysignificand!=0 || mysignificand2 !=0)) {
3114     // exponent meaningless
3115     category = fcNaN;
3116     significandParts()[0] = mysignificand;
3117     significandParts()[1] = mysignificand2;
3118   } else {
3119     category = fcNormal;
3120     exponent = myexponent - 16383;
3121     significandParts()[0] = mysignificand;
3122     significandParts()[1] = mysignificand2;
3123     if (myexponent==0)          // denormal
3124       exponent = -16382;
3125     else
3126       significandParts()[1] |= 0x1000000000000LL;  // integer bit
3127   }
3128 }
3129
3130 void
3131 APFloat::initFromDoubleAPInt(const APInt &api)
3132 {
3133   assert(api.getBitWidth()==64);
3134   uint64_t i = *api.getRawData();
3135   uint64_t myexponent = (i >> 52) & 0x7ff;
3136   uint64_t mysignificand = i & 0xfffffffffffffLL;
3137
3138   initialize(&APFloat::IEEEdouble);
3139   assert(partCount()==1);
3140
3141   sign = static_cast<unsigned int>(i>>63);
3142   if (myexponent==0 && mysignificand==0) {
3143     // exponent, significand meaningless
3144     category = fcZero;
3145   } else if (myexponent==0x7ff && mysignificand==0) {
3146     // exponent, significand meaningless
3147     category = fcInfinity;
3148   } else if (myexponent==0x7ff && mysignificand!=0) {
3149     // exponent meaningless
3150     category = fcNaN;
3151     *significandParts() = mysignificand;
3152   } else {
3153     category = fcNormal;
3154     exponent = myexponent - 1023;
3155     *significandParts() = mysignificand;
3156     if (myexponent==0)          // denormal
3157       exponent = -1022;
3158     else
3159       *significandParts() |= 0x10000000000000LL;  // integer bit
3160   }
3161 }
3162
3163 void
3164 APFloat::initFromFloatAPInt(const APInt & api)
3165 {
3166   assert(api.getBitWidth()==32);
3167   uint32_t i = (uint32_t)*api.getRawData();
3168   uint32_t myexponent = (i >> 23) & 0xff;
3169   uint32_t mysignificand = i & 0x7fffff;
3170
3171   initialize(&APFloat::IEEEsingle);
3172   assert(partCount()==1);
3173
3174   sign = i >> 31;
3175   if (myexponent==0 && mysignificand==0) {
3176     // exponent, significand meaningless
3177     category = fcZero;
3178   } else if (myexponent==0xff && mysignificand==0) {
3179     // exponent, significand meaningless
3180     category = fcInfinity;
3181   } else if (myexponent==0xff && mysignificand!=0) {
3182     // sign, exponent, significand meaningless
3183     category = fcNaN;
3184     *significandParts() = mysignificand;
3185   } else {
3186     category = fcNormal;
3187     exponent = myexponent - 127;  //bias
3188     *significandParts() = mysignificand;
3189     if (myexponent==0)    // denormal
3190       exponent = -126;
3191     else
3192       *significandParts() |= 0x800000; // integer bit
3193   }
3194 }
3195
3196 void
3197 APFloat::initFromHalfAPInt(const APInt & api)
3198 {
3199   assert(api.getBitWidth()==16);
3200   uint32_t i = (uint32_t)*api.getRawData();
3201   uint32_t myexponent = (i >> 10) & 0x1f;
3202   uint32_t mysignificand = i & 0x3ff;
3203
3204   initialize(&APFloat::IEEEhalf);
3205   assert(partCount()==1);
3206
3207   sign = i >> 15;
3208   if (myexponent==0 && mysignificand==0) {
3209     // exponent, significand meaningless
3210     category = fcZero;
3211   } else if (myexponent==0x1f && mysignificand==0) {
3212     // exponent, significand meaningless
3213     category = fcInfinity;
3214   } else if (myexponent==0x1f && mysignificand!=0) {
3215     // sign, exponent, significand meaningless
3216     category = fcNaN;
3217     *significandParts() = mysignificand;
3218   } else {
3219     category = fcNormal;
3220     exponent = myexponent - 15;  //bias
3221     *significandParts() = mysignificand;
3222     if (myexponent==0)    // denormal
3223       exponent = -14;
3224     else
3225       *significandParts() |= 0x400; // integer bit
3226   }
3227 }
3228
3229 /// Treat api as containing the bits of a floating point number.  Currently
3230 /// we infer the floating point type from the size of the APInt.  The
3231 /// isIEEE argument distinguishes between PPC128 and IEEE128 (not meaningful
3232 /// when the size is anything else).
3233 void
3234 APFloat::initFromAPInt(const APInt& api, bool isIEEE)
3235 {
3236   if (api.getBitWidth() == 16)
3237     return initFromHalfAPInt(api);
3238   else if (api.getBitWidth() == 32)
3239     return initFromFloatAPInt(api);
3240   else if (api.getBitWidth()==64)
3241     return initFromDoubleAPInt(api);
3242   else if (api.getBitWidth()==80)
3243     return initFromF80LongDoubleAPInt(api);
3244   else if (api.getBitWidth()==128)
3245     return (isIEEE ?
3246             initFromQuadrupleAPInt(api) : initFromPPCDoubleDoubleAPInt(api));
3247   else
3248     llvm_unreachable(0);
3249 }
3250
3251 APFloat
3252 APFloat::getAllOnesValue(unsigned BitWidth, bool isIEEE)
3253 {
3254   return APFloat(APInt::getAllOnesValue(BitWidth), isIEEE);
3255 }
3256
3257 APFloat APFloat::getLargest(const fltSemantics &Sem, bool Negative) {
3258   APFloat Val(Sem, fcNormal, Negative);
3259
3260   // We want (in interchange format):
3261   //   sign = {Negative}
3262   //   exponent = 1..10
3263   //   significand = 1..1
3264
3265   Val.exponent = Sem.maxExponent; // unbiased
3266
3267   // 1-initialize all bits....
3268   Val.zeroSignificand();
3269   integerPart *significand = Val.significandParts();
3270   unsigned N = partCountForBits(Sem.precision);
3271   for (unsigned i = 0; i != N; ++i)
3272     significand[i] = ~((integerPart) 0);
3273
3274   // ...and then clear the top bits for internal consistency.
3275   if (Sem.precision % integerPartWidth != 0)
3276     significand[N-1] &=
3277       (((integerPart) 1) << (Sem.precision % integerPartWidth)) - 1;
3278
3279   return Val;
3280 }
3281
3282 APFloat APFloat::getSmallest(const fltSemantics &Sem, bool Negative) {
3283   APFloat Val(Sem, fcNormal, Negative);
3284
3285   // We want (in interchange format):
3286   //   sign = {Negative}
3287   //   exponent = 0..0
3288   //   significand = 0..01
3289
3290   Val.exponent = Sem.minExponent; // unbiased
3291   Val.zeroSignificand();
3292   Val.significandParts()[0] = 1;
3293   return Val;
3294 }
3295
3296 APFloat APFloat::getSmallestNormalized(const fltSemantics &Sem, bool Negative) {
3297   APFloat Val(Sem, fcNormal, Negative);
3298
3299   // We want (in interchange format):
3300   //   sign = {Negative}
3301   //   exponent = 0..0
3302   //   significand = 10..0
3303
3304   Val.exponent = Sem.minExponent;
3305   Val.zeroSignificand();
3306   Val.significandParts()[partCountForBits(Sem.precision)-1] |=
3307     (((integerPart) 1) << ((Sem.precision - 1) % integerPartWidth));
3308
3309   return Val;
3310 }
3311
3312 APFloat::APFloat(const APInt& api, bool isIEEE) : exponent2(0), sign2(0) {
3313   initFromAPInt(api, isIEEE);
3314 }
3315
3316 APFloat::APFloat(float f) : exponent2(0), sign2(0) {
3317   initFromAPInt(APInt::floatToBits(f));
3318 }
3319
3320 APFloat::APFloat(double d) : exponent2(0), sign2(0) {
3321   initFromAPInt(APInt::doubleToBits(d));
3322 }
3323
3324 namespace {
3325   void append(SmallVectorImpl<char> &Buffer, StringRef Str) {
3326     Buffer.append(Str.begin(), Str.end());
3327   }
3328
3329   /// Removes data from the given significand until it is no more
3330   /// precise than is required for the desired precision.
3331   void AdjustToPrecision(APInt &significand,
3332                          int &exp, unsigned FormatPrecision) {
3333     unsigned bits = significand.getActiveBits();
3334
3335     // 196/59 is a very slight overestimate of lg_2(10).
3336     unsigned bitsRequired = (FormatPrecision * 196 + 58) / 59;
3337
3338     if (bits <= bitsRequired) return;
3339
3340     unsigned tensRemovable = (bits - bitsRequired) * 59 / 196;
3341     if (!tensRemovable) return;
3342
3343     exp += tensRemovable;
3344
3345     APInt divisor(significand.getBitWidth(), 1);
3346     APInt powten(significand.getBitWidth(), 10);
3347     while (true) {
3348       if (tensRemovable & 1)
3349         divisor *= powten;
3350       tensRemovable >>= 1;
3351       if (!tensRemovable) break;
3352       powten *= powten;
3353     }
3354
3355     significand = significand.udiv(divisor);
3356
3357     // Truncate the significand down to its active bit count, but
3358     // don't try to drop below 32.
3359     unsigned newPrecision = std::max(32U, significand.getActiveBits());
3360     significand = significand.trunc(newPrecision);
3361   }
3362
3363
3364   void AdjustToPrecision(SmallVectorImpl<char> &buffer,
3365                          int &exp, unsigned FormatPrecision) {
3366     unsigned N = buffer.size();
3367     if (N <= FormatPrecision) return;
3368
3369     // The most significant figures are the last ones in the buffer.
3370     unsigned FirstSignificant = N - FormatPrecision;
3371
3372     // Round.
3373     // FIXME: this probably shouldn't use 'round half up'.
3374
3375     // Rounding down is just a truncation, except we also want to drop
3376     // trailing zeros from the new result.
3377     if (buffer[FirstSignificant - 1] < '5') {
3378       while (FirstSignificant < N && buffer[FirstSignificant] == '0')
3379         FirstSignificant++;
3380
3381       exp += FirstSignificant;
3382       buffer.erase(&buffer[0], &buffer[FirstSignificant]);
3383       return;
3384     }
3385
3386     // Rounding up requires a decimal add-with-carry.  If we continue
3387     // the carry, the newly-introduced zeros will just be truncated.
3388     for (unsigned I = FirstSignificant; I != N; ++I) {
3389       if (buffer[I] == '9') {
3390         FirstSignificant++;
3391       } else {
3392         buffer[I]++;
3393         break;
3394       }
3395     }
3396
3397     // If we carried through, we have exactly one digit of precision.
3398     if (FirstSignificant == N) {
3399       exp += FirstSignificant;
3400       buffer.clear();
3401       buffer.push_back('1');
3402       return;
3403     }
3404
3405     exp += FirstSignificant;
3406     buffer.erase(&buffer[0], &buffer[FirstSignificant]);
3407   }
3408 }
3409
3410 void APFloat::toString(SmallVectorImpl<char> &Str,
3411                        unsigned FormatPrecision,
3412                        unsigned FormatMaxPadding) const {
3413   switch (category) {
3414   case fcInfinity:
3415     if (isNegative())
3416       return append(Str, "-Inf");
3417     else
3418       return append(Str, "+Inf");
3419
3420   case fcNaN: return append(Str, "NaN");
3421
3422   case fcZero:
3423     if (isNegative())
3424       Str.push_back('-');
3425
3426     if (!FormatMaxPadding)
3427       append(Str, "0.0E+0");
3428     else
3429       Str.push_back('0');
3430     return;
3431
3432   case fcNormal:
3433     break;
3434   }
3435
3436   if (isNegative())
3437     Str.push_back('-');
3438
3439   // Decompose the number into an APInt and an exponent.
3440   int exp = exponent - ((int) semantics->precision - 1);
3441   APInt significand(semantics->precision,
3442                     makeArrayRef(significandParts(),
3443                                  partCountForBits(semantics->precision)));
3444
3445   // Set FormatPrecision if zero.  We want to do this before we
3446   // truncate trailing zeros, as those are part of the precision.
3447   if (!FormatPrecision) {
3448     // It's an interesting question whether to use the nominal
3449     // precision or the active precision here for denormals.
3450
3451     // FormatPrecision = ceil(significandBits / lg_2(10))
3452     FormatPrecision = (semantics->precision * 59 + 195) / 196;
3453   }
3454
3455   // Ignore trailing binary zeros.
3456   int trailingZeros = significand.countTrailingZeros();
3457   exp += trailingZeros;
3458   significand = significand.lshr(trailingZeros);
3459
3460   // Change the exponent from 2^e to 10^e.
3461   if (exp == 0) {
3462     // Nothing to do.
3463   } else if (exp > 0) {
3464     // Just shift left.
3465     significand = significand.zext(semantics->precision + exp);
3466     significand <<= exp;
3467     exp = 0;
3468   } else { /* exp < 0 */
3469     int texp = -exp;
3470
3471     // We transform this using the identity:
3472     //   (N)(2^-e) == (N)(5^e)(10^-e)
3473     // This means we have to multiply N (the significand) by 5^e.
3474     // To avoid overflow, we have to operate on numbers large
3475     // enough to store N * 5^e:
3476     //   log2(N * 5^e) == log2(N) + e * log2(5)
3477     //                 <= semantics->precision + e * 137 / 59
3478     //   (log_2(5) ~ 2.321928 < 2.322034 ~ 137/59)
3479
3480     unsigned precision = semantics->precision + (137 * texp + 136) / 59;
3481
3482     // Multiply significand by 5^e.
3483     //   N * 5^0101 == N * 5^(1*1) * 5^(0*2) * 5^(1*4) * 5^(0*8)
3484     significand = significand.zext(precision);
3485     APInt five_to_the_i(precision, 5);
3486     while (true) {
3487       if (texp & 1) significand *= five_to_the_i;
3488
3489       texp >>= 1;
3490       if (!texp) break;
3491       five_to_the_i *= five_to_the_i;
3492     }
3493   }
3494
3495   AdjustToPrecision(significand, exp, FormatPrecision);
3496
3497   llvm::SmallVector<char, 256> buffer;
3498
3499   // Fill the buffer.
3500   unsigned precision = significand.getBitWidth();
3501   APInt ten(precision, 10);
3502   APInt digit(precision, 0);
3503
3504   bool inTrail = true;
3505   while (significand != 0) {
3506     // digit <- significand % 10
3507     // significand <- significand / 10
3508     APInt::udivrem(significand, ten, significand, digit);
3509
3510     unsigned d = digit.getZExtValue();
3511
3512     // Drop trailing zeros.
3513     if (inTrail && !d) exp++;
3514     else {
3515       buffer.push_back((char) ('0' + d));
3516       inTrail = false;
3517     }
3518   }
3519
3520   assert(!buffer.empty() && "no characters in buffer!");
3521
3522   // Drop down to FormatPrecision.
3523   // TODO: don't do more precise calculations above than are required.
3524   AdjustToPrecision(buffer, exp, FormatPrecision);
3525
3526   unsigned NDigits = buffer.size();
3527
3528   // Check whether we should use scientific notation.
3529   bool FormatScientific;
3530   if (!FormatMaxPadding)
3531     FormatScientific = true;
3532   else {
3533     if (exp >= 0) {
3534       // 765e3 --> 765000
3535       //              ^^^
3536       // But we shouldn't make the number look more precise than it is.
3537       FormatScientific = ((unsigned) exp > FormatMaxPadding ||
3538                           NDigits + (unsigned) exp > FormatPrecision);
3539     } else {
3540       // Power of the most significant digit.
3541       int MSD = exp + (int) (NDigits - 1);
3542       if (MSD >= 0) {
3543         // 765e-2 == 7.65
3544         FormatScientific = false;
3545       } else {
3546         // 765e-5 == 0.00765
3547         //           ^ ^^
3548         FormatScientific = ((unsigned) -MSD) > FormatMaxPadding;
3549       }
3550     }
3551   }
3552
3553   // Scientific formatting is pretty straightforward.
3554   if (FormatScientific) {
3555     exp += (NDigits - 1);
3556
3557     Str.push_back(buffer[NDigits-1]);
3558     Str.push_back('.');
3559     if (NDigits == 1)
3560       Str.push_back('0');
3561     else
3562       for (unsigned I = 1; I != NDigits; ++I)
3563         Str.push_back(buffer[NDigits-1-I]);
3564     Str.push_back('E');
3565
3566     Str.push_back(exp >= 0 ? '+' : '-');
3567     if (exp < 0) exp = -exp;
3568     SmallVector<char, 6> expbuf;
3569     do {
3570       expbuf.push_back((char) ('0' + (exp % 10)));
3571       exp /= 10;
3572     } while (exp);
3573     for (unsigned I = 0, E = expbuf.size(); I != E; ++I)
3574       Str.push_back(expbuf[E-1-I]);
3575     return;
3576   }
3577
3578   // Non-scientific, positive exponents.
3579   if (exp >= 0) {
3580     for (unsigned I = 0; I != NDigits; ++I)
3581       Str.push_back(buffer[NDigits-1-I]);
3582     for (unsigned I = 0; I != (unsigned) exp; ++I)
3583       Str.push_back('0');
3584     return;
3585   }
3586
3587   // Non-scientific, negative exponents.
3588
3589   // The number of digits to the left of the decimal point.
3590   int NWholeDigits = exp + (int) NDigits;
3591
3592   unsigned I = 0;
3593   if (NWholeDigits > 0) {
3594     for (; I != (unsigned) NWholeDigits; ++I)
3595       Str.push_back(buffer[NDigits-I-1]);
3596     Str.push_back('.');
3597   } else {
3598     unsigned NZeros = 1 + (unsigned) -NWholeDigits;
3599
3600     Str.push_back('0');
3601     Str.push_back('.');
3602     for (unsigned Z = 1; Z != NZeros; ++Z)
3603       Str.push_back('0');
3604   }
3605
3606   for (; I != NDigits; ++I)
3607     Str.push_back(buffer[NDigits-I-1]);
3608 }
3609
3610 bool APFloat::getExactInverse(APFloat *inv) const {
3611   // We can only guarantee the existence of an exact inverse for IEEE floats.
3612   if (semantics != &IEEEhalf && semantics != &IEEEsingle &&
3613       semantics != &IEEEdouble && semantics != &IEEEquad)
3614     return false;
3615
3616   // Special floats and denormals have no exact inverse.
3617   if (category != fcNormal)
3618     return false;
3619
3620   // Check that the number is a power of two by making sure that only the
3621   // integer bit is set in the significand.
3622   if (significandLSB() != semantics->precision - 1)
3623     return false;
3624
3625   // Get the inverse.
3626   APFloat reciprocal(*semantics, 1ULL);
3627   if (reciprocal.divide(*this, rmNearestTiesToEven) != opOK)
3628     return false;
3629
3630   // Avoid multiplication with a denormal, it is not safe on all platforms and
3631   // may be slower than a normal division.
3632   if (reciprocal.significandMSB() + 1 < reciprocal.semantics->precision)
3633     return false;
3634
3635   assert(reciprocal.category == fcNormal &&
3636          reciprocal.significandLSB() == reciprocal.semantics->precision - 1);
3637
3638   if (inv)
3639     *inv = reciprocal;
3640
3641   return true;
3642 }