APFloat.cpp revision e15c2db9935eee66a8008f1bd09882aff2ed3aae 
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
19using namespace llvm;
20
21#define convolve(lhs, rhs) ((lhs) * 4 + (rhs))
22
23/* Assumed in hexadecimal significand parsing.  */
24COMPILE_TIME_ASSERT(integerPartWidth % 4 == 0);
25
26namespace 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.  */
54namespace {
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.  */
224void
225APFloat::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
235void
236APFloat::freeSignificand()
237{
238  if(partCount() > 1)
239    delete [] significand.parts;
240}
241
242void
243APFloat::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
254void
255APFloat::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
264APFloat &
265APFloat::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
278bool
279APFloat::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
302APFloat::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
312APFloat::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
322APFloat::APFloat(const fltSemantics &ourSemantics, const char *text)
323{
324  initialize(&ourSemantics);
325  convertFromString(text, rmNearestTiesToEven);
326}
327
328APFloat::APFloat(const APFloat &rhs)
329{
330  initialize(rhs.semantics);
331  assign(rhs);
332}
333
334APFloat::~APFloat()
335{
336  freeSignificand();
337}
338
339unsigned int
340APFloat::partCount() const
341{
342  return partCountForBits(semantics->precision + 1);
343}
344
345unsigned int
346APFloat::semanticsPrecision(const fltSemantics &semantics)
347{
348  return semantics.precision;
349}
350
351const integerPart *
352APFloat::significandParts() const
353{
354  return const_cast<APFloat *>(this)->significandParts();
355}
356
357integerPart *
358APFloat::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.  */
369lostFraction
370APFloat::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
383void
384APFloat::zeroSignificand()
385{
386  category = fcNormal;
387  APInt::tcSet(significandParts(), 0, partCount());
388}
389
390/* Increment an fcNormal floating point number's significand.  */
391void
392APFloat::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.  */
403integerPart
404APFloat::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.  */
418integerPart
419APFloat::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.  */
435lostFraction
436APFloat::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.  */
527lostFraction
528APFloat::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
608unsigned int
609APFloat::significandMSB() const
610{
611  return APInt::tcMSB(significandParts(), partCount());
612}
613
614unsigned int
615APFloat::significandLSB() const
616{
617  return APInt::tcLSB(significandParts(), partCount());
618}
619
620/* Note that a zero result is NOT normalized to fcZero.  */
621lostFraction
622APFloat::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.  */
633void
634APFloat::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
648APFloat::cmpResult
649APFloat::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.  */
675APFloat::opStatus
676APFloat::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.  */
699bool
700APFloat::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
737APFloat::opStatus
738APFloat::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
846APFloat::opStatus
847APFloat::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.  */
902lostFraction
903APFloat::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
972APFloat::opStatus
973APFloat::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
1016APFloat::opStatus
1017APFloat::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.  */
1061void
1062APFloat::changeSign()
1063{
1064  /* Look mummy, this one's easy.  */
1065  sign = !sign;
1066}
1067
1068void
1069APFloat::clearSign()
1070{
1071  /* So is this one. */
1072  sign = 0;
1073}
1074
1075void
1076APFloat::copySign(const APFloat &rhs)
1077{
1078  /* And this one. */
1079  sign = rhs.sign;
1080}
1081
1082/* Normalized addition or subtraction.  */
1083APFloat::opStatus
1084APFloat::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.  */
1114APFloat::opStatus
1115APFloat::add(const APFloat &rhs, roundingMode rounding_mode)
1116{
1117  return addOrSubtract(rhs, rounding_mode, false);
1118}
1119
1120/* Normalized subtraction.  */
1121APFloat::opStatus
1122APFloat::subtract(const APFloat &rhs, roundingMode rounding_mode)
1123{
1124  return addOrSubtract(rhs, rounding_mode, true);
1125}
1126
1127/* Normalized multiply.  */
1128APFloat::opStatus
1129APFloat::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.  */
1147APFloat::opStatus
1148APFloat::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. */
1166APFloat::opStatus
1167APFloat::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.  */
1190APFloat::opStatus
1191APFloat::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.  */
1235APFloat::cmpResult
1236APFloat::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
1307APFloat::opStatus
1308APFloat::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.  */
1350APFloat::opStatus
1351APFloat::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
1413APFloat::opStatus
1414APFloat::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
1441APFloat::opStatus
1442APFloat::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
1467APFloat::opStatus
1468APFloat::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
1547APFloat::opStatus
1548APFloat::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.
1564uint32_t
1565APFloat::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
1584double
1585APFloat::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
1611float
1612APFloat::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
1637APFloat::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
1663APFloat::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}
1688