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