Oops, should be part of 41664; won't work very well without this piece.
[oota-llvm.git] / lib / Support / APFloat.cpp
1 //===-- APFloat.cpp - Implement APFloat class -----------------------------===//
2 //
3 //                     The LLVM Compiler Infrastructure
4 //
5 // This file was developed by Neil Booth and is distributed under the
6 // University of Illinois Open Source 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 <cassert>
16 #include "llvm/ADT/APFloat.h"
17 #include "llvm/Support/MathExtras.h"
18
19 using namespace llvm;
20
21 #define convolve(lhs, rhs) ((lhs) * 4 + (rhs))
22
23 /* Assumed in hexadecimal significand parsing.  */
24 COMPILE_TIME_ASSERT(integerPartWidth % 4 == 0);
25
26 namespace llvm {
27
28   /* Represents floating point arithmetic semantics.  */
29   struct fltSemantics {
30     /* The largest E such that 2^E is representable; this matches the
31        definition of IEEE 754.  */
32     exponent_t maxExponent;
33
34     /* The smallest E such that 2^E is a normalized number; this
35        matches the definition of IEEE 754.  */
36     exponent_t minExponent;
37
38     /* Number of bits in the significand.  This includes the integer
39        bit.  */
40     unsigned char precision;
41
42     /* If the target format has an implicit integer bit.  */
43     bool implicitIntegerBit;
44   };
45
46   const fltSemantics APFloat::IEEEsingle = { 127, -126, 24, true };
47   const fltSemantics APFloat::IEEEdouble = { 1023, -1022, 53, true };
48   const fltSemantics APFloat::IEEEquad = { 16383, -16382, 113, true };
49   const fltSemantics APFloat::x87DoubleExtended = { 16383, -16382, 64, false };
50   const fltSemantics APFloat::Bogus = { 0, 0, 0, false };
51 }
52
53 /* Put a bunch of private, handy routines in an anonymous namespace.  */
54 namespace {
55
56   inline unsigned int
57   partCountForBits(unsigned int bits)
58   {
59     return ((bits) + integerPartWidth - 1) / integerPartWidth;
60   }
61
62   unsigned int
63   digitValue(unsigned int c)
64   {
65     unsigned int r;
66
67     r = c - '0';
68     if(r <= 9)
69       return r;
70
71     return -1U;
72   }
73
74   unsigned int
75   hexDigitValue (unsigned int c)
76   {
77     unsigned int r;
78
79     r = c - '0';
80     if(r <= 9)
81       return r;
82
83     r = c - 'A';
84     if(r <= 5)
85       return r + 10;
86
87     r = c - 'a';
88     if(r <= 5)
89       return r + 10;
90
91     return -1U;
92   }
93
94   /* This is ugly and needs cleaning up, but I don't immediately see
95      how whilst remaining safe.  */
96   static int
97   totalExponent(const char *p, int exponentAdjustment)
98   {
99     integerPart unsignedExponent;
100     bool negative, overflow;
101     long exponent;
102
103     /* Move past the exponent letter and sign to the digits.  */
104     p++;
105     negative = *p == '-';
106     if(*p == '-' || *p == '+')
107       p++;
108
109     unsignedExponent = 0;
110     overflow = false;
111     for(;;) {
112       unsigned int value;
113
114       value = digitValue(*p);
115       if(value == -1U)
116         break;
117
118       p++;
119       unsignedExponent = unsignedExponent * 10 + value;
120       if(unsignedExponent > 65535)
121         overflow = true;
122     }
123
124     if(exponentAdjustment > 65535 || exponentAdjustment < -65536)
125       overflow = true;
126
127     if(!overflow) {
128       exponent = unsignedExponent;
129       if(negative)
130         exponent = -exponent;
131       exponent += exponentAdjustment;
132       if(exponent > 65535 || exponent < -65536)
133         overflow = true;
134     }
135
136     if(overflow)
137       exponent = negative ? -65536: 65535;
138
139     return exponent;
140   }
141
142   const char *
143   skipLeadingZeroesAndAnyDot(const char *p, const char **dot)
144   {
145     *dot = 0;
146     while(*p == '0')
147       p++;
148
149     if(*p == '.') {
150       *dot = p++;
151       while(*p == '0')
152         p++;
153     }
154
155     return p;
156   }
157
158   /* Return the trailing fraction of a hexadecimal number.
159      DIGITVALUE is the first hex digit of the fraction, P points to
160      the next digit.  */
161   lostFraction
162   trailingHexadecimalFraction(const char *p, unsigned int digitValue)
163   {
164     unsigned int hexDigit;
165
166     /* If the first trailing digit isn't 0 or 8 we can work out the
167        fraction immediately.  */
168     if(digitValue > 8)
169       return lfMoreThanHalf;
170     else if(digitValue < 8 && digitValue > 0)
171       return lfLessThanHalf;
172
173     /* Otherwise we need to find the first non-zero digit.  */
174     while(*p == '0')
175       p++;
176
177     hexDigit = hexDigitValue(*p);
178
179     /* If we ran off the end it is exactly zero or one-half, otherwise
180        a little more.  */
181     if(hexDigit == -1U)
182       return digitValue == 0 ? lfExactlyZero: lfExactlyHalf;
183     else
184       return digitValue == 0 ? lfLessThanHalf: lfMoreThanHalf;
185   }
186
187   /* Return the fraction lost were a bignum truncated.  */
188   lostFraction
189   lostFractionThroughTruncation(integerPart *parts,
190                                 unsigned int partCount,
191                                 unsigned int bits)
192   {
193     unsigned int lsb;
194
195     lsb = APInt::tcLSB(parts, partCount);
196
197     /* Note this is guaranteed true if bits == 0, or LSB == -1U.  */
198     if(bits <= lsb)
199       return lfExactlyZero;
200     if(bits == lsb + 1)
201       return lfExactlyHalf;
202     if(bits <= partCount * integerPartWidth
203        && APInt::tcExtractBit(parts, bits - 1))
204       return lfMoreThanHalf;
205
206     return lfLessThanHalf;
207   }
208
209   /* Shift DST right BITS bits noting lost fraction.  */
210   lostFraction
211   shiftRight(integerPart *dst, unsigned int parts, unsigned int bits)
212   {
213     lostFraction lost_fraction;
214
215     lost_fraction = lostFractionThroughTruncation(dst, parts, bits);
216
217     APInt::tcShiftRight(dst, parts, bits);
218
219     return lost_fraction;
220   }
221 }
222
223 /* Constructors.  */
224 void
225 APFloat::initialize(const fltSemantics *ourSemantics)
226 {
227   unsigned int count;
228
229   semantics = ourSemantics;
230   count = partCount();
231   if(count > 1)
232     significand.parts = new integerPart[count];
233 }
234
235 void
236 APFloat::freeSignificand()
237 {
238   if(partCount() > 1)
239     delete [] significand.parts;
240 }
241
242 void
243 APFloat::assign(const APFloat &rhs)
244 {
245   assert(semantics == rhs.semantics);
246
247   sign = rhs.sign;
248   category = rhs.category;
249   exponent = rhs.exponent;
250   if(category == fcNormal || category == fcNaN)
251     copySignificand(rhs);
252 }
253
254 void
255 APFloat::copySignificand(const APFloat &rhs)
256 {
257   assert(category == fcNormal || category == fcNaN);
258   assert(rhs.partCount() >= partCount());
259
260   APInt::tcAssign(significandParts(), rhs.significandParts(),
261                   partCount());
262 }
263
264 APFloat &
265 APFloat::operator=(const APFloat &rhs)
266 {
267   if(this != &rhs) {
268     if(semantics != rhs.semantics) {
269       freeSignificand();
270       initialize(rhs.semantics);
271     }
272     assign(rhs);
273   }
274
275   return *this;
276 }
277
278 bool
279 APFloat::bitwiseIsEqual(const APFloat &rhs) const {
280   if (this == &rhs)
281     return true;
282   if (semantics != rhs.semantics ||
283       category != rhs.category ||
284       sign != rhs.sign)
285     return false;
286   if (category==fcZero || category==fcInfinity)
287     return true;
288   else if (category==fcNormal && exponent!=rhs.exponent)
289     return false;
290   else {
291     int i= partCount();
292     const integerPart* p=significandParts();
293     const integerPart* q=rhs.significandParts();
294     for (; i>0; i--, p++, q++) {
295       if (*p != *q)
296         return false;
297     }
298     return true;
299   }
300 }
301
302 APFloat::APFloat(const fltSemantics &ourSemantics, integerPart value)
303 {
304   initialize(&ourSemantics);
305   sign = 0;
306   zeroSignificand();
307   exponent = ourSemantics.precision - 1;
308   significandParts()[0] = value;
309   normalize(rmNearestTiesToEven, lfExactlyZero);
310 }
311
312 APFloat::APFloat(const fltSemantics &ourSemantics,
313                  fltCategory ourCategory, bool negative)
314 {
315   initialize(&ourSemantics);
316   category = ourCategory;
317   sign = negative;
318   if(category == fcNormal)
319     category = fcZero;
320 }
321
322 APFloat::APFloat(const fltSemantics &ourSemantics, const char *text)
323 {
324   initialize(&ourSemantics);
325   convertFromString(text, rmNearestTiesToEven);
326 }
327
328 APFloat::APFloat(const APFloat &rhs)
329 {
330   initialize(rhs.semantics);
331   assign(rhs);
332 }
333
334 APFloat::~APFloat()
335 {
336   freeSignificand();
337 }
338
339 unsigned int
340 APFloat::partCount() const
341 {
342   return partCountForBits(semantics->precision + 1);
343 }
344
345 unsigned int
346 APFloat::semanticsPrecision(const fltSemantics &semantics)
347 {
348   return semantics.precision;
349 }
350
351 const integerPart *
352 APFloat::significandParts() const
353 {
354   return const_cast<APFloat *>(this)->significandParts();
355 }
356
357 integerPart *
358 APFloat::significandParts()
359 {
360   assert(category == fcNormal || category == fcNaN);
361
362   if(partCount() > 1)
363     return significand.parts;
364   else
365     return &significand.part;
366 }
367
368 /* Combine the effect of two lost fractions.  */
369 lostFraction
370 APFloat::combineLostFractions(lostFraction moreSignificant,
371                               lostFraction lessSignificant)
372 {
373   if(lessSignificant != lfExactlyZero) {
374     if(moreSignificant == lfExactlyZero)
375       moreSignificant = lfLessThanHalf;
376     else if(moreSignificant == lfExactlyHalf)
377       moreSignificant = lfMoreThanHalf;
378   }
379
380   return moreSignificant;
381 }
382
383 void
384 APFloat::zeroSignificand()
385 {
386   category = fcNormal;
387   APInt::tcSet(significandParts(), 0, partCount());
388 }
389
390 /* Increment an fcNormal floating point number's significand.  */
391 void
392 APFloat::incrementSignificand()
393 {
394   integerPart carry;
395
396   carry = APInt::tcIncrement(significandParts(), partCount());
397
398   /* Our callers should never cause us to overflow.  */
399   assert(carry == 0);
400 }
401
402 /* Add the significand of the RHS.  Returns the carry flag.  */
403 integerPart
404 APFloat::addSignificand(const APFloat &rhs)
405 {
406   integerPart *parts;
407
408   parts = significandParts();
409
410   assert(semantics == rhs.semantics);
411   assert(exponent == rhs.exponent);
412
413   return APInt::tcAdd(parts, rhs.significandParts(), 0, partCount());
414 }
415
416 /* Subtract the significand of the RHS with a borrow flag.  Returns
417    the borrow flag.  */
418 integerPart
419 APFloat::subtractSignificand(const APFloat &rhs, integerPart borrow)
420 {
421   integerPart *parts;
422
423   parts = significandParts();
424
425   assert(semantics == rhs.semantics);
426   assert(exponent == rhs.exponent);
427
428   return APInt::tcSubtract(parts, rhs.significandParts(), borrow,
429                            partCount());
430 }
431
432 /* Multiply the significand of the RHS.  If ADDEND is non-NULL, add it
433    on to the full-precision result of the multiplication.  Returns the
434    lost fraction.  */
435 lostFraction
436 APFloat::multiplySignificand(const APFloat &rhs, const APFloat *addend)
437 {
438   unsigned int omsb;    // One, not zero, based MSB.
439   unsigned int partsCount, newPartsCount, precision;
440   integerPart *lhsSignificand;
441   integerPart scratch[4];
442   integerPart *fullSignificand;
443   lostFraction lost_fraction;
444
445   assert(semantics == rhs.semantics);
446
447   precision = semantics->precision;
448   newPartsCount = partCountForBits(precision * 2);
449
450   if(newPartsCount > 4)
451     fullSignificand = new integerPart[newPartsCount];
452   else
453     fullSignificand = scratch;
454
455   lhsSignificand = significandParts();
456   partsCount = partCount();
457
458   APInt::tcFullMultiply(fullSignificand, lhsSignificand,
459                         rhs.significandParts(), partsCount);
460
461   lost_fraction = lfExactlyZero;
462   omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
463   exponent += rhs.exponent;
464
465   if(addend) {
466     Significand savedSignificand = significand;
467     const fltSemantics *savedSemantics = semantics;
468     fltSemantics extendedSemantics;
469     opStatus status;
470     unsigned int extendedPrecision;
471
472     /* Normalize our MSB.  */
473     extendedPrecision = precision + precision - 1;
474     if(omsb != extendedPrecision)
475       {
476         APInt::tcShiftLeft(fullSignificand, newPartsCount,
477                            extendedPrecision - omsb);
478         exponent -= extendedPrecision - omsb;
479       }
480
481     /* Create new semantics.  */
482     extendedSemantics = *semantics;
483     extendedSemantics.precision = extendedPrecision;
484
485     if(newPartsCount == 1)
486       significand.part = fullSignificand[0];
487     else
488       significand.parts = fullSignificand;
489     semantics = &extendedSemantics;
490
491     APFloat extendedAddend(*addend);
492     status = extendedAddend.convert(extendedSemantics, rmTowardZero);
493     assert(status == opOK);
494     lost_fraction = addOrSubtractSignificand(extendedAddend, false);
495
496     /* Restore our state.  */
497     if(newPartsCount == 1)
498       fullSignificand[0] = significand.part;
499     significand = savedSignificand;
500     semantics = savedSemantics;
501
502     omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
503   }
504
505   exponent -= (precision - 1);
506
507   if(omsb > precision) {
508     unsigned int bits, significantParts;
509     lostFraction lf;
510
511     bits = omsb - precision;
512     significantParts = partCountForBits(omsb);
513     lf = shiftRight(fullSignificand, significantParts, bits);
514     lost_fraction = combineLostFractions(lf, lost_fraction);
515     exponent += bits;
516   }
517
518   APInt::tcAssign(lhsSignificand, fullSignificand, partsCount);
519
520   if(newPartsCount > 4)
521     delete [] fullSignificand;
522
523   return lost_fraction;
524 }
525
526 /* Multiply the significands of LHS and RHS to DST.  */
527 lostFraction
528 APFloat::divideSignificand(const APFloat &rhs)
529 {
530   unsigned int bit, i, partsCount;
531   const integerPart *rhsSignificand;
532   integerPart *lhsSignificand, *dividend, *divisor;
533   integerPart scratch[4];
534   lostFraction lost_fraction;
535
536   assert(semantics == rhs.semantics);
537
538   lhsSignificand = significandParts();
539   rhsSignificand = rhs.significandParts();
540   partsCount = partCount();
541
542   if(partsCount > 2)
543     dividend = new integerPart[partsCount * 2];
544   else
545     dividend = scratch;
546
547   divisor = dividend + partsCount;
548
549   /* Copy the dividend and divisor as they will be modified in-place.  */
550   for(i = 0; i < partsCount; i++) {
551     dividend[i] = lhsSignificand[i];
552     divisor[i] = rhsSignificand[i];
553     lhsSignificand[i] = 0;
554   }
555
556   exponent -= rhs.exponent;
557
558   unsigned int precision = semantics->precision;
559
560   /* Normalize the divisor.  */
561   bit = precision - APInt::tcMSB(divisor, partsCount) - 1;
562   if(bit) {
563     exponent += bit;
564     APInt::tcShiftLeft(divisor, partsCount, bit);
565   }
566
567   /* Normalize the dividend.  */
568   bit = precision - APInt::tcMSB(dividend, partsCount) - 1;
569   if(bit) {
570     exponent -= bit;
571     APInt::tcShiftLeft(dividend, partsCount, bit);
572   }
573
574   if(APInt::tcCompare(dividend, divisor, partsCount) < 0) {
575     exponent--;
576     APInt::tcShiftLeft(dividend, partsCount, 1);
577     assert(APInt::tcCompare(dividend, divisor, partsCount) >= 0);
578   }
579
580   /* Long division.  */
581   for(bit = precision; bit; bit -= 1) {
582     if(APInt::tcCompare(dividend, divisor, partsCount) >= 0) {
583       APInt::tcSubtract(dividend, divisor, 0, partsCount);
584       APInt::tcSetBit(lhsSignificand, bit - 1);
585     }
586
587     APInt::tcShiftLeft(dividend, partsCount, 1);
588   }
589
590   /* Figure out the lost fraction.  */
591   int cmp = APInt::tcCompare(dividend, divisor, partsCount);
592
593   if(cmp > 0)
594     lost_fraction = lfMoreThanHalf;
595   else if(cmp == 0)
596     lost_fraction = lfExactlyHalf;
597   else if(APInt::tcIsZero(dividend, partsCount))
598     lost_fraction = lfExactlyZero;
599   else
600     lost_fraction = lfLessThanHalf;
601
602   if(partsCount > 2)
603     delete [] dividend;
604
605   return lost_fraction;
606 }
607
608 unsigned int
609 APFloat::significandMSB() const
610 {
611   return APInt::tcMSB(significandParts(), partCount());
612 }
613
614 unsigned int
615 APFloat::significandLSB() const
616 {
617   return APInt::tcLSB(significandParts(), partCount());
618 }
619
620 /* Note that a zero result is NOT normalized to fcZero.  */
621 lostFraction
622 APFloat::shiftSignificandRight(unsigned int bits)
623 {
624   /* Our exponent should not overflow.  */
625   assert((exponent_t) (exponent + bits) >= exponent);
626
627   exponent += bits;
628
629   return shiftRight(significandParts(), partCount(), bits);
630 }
631
632 /* Shift the significand left BITS bits, subtract BITS from its exponent.  */
633 void
634 APFloat::shiftSignificandLeft(unsigned int bits)
635 {
636   assert(bits < semantics->precision);
637
638   if(bits) {
639     unsigned int partsCount = partCount();
640
641     APInt::tcShiftLeft(significandParts(), partsCount, bits);
642     exponent -= bits;
643
644     assert(!APInt::tcIsZero(significandParts(), partsCount));
645   }
646 }
647
648 APFloat::cmpResult
649 APFloat::compareAbsoluteValue(const APFloat &rhs) const
650 {
651   int compare;
652
653   assert(semantics == rhs.semantics);
654   assert(category == fcNormal);
655   assert(rhs.category == fcNormal);
656
657   compare = exponent - rhs.exponent;
658
659   /* If exponents are equal, do an unsigned bignum comparison of the
660      significands.  */
661   if(compare == 0)
662     compare = APInt::tcCompare(significandParts(), rhs.significandParts(),
663                                partCount());
664
665   if(compare > 0)
666     return cmpGreaterThan;
667   else if(compare < 0)
668     return cmpLessThan;
669   else
670     return cmpEqual;
671 }
672
673 /* Handle overflow.  Sign is preserved.  We either become infinity or
674    the largest finite number.  */
675 APFloat::opStatus
676 APFloat::handleOverflow(roundingMode rounding_mode)
677 {
678   /* Infinity?  */
679   if(rounding_mode == rmNearestTiesToEven
680      || rounding_mode == rmNearestTiesToAway
681      || (rounding_mode == rmTowardPositive && !sign)
682      || (rounding_mode == rmTowardNegative && sign))
683     {
684       category = fcInfinity;
685       return (opStatus) (opOverflow | opInexact);
686     }
687
688   /* Otherwise we become the largest finite number.  */
689   category = fcNormal;
690   exponent = semantics->maxExponent;
691   APInt::tcSetLeastSignificantBits(significandParts(), partCount(),
692                                    semantics->precision);
693
694   return opInexact;
695 }
696
697 /* This routine must work for fcZero of both signs, and fcNormal
698    numbers.  */
699 bool
700 APFloat::roundAwayFromZero(roundingMode rounding_mode,
701                            lostFraction lost_fraction)
702 {
703   /* NaNs and infinities should not have lost fractions.  */
704   assert(category == fcNormal || category == fcZero);
705
706   /* Our caller has already handled this case.  */
707   assert(lost_fraction != lfExactlyZero);
708
709   switch(rounding_mode) {
710   default:
711     assert(0);
712
713   case rmNearestTiesToAway:
714     return lost_fraction == lfExactlyHalf || lost_fraction == lfMoreThanHalf;
715
716   case rmNearestTiesToEven:
717     if(lost_fraction == lfMoreThanHalf)
718       return true;
719
720     /* Our zeroes don't have a significand to test.  */
721     if(lost_fraction == lfExactlyHalf && category != fcZero)
722       return significandParts()[0] & 1;
723
724     return false;
725
726   case rmTowardZero:
727     return false;
728
729   case rmTowardPositive:
730     return sign == false;
731
732   case rmTowardNegative:
733     return sign == true;
734   }
735 }
736
737 APFloat::opStatus
738 APFloat::normalize(roundingMode rounding_mode,
739                    lostFraction lost_fraction)
740 {
741   unsigned int omsb;            /* One, not zero, based MSB.  */
742   int exponentChange;
743
744   if(category != fcNormal)
745     return opOK;
746
747   /* Before rounding normalize the exponent of fcNormal numbers.  */
748   omsb = significandMSB() + 1;
749
750   if(omsb) {
751     /* OMSB is numbered from 1.  We want to place it in the integer
752        bit numbered PRECISON if possible, with a compensating change in
753        the exponent.  */
754     exponentChange = omsb - semantics->precision;
755
756     /* If the resulting exponent is too high, overflow according to
757        the rounding mode.  */
758     if(exponent + exponentChange > semantics->maxExponent)
759       return handleOverflow(rounding_mode);
760
761     /* Subnormal numbers have exponent minExponent, and their MSB
762        is forced based on that.  */
763     if(exponent + exponentChange < semantics->minExponent)
764       exponentChange = semantics->minExponent - exponent;
765
766     /* Shifting left is easy as we don't lose precision.  */
767     if(exponentChange < 0) {
768       assert(lost_fraction == lfExactlyZero);
769
770       shiftSignificandLeft(-exponentChange);
771
772       return opOK;
773     }
774
775     if(exponentChange > 0) {
776       lostFraction lf;
777
778       /* Shift right and capture any new lost fraction.  */
779       lf = shiftSignificandRight(exponentChange);
780
781       lost_fraction = combineLostFractions(lf, lost_fraction);
782
783       /* Keep OMSB up-to-date.  */
784       if(omsb > (unsigned) exponentChange)
785         omsb -= (unsigned) exponentChange;
786       else
787         omsb = 0;
788     }
789   }
790
791   /* Now round the number according to rounding_mode given the lost
792      fraction.  */
793
794   /* As specified in IEEE 754, since we do not trap we do not report
795      underflow for exact results.  */
796   if(lost_fraction == lfExactlyZero) {
797     /* Canonicalize zeroes.  */
798     if(omsb == 0)
799       category = fcZero;
800
801     return opOK;
802   }
803
804   /* Increment the significand if we're rounding away from zero.  */
805   if(roundAwayFromZero(rounding_mode, lost_fraction)) {
806     if(omsb == 0)
807       exponent = semantics->minExponent;
808
809     incrementSignificand();
810     omsb = significandMSB() + 1;
811
812     /* Did the significand increment overflow?  */
813     if(omsb == (unsigned) semantics->precision + 1) {
814       /* Renormalize by incrementing the exponent and shifting our
815          significand right one.  However if we already have the
816          maximum exponent we overflow to infinity.  */
817       if(exponent == semantics->maxExponent) {
818         category = fcInfinity;
819
820         return (opStatus) (opOverflow | opInexact);
821       }
822
823       shiftSignificandRight(1);
824
825       return opInexact;
826     }
827   }
828
829   /* The normal case - we were and are not denormal, and any
830      significand increment above didn't overflow.  */
831   if(omsb == semantics->precision)
832     return opInexact;
833
834   /* We have a non-zero denormal.  */
835   assert(omsb < semantics->precision);
836   assert(exponent == semantics->minExponent);
837
838   /* Canonicalize zeroes.  */
839   if(omsb == 0)
840     category = fcZero;
841
842   /* The fcZero case is a denormal that underflowed to zero.  */
843   return (opStatus) (opUnderflow | opInexact);
844 }
845
846 APFloat::opStatus
847 APFloat::addOrSubtractSpecials(const APFloat &rhs, bool subtract)
848 {
849   switch(convolve(category, rhs.category)) {
850   default:
851     assert(0);
852
853   case convolve(fcNaN, fcZero):
854   case convolve(fcNaN, fcNormal):
855   case convolve(fcNaN, fcInfinity):
856   case convolve(fcNaN, fcNaN):
857   case convolve(fcNormal, fcZero):
858   case convolve(fcInfinity, fcNormal):
859   case convolve(fcInfinity, fcZero):
860     return opOK;
861
862   case convolve(fcZero, fcNaN):
863   case convolve(fcNormal, fcNaN):
864   case convolve(fcInfinity, fcNaN):
865     category = fcNaN;
866     copySignificand(rhs);
867     return opOK;
868
869   case convolve(fcNormal, fcInfinity):
870   case convolve(fcZero, fcInfinity):
871     category = fcInfinity;
872     sign = rhs.sign ^ subtract;
873     return opOK;
874
875   case convolve(fcZero, fcNormal):
876     assign(rhs);
877     sign = rhs.sign ^ subtract;
878     return opOK;
879
880   case convolve(fcZero, fcZero):
881     /* Sign depends on rounding mode; handled by caller.  */
882     return opOK;
883
884   case convolve(fcInfinity, fcInfinity):
885     /* Differently signed infinities can only be validly
886        subtracted.  */
887     if(sign ^ rhs.sign != subtract) {
888       category = fcNaN;
889       // Arbitrary but deterministic value for significand
890       APInt::tcSet(significandParts(), ~0U, partCount());
891       return opInvalidOp;
892     }
893
894     return opOK;
895
896   case convolve(fcNormal, fcNormal):
897     return opDivByZero;
898   }
899 }
900
901 /* Add or subtract two normal numbers.  */
902 lostFraction
903 APFloat::addOrSubtractSignificand(const APFloat &rhs, bool subtract)
904 {
905   integerPart carry;
906   lostFraction lost_fraction;
907   int bits;
908
909   /* Determine if the operation on the absolute values is effectively
910      an addition or subtraction.  */
911   subtract ^= (sign ^ rhs.sign);
912
913   /* Are we bigger exponent-wise than the RHS?  */
914   bits = exponent - rhs.exponent;
915
916   /* Subtraction is more subtle than one might naively expect.  */
917   if(subtract) {
918     APFloat temp_rhs(rhs);
919     bool reverse;
920
921     if (bits == 0) {
922       reverse = compareAbsoluteValue(temp_rhs) == cmpLessThan;
923       lost_fraction = lfExactlyZero;
924     } else if (bits > 0) {
925       lost_fraction = temp_rhs.shiftSignificandRight(bits - 1);
926       shiftSignificandLeft(1);
927       reverse = false;
928     } else {
929       lost_fraction = shiftSignificandRight(-bits - 1);
930       temp_rhs.shiftSignificandLeft(1);
931       reverse = true;
932     }
933
934     if (reverse) {
935       carry = temp_rhs.subtractSignificand
936         (*this, lost_fraction != lfExactlyZero);
937       copySignificand(temp_rhs);
938       sign = !sign;
939     } else {
940       carry = subtractSignificand
941         (temp_rhs, lost_fraction != lfExactlyZero);
942     }
943
944     /* Invert the lost fraction - it was on the RHS and
945        subtracted.  */
946     if(lost_fraction == lfLessThanHalf)
947       lost_fraction = lfMoreThanHalf;
948     else if(lost_fraction == lfMoreThanHalf)
949       lost_fraction = lfLessThanHalf;
950
951     /* The code above is intended to ensure that no borrow is
952        necessary.  */
953     assert(!carry);
954   } else {
955     if(bits > 0) {
956       APFloat temp_rhs(rhs);
957
958       lost_fraction = temp_rhs.shiftSignificandRight(bits);
959       carry = addSignificand(temp_rhs);
960     } else {
961       lost_fraction = shiftSignificandRight(-bits);
962       carry = addSignificand(rhs);
963     }
964
965     /* We have a guard bit; generating a carry cannot happen.  */
966     assert(!carry);
967   }
968
969   return lost_fraction;
970 }
971
972 APFloat::opStatus
973 APFloat::multiplySpecials(const APFloat &rhs)
974 {
975   switch(convolve(category, rhs.category)) {
976   default:
977     assert(0);
978
979   case convolve(fcNaN, fcZero):
980   case convolve(fcNaN, fcNormal):
981   case convolve(fcNaN, fcInfinity):
982   case convolve(fcNaN, fcNaN):
983     return opOK;
984
985   case convolve(fcZero, fcNaN):
986   case convolve(fcNormal, fcNaN):
987   case convolve(fcInfinity, fcNaN):
988     category = fcNaN;
989     copySignificand(rhs);
990     return opOK;
991
992   case convolve(fcNormal, fcInfinity):
993   case convolve(fcInfinity, fcNormal):
994   case convolve(fcInfinity, fcInfinity):
995     category = fcInfinity;
996     return opOK;
997
998   case convolve(fcZero, fcNormal):
999   case convolve(fcNormal, fcZero):
1000   case convolve(fcZero, fcZero):
1001     category = fcZero;
1002     return opOK;
1003
1004   case convolve(fcZero, fcInfinity):
1005   case convolve(fcInfinity, fcZero):
1006     category = fcNaN;
1007     // Arbitrary but deterministic value for significand
1008     APInt::tcSet(significandParts(), ~0U, partCount());
1009     return opInvalidOp;
1010
1011   case convolve(fcNormal, fcNormal):
1012     return opOK;
1013   }
1014 }
1015
1016 APFloat::opStatus
1017 APFloat::divideSpecials(const APFloat &rhs)
1018 {
1019   switch(convolve(category, rhs.category)) {
1020   default:
1021     assert(0);
1022
1023   case convolve(fcNaN, fcZero):
1024   case convolve(fcNaN, fcNormal):
1025   case convolve(fcNaN, fcInfinity):
1026   case convolve(fcNaN, fcNaN):
1027   case convolve(fcInfinity, fcZero):
1028   case convolve(fcInfinity, fcNormal):
1029   case convolve(fcZero, fcInfinity):
1030   case convolve(fcZero, fcNormal):
1031     return opOK;
1032
1033   case convolve(fcZero, fcNaN):
1034   case convolve(fcNormal, fcNaN):
1035   case convolve(fcInfinity, fcNaN):
1036     category = fcNaN;
1037     copySignificand(rhs);
1038     return opOK;
1039
1040   case convolve(fcNormal, fcInfinity):
1041     category = fcZero;
1042     return opOK;
1043
1044   case convolve(fcNormal, fcZero):
1045     category = fcInfinity;
1046     return opDivByZero;
1047
1048   case convolve(fcInfinity, fcInfinity):
1049   case convolve(fcZero, fcZero):
1050     category = fcNaN;
1051     // Arbitrary but deterministic value for significand
1052     APInt::tcSet(significandParts(), ~0U, partCount());
1053     return opInvalidOp;
1054
1055   case convolve(fcNormal, fcNormal):
1056     return opOK;
1057   }
1058 }
1059
1060 /* Change sign.  */
1061 void
1062 APFloat::changeSign()
1063 {
1064   /* Look mummy, this one's easy.  */
1065   sign = !sign;
1066 }
1067
1068 void
1069 APFloat::clearSign()
1070 {
1071   /* So is this one. */
1072   sign = 0;
1073 }
1074
1075 void
1076 APFloat::copySign(const APFloat &rhs)
1077 {
1078   /* And this one. */
1079   sign = rhs.sign;
1080 }
1081
1082 /* Normalized addition or subtraction.  */
1083 APFloat::opStatus
1084 APFloat::addOrSubtract(const APFloat &rhs, roundingMode rounding_mode,
1085                        bool subtract)
1086 {
1087   opStatus fs;
1088
1089   fs = addOrSubtractSpecials(rhs, subtract);
1090
1091   /* This return code means it was not a simple case.  */
1092   if(fs == opDivByZero) {
1093     lostFraction lost_fraction;
1094
1095     lost_fraction = addOrSubtractSignificand(rhs, subtract);
1096     fs = normalize(rounding_mode, lost_fraction);
1097
1098     /* Can only be zero if we lost no fraction.  */
1099     assert(category != fcZero || lost_fraction == lfExactlyZero);
1100   }
1101
1102   /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1103      positive zero unless rounding to minus infinity, except that
1104      adding two like-signed zeroes gives that zero.  */
1105   if(category == fcZero) {
1106     if(rhs.category != fcZero || (sign == rhs.sign) == subtract)
1107       sign = (rounding_mode == rmTowardNegative);
1108   }
1109
1110   return fs;
1111 }
1112
1113 /* Normalized addition.  */
1114 APFloat::opStatus
1115 APFloat::add(const APFloat &rhs, roundingMode rounding_mode)
1116 {
1117   return addOrSubtract(rhs, rounding_mode, false);
1118 }
1119
1120 /* Normalized subtraction.  */
1121 APFloat::opStatus
1122 APFloat::subtract(const APFloat &rhs, roundingMode rounding_mode)
1123 {
1124   return addOrSubtract(rhs, rounding_mode, true);
1125 }
1126
1127 /* Normalized multiply.  */
1128 APFloat::opStatus
1129 APFloat::multiply(const APFloat &rhs, roundingMode rounding_mode)
1130 {
1131   opStatus fs;
1132
1133   sign ^= rhs.sign;
1134   fs = multiplySpecials(rhs);
1135
1136   if(category == fcNormal) {
1137     lostFraction lost_fraction = multiplySignificand(rhs, 0);
1138     fs = normalize(rounding_mode, lost_fraction);
1139     if(lost_fraction != lfExactlyZero)
1140       fs = (opStatus) (fs | opInexact);
1141   }
1142
1143   return fs;
1144 }
1145
1146 /* Normalized divide.  */
1147 APFloat::opStatus
1148 APFloat::divide(const APFloat &rhs, roundingMode rounding_mode)
1149 {
1150   opStatus fs;
1151
1152   sign ^= rhs.sign;
1153   fs = divideSpecials(rhs);
1154
1155   if(category == fcNormal) {
1156     lostFraction lost_fraction = divideSignificand(rhs);
1157     fs = normalize(rounding_mode, lost_fraction);
1158     if(lost_fraction != lfExactlyZero)
1159       fs = (opStatus) (fs | opInexact);
1160   }
1161
1162   return fs;
1163 }
1164
1165 /* Normalized remainder. */
1166 APFloat::opStatus
1167 APFloat::mod(const APFloat &rhs, roundingMode rounding_mode)
1168 {
1169   opStatus fs;
1170   APFloat V = *this;
1171   fs = V.divide(rhs, rmNearestTiesToEven);
1172   if (fs == opDivByZero)
1173     return fs;
1174
1175   integerPart x;
1176   fs = V.convertToInteger(&x, integerPartWidth, true, rmNearestTiesToEven);
1177   if (fs==opInvalidOp)
1178     return fs;
1179
1180   fs = V.convertFromInteger(&x, integerPartWidth, true, rmNearestTiesToEven);
1181   assert(fs==opOK);   // should always work
1182   fs = V.multiply(rhs, rounding_mode);
1183   assert(fs==opOK);   // should not overflow or underflow
1184   fs = subtract(V, rounding_mode);
1185   assert(fs==opOK);
1186   return fs;
1187 }
1188
1189 /* Normalized fused-multiply-add.  */
1190 APFloat::opStatus
1191 APFloat::fusedMultiplyAdd(const APFloat &multiplicand,
1192                           const APFloat &addend,
1193                           roundingMode rounding_mode)
1194 {
1195   opStatus fs;
1196
1197   /* Post-multiplication sign, before addition.  */
1198   sign ^= multiplicand.sign;
1199
1200   /* If and only if all arguments are normal do we need to do an
1201      extended-precision calculation.  */
1202   if(category == fcNormal
1203      && multiplicand.category == fcNormal
1204      && addend.category == fcNormal) {
1205     lostFraction lost_fraction;
1206
1207     lost_fraction = multiplySignificand(multiplicand, &addend);
1208     fs = normalize(rounding_mode, lost_fraction);
1209     if(lost_fraction != lfExactlyZero)
1210       fs = (opStatus) (fs | opInexact);
1211
1212     /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1213        positive zero unless rounding to minus infinity, except that
1214        adding two like-signed zeroes gives that zero.  */
1215     if(category == fcZero && sign != addend.sign)
1216       sign = (rounding_mode == rmTowardNegative);
1217   } else {
1218     fs = multiplySpecials(multiplicand);
1219
1220     /* FS can only be opOK or opInvalidOp.  There is no more work
1221        to do in the latter case.  The IEEE-754R standard says it is
1222        implementation-defined in this case whether, if ADDEND is a
1223        quiet NaN, we raise invalid op; this implementation does so.
1224
1225        If we need to do the addition we can do so with normal
1226        precision.  */
1227     if(fs == opOK)
1228       fs = addOrSubtract(addend, rounding_mode, false);
1229   }
1230
1231   return fs;
1232 }
1233
1234 /* Comparison requires normalized numbers.  */
1235 APFloat::cmpResult
1236 APFloat::compare(const APFloat &rhs) const
1237 {
1238   cmpResult result;
1239
1240   assert(semantics == rhs.semantics);
1241
1242   switch(convolve(category, rhs.category)) {
1243   default:
1244     assert(0);
1245
1246   case convolve(fcNaN, fcZero):
1247   case convolve(fcNaN, fcNormal):
1248   case convolve(fcNaN, fcInfinity):
1249   case convolve(fcNaN, fcNaN):
1250   case convolve(fcZero, fcNaN):
1251   case convolve(fcNormal, fcNaN):
1252   case convolve(fcInfinity, fcNaN):
1253     return cmpUnordered;
1254
1255   case convolve(fcInfinity, fcNormal):
1256   case convolve(fcInfinity, fcZero):
1257   case convolve(fcNormal, fcZero):
1258     if(sign)
1259       return cmpLessThan;
1260     else
1261       return cmpGreaterThan;
1262
1263   case convolve(fcNormal, fcInfinity):
1264   case convolve(fcZero, fcInfinity):
1265   case convolve(fcZero, fcNormal):
1266     if(rhs.sign)
1267       return cmpGreaterThan;
1268     else
1269       return cmpLessThan;
1270
1271   case convolve(fcInfinity, fcInfinity):
1272     if(sign == rhs.sign)
1273       return cmpEqual;
1274     else if(sign)
1275       return cmpLessThan;
1276     else
1277       return cmpGreaterThan;
1278
1279   case convolve(fcZero, fcZero):
1280     return cmpEqual;
1281
1282   case convolve(fcNormal, fcNormal):
1283     break;
1284   }
1285
1286   /* Two normal numbers.  Do they have the same sign?  */
1287   if(sign != rhs.sign) {
1288     if(sign)
1289       result = cmpLessThan;
1290     else
1291       result = cmpGreaterThan;
1292   } else {
1293     /* Compare absolute values; invert result if negative.  */
1294     result = compareAbsoluteValue(rhs);
1295
1296     if(sign) {
1297       if(result == cmpLessThan)
1298         result = cmpGreaterThan;
1299       else if(result == cmpGreaterThan)
1300         result = cmpLessThan;
1301     }
1302   }
1303
1304   return result;
1305 }
1306
1307 APFloat::opStatus
1308 APFloat::convert(const fltSemantics &toSemantics,
1309                  roundingMode rounding_mode)
1310 {
1311   unsigned int newPartCount;
1312   opStatus fs;
1313
1314   newPartCount = partCountForBits(toSemantics.precision + 1);
1315
1316   /* If our new form is wider, re-allocate our bit pattern into wider
1317      storage.  */
1318   if(newPartCount > partCount()) {
1319     integerPart *newParts;
1320
1321     newParts = new integerPart[newPartCount];
1322     APInt::tcSet(newParts, 0, newPartCount);
1323     APInt::tcAssign(newParts, significandParts(), partCount());
1324     freeSignificand();
1325     significand.parts = newParts;
1326   }
1327
1328   if(category == fcNormal) {
1329     /* Re-interpret our bit-pattern.  */
1330     exponent += toSemantics.precision - semantics->precision;
1331     semantics = &toSemantics;
1332     fs = normalize(rounding_mode, lfExactlyZero);
1333   } else {
1334     semantics = &toSemantics;
1335     fs = opOK;
1336   }
1337
1338   return fs;
1339 }
1340
1341 /* Convert a floating point number to an integer according to the
1342    rounding mode.  If the rounded integer value is out of range this
1343    returns an invalid operation exception.  If the rounded value is in
1344    range but the floating point number is not the exact integer, the C
1345    standard doesn't require an inexact exception to be raised.  IEEE
1346    854 does require it so we do that.
1347
1348    Note that for conversions to integer type the C standard requires
1349    round-to-zero to always be used.  */
1350 APFloat::opStatus
1351 APFloat::convertToInteger(integerPart *parts, unsigned int width,
1352                           bool isSigned,
1353                           roundingMode rounding_mode) const
1354 {
1355   lostFraction lost_fraction;
1356   unsigned int msb, partsCount;
1357   int bits;
1358
1359   /* Handle the three special cases first.  */
1360   if(category == fcInfinity || category == fcNaN)
1361     return opInvalidOp;
1362
1363   partsCount = partCountForBits(width);
1364
1365   if(category == fcZero) {
1366     APInt::tcSet(parts, 0, partsCount);
1367     return opOK;
1368   }
1369
1370   /* Shift the bit pattern so the fraction is lost.  */
1371   APFloat tmp(*this);
1372
1373   bits = (int) semantics->precision - 1 - exponent;
1374
1375   if(bits > 0) {
1376     lost_fraction = tmp.shiftSignificandRight(bits);
1377   } else {
1378     tmp.shiftSignificandLeft(-bits);
1379     lost_fraction = lfExactlyZero;
1380   }
1381
1382   if(lost_fraction != lfExactlyZero
1383      && tmp.roundAwayFromZero(rounding_mode, lost_fraction))
1384     tmp.incrementSignificand();
1385
1386   msb = tmp.significandMSB();
1387
1388   /* Negative numbers cannot be represented as unsigned.  */
1389   if(!isSigned && tmp.sign && msb != -1U)
1390     return opInvalidOp;
1391
1392   /* It takes exponent + 1 bits to represent the truncated floating
1393      point number without its sign.  We lose a bit for the sign, but
1394      the maximally negative integer is a special case.  */
1395   if(msb + 1 > width)           /* !! Not same as msb >= width !! */
1396     return opInvalidOp;
1397
1398   if(isSigned && msb + 1 == width
1399      && (!tmp.sign || tmp.significandLSB() != msb))
1400     return opInvalidOp;
1401
1402   APInt::tcAssign(parts, tmp.significandParts(), partsCount);
1403
1404   if(tmp.sign)
1405     APInt::tcNegate(parts, partsCount);
1406
1407   if(lost_fraction == lfExactlyZero)
1408     return opOK;
1409   else
1410     return opInexact;
1411 }
1412
1413 APFloat::opStatus
1414 APFloat::convertFromUnsignedInteger(integerPart *parts,
1415                                     unsigned int partCount,
1416                                     roundingMode rounding_mode)
1417 {
1418   unsigned int msb, precision;
1419   lostFraction lost_fraction;
1420
1421   msb = APInt::tcMSB(parts, partCount) + 1;
1422   precision = semantics->precision;
1423
1424   category = fcNormal;
1425   exponent = precision - 1;
1426
1427   if(msb > precision) {
1428     exponent += (msb - precision);
1429     lost_fraction = shiftRight(parts, partCount, msb - precision);
1430     msb = precision;
1431   } else
1432     lost_fraction = lfExactlyZero;
1433
1434   /* Copy the bit image.  */
1435   zeroSignificand();
1436   APInt::tcAssign(significandParts(), parts, partCountForBits(msb));
1437
1438   return normalize(rounding_mode, lost_fraction);
1439 }
1440
1441 APFloat::opStatus
1442 APFloat::convertFromInteger(const integerPart *parts,
1443                             unsigned int partCount, bool isSigned,
1444                             roundingMode rounding_mode)
1445 {
1446   unsigned int width;
1447   opStatus status;
1448   integerPart *copy;
1449
1450   copy = new integerPart[partCount];
1451   APInt::tcAssign(copy, parts, partCount);
1452
1453   width = partCount * integerPartWidth;
1454
1455   sign = false;
1456   if(isSigned && APInt::tcExtractBit(parts, width - 1)) {
1457     sign = true;
1458     APInt::tcNegate(copy, partCount);
1459   }
1460
1461   status = convertFromUnsignedInteger(copy, partCount, rounding_mode);
1462   delete [] copy;
1463
1464   return status;
1465 }
1466
1467 APFloat::opStatus
1468 APFloat::convertFromHexadecimalString(const char *p,
1469                                       roundingMode rounding_mode)
1470 {
1471   lostFraction lost_fraction;
1472   integerPart *significand;
1473   unsigned int bitPos, partsCount;
1474   const char *dot, *firstSignificantDigit;
1475
1476   zeroSignificand();
1477   exponent = 0;
1478   category = fcNormal;
1479
1480   significand = significandParts();
1481   partsCount = partCount();
1482   bitPos = partsCount * integerPartWidth;
1483
1484   /* Skip leading zeroes and any(hexa)decimal point.  */
1485   p = skipLeadingZeroesAndAnyDot(p, &dot);
1486   firstSignificantDigit = p;
1487
1488   for(;;) {
1489     integerPart hex_value;
1490
1491     if(*p == '.') {
1492       assert(dot == 0);
1493       dot = p++;
1494     }
1495
1496     hex_value = hexDigitValue(*p);
1497     if(hex_value == -1U) {
1498       lost_fraction = lfExactlyZero;
1499       break;
1500     }
1501
1502     p++;
1503
1504     /* Store the number whilst 4-bit nibbles remain.  */
1505     if(bitPos) {
1506       bitPos -= 4;
1507       hex_value <<= bitPos % integerPartWidth;
1508       significand[bitPos / integerPartWidth] |= hex_value;
1509     } else {
1510       lost_fraction = trailingHexadecimalFraction(p, hex_value);
1511       while(hexDigitValue(*p) != -1U)
1512         p++;
1513       break;
1514     }
1515   }
1516
1517   /* Hex floats require an exponent but not a hexadecimal point.  */
1518   assert(*p == 'p' || *p == 'P');
1519
1520   /* Ignore the exponent if we are zero.  */
1521   if(p != firstSignificantDigit) {
1522     int expAdjustment;
1523
1524     /* Implicit hexadecimal point?  */
1525     if(!dot)
1526       dot = p;
1527
1528     /* Calculate the exponent adjustment implicit in the number of
1529        significant digits.  */
1530     expAdjustment = dot - firstSignificantDigit;
1531     if(expAdjustment < 0)
1532       expAdjustment++;
1533     expAdjustment = expAdjustment * 4 - 1;
1534
1535     /* Adjust for writing the significand starting at the most
1536        significant nibble.  */
1537     expAdjustment += semantics->precision;
1538     expAdjustment -= partsCount * integerPartWidth;
1539
1540     /* Adjust for the given exponent.  */
1541     exponent = totalExponent(p, expAdjustment);
1542   }
1543
1544   return normalize(rounding_mode, lost_fraction);
1545 }
1546
1547 APFloat::opStatus
1548 APFloat::convertFromString(const char *p, roundingMode rounding_mode) {
1549   /* Handle a leading minus sign.  */
1550   if(*p == '-')
1551     sign = 1, p++;
1552   else
1553     sign = 0;
1554
1555   if(p[0] == '0' && (p[1] == 'x' || p[1] == 'X'))
1556     return convertFromHexadecimalString(p + 2, rounding_mode);
1557
1558   assert(0 && "Decimal to binary conversions not yet implemented");
1559   abort();
1560 }
1561
1562 // For good performance it is desirable for different APFloats
1563 // to produce different integers.
1564 uint32_t
1565 APFloat::getHashValue() const { 
1566   if (category==fcZero) return sign<<8 | semantics->precision ;
1567   else if (category==fcInfinity) return sign<<9 | semantics->precision;
1568   else if (category==fcNaN) return 1<<10 | semantics->precision;
1569   else {
1570     uint32_t hash = sign<<11 | semantics->precision | exponent<<12;
1571     const integerPart* p = significandParts();
1572     for (int i=partCount(); i>0; i--, p++)
1573       hash ^= ((uint32_t)*p) ^ (*p)>>32;
1574     return hash;
1575   }
1576 }
1577
1578 // Conversion from APFloat to/from host float/double.  It may eventually be
1579 // possible to eliminate these and have everybody deal with APFloats, but that
1580 // will take a while.  This approach will not easily extend to long double.
1581 // Current implementation requires partCount()==1, which is correct at the
1582 // moment but could be made more general.
1583
1584 double
1585 APFloat::convertToDouble() const {
1586   assert(semantics == (const llvm::fltSemantics* const)&IEEEdouble);
1587   assert (partCount()==1);
1588
1589   uint64_t myexponent, mysignificand;
1590
1591   if (category==fcNormal) {
1592     mysignificand = *significandParts();
1593     myexponent = exponent+1023; //bias
1594   } else if (category==fcZero) {
1595     myexponent = 0;
1596     mysignificand = 0;
1597   } else if (category==fcInfinity) {
1598     myexponent = 0x7ff;
1599     mysignificand = 0;
1600   } else if (category==fcNaN) {
1601     myexponent = 0x7ff;
1602     mysignificand = *significandParts();
1603   } else
1604     assert(0);
1605
1606   return BitsToDouble((((uint64_t)sign & 1) << 63) | 
1607         ((myexponent & 0x7ff) <<  52) | 
1608         (mysignificand & 0xfffffffffffffLL));
1609 }
1610
1611 float
1612 APFloat::convertToFloat() const {
1613   assert(semantics == (const llvm::fltSemantics* const)&IEEEsingle);
1614   assert (partCount()==1);
1615
1616   uint32_t myexponent, mysignificand;
1617
1618   if (category==fcNormal) {
1619     myexponent = exponent+127; //bias
1620     mysignificand = *significandParts();
1621   } else if (category==fcZero) {
1622     myexponent = 0;
1623     mysignificand = 0;
1624   } else if (category==fcInfinity) {
1625     myexponent = 0xff;
1626     mysignificand = 0;
1627   } else if (category==fcNaN) {
1628     myexponent = 0x7ff;
1629     mysignificand = *significandParts();
1630   } else
1631     assert(0);
1632
1633   return BitsToFloat(((sign&1) << 31) | ((myexponent&0xff) << 23) | 
1634         (mysignificand & 0x7fffff));
1635 }
1636
1637 APFloat::APFloat(double d) {
1638   uint64_t i = DoubleToBits(d);
1639   uint64_t myexponent = (i >> 52) & 0x7ff;
1640   uint64_t mysignificand = i & 0xfffffffffffffLL;
1641
1642   initialize(&APFloat::IEEEdouble);
1643   assert(partCount()==1);
1644
1645   sign = i>>63;
1646   if (myexponent==0 && mysignificand==0) {
1647     // exponent, significand meaningless
1648     category = fcZero;
1649   } else if (myexponent==0x7ff && mysignificand==0) {
1650     // exponent, significand meaningless
1651     category = fcInfinity;
1652   } else if (myexponent==0x7ff && mysignificand!=0) {
1653     // exponent meaningless
1654     category = fcNaN;
1655     *significandParts() = mysignificand;
1656   } else {
1657     category = fcNormal;
1658     exponent = myexponent - 1023;
1659     *significandParts() = mysignificand | 0x10000000000000LL; 
1660  }
1661 }
1662
1663 APFloat::APFloat(float f) {
1664   uint32_t i = FloatToBits(f);
1665   uint32_t myexponent = (i >> 23) & 0xff;
1666   uint32_t mysignificand = i & 0x7fffff;
1667
1668   initialize(&APFloat::IEEEsingle);
1669   assert(partCount()==1);
1670
1671   sign = i >> 31;
1672   if (myexponent==0 && mysignificand==0) {
1673     // exponent, significand meaningless
1674     category = fcZero;
1675   } else if (myexponent==0xff && mysignificand==0) {
1676     // exponent, significand meaningless
1677     category = fcInfinity;
1678   } else if (myexponent==0xff && (mysignificand & 0x400000)) {
1679     // sign, exponent, significand meaningless
1680     category = fcNaN;
1681     *significandParts() = mysignificand;
1682   } else {
1683     category = fcNormal;
1684     exponent = myexponent - 127;  //bias
1685     *significandParts() = mysignificand | 0x800000; // integer bit
1686   }
1687 }