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