APFloat.cpp revision 902ff94aff34acb2b11f4651f4a4006cb788301a 
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  unsigned int origSign = sign;
1172  fs = V.divide(rhs, rmNearestTiesToEven);
1173  if (fs == opDivByZero)
1174    return fs;
1175
1176  int parts = partCount();
1177  integerPart *x = new integerPart[parts];
1178  fs = V.convertToInteger(x, parts * integerPartWidth, true,
1179                          rmNearestTiesToEven);
1180  if (fs==opInvalidOp)
1181    return fs;
1182
1183  fs = V.convertFromInteger(x, parts * integerPartWidth, true,
1184                            rmNearestTiesToEven);
1185  assert(fs==opOK);   // should always work
1186
1187  fs = V.multiply(rhs, rounding_mode);
1188  assert(fs==opOK || fs==opInexact);   // should not overflow or underflow
1189
1190  fs = subtract(V, rounding_mode);
1191  assert(fs==opOK || fs==opInexact);   // likewise
1192
1193  if (isZero())
1194    sign = origSign;    // IEEE754 requires this
1195  delete[] x;
1196  return fs;
1197}
1198
1199/* Normalized fused-multiply-add.  */
1200APFloat::opStatus
1201APFloat::fusedMultiplyAdd(const APFloat &multiplicand,
1202			  const APFloat &addend,
1203			  roundingMode rounding_mode)
1204{
1205  opStatus fs;
1206
1207  /* Post-multiplication sign, before addition.  */
1208  sign ^= multiplicand.sign;
1209
1210  /* If and only if all arguments are normal do we need to do an
1211     extended-precision calculation.  */
1212  if(category == fcNormal
1213     && multiplicand.category == fcNormal
1214     && addend.category == fcNormal) {
1215    lostFraction lost_fraction;
1216
1217    lost_fraction = multiplySignificand(multiplicand, &addend);
1218    fs = normalize(rounding_mode, lost_fraction);
1219    if(lost_fraction != lfExactlyZero)
1220      fs = (opStatus) (fs | opInexact);
1221
1222    /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1223       positive zero unless rounding to minus infinity, except that
1224       adding two like-signed zeroes gives that zero.  */
1225    if(category == fcZero && sign != addend.sign)
1226      sign = (rounding_mode == rmTowardNegative);
1227  } else {
1228    fs = multiplySpecials(multiplicand);
1229
1230    /* FS can only be opOK or opInvalidOp.  There is no more work
1231       to do in the latter case.  The IEEE-754R standard says it is
1232       implementation-defined in this case whether, if ADDEND is a
1233       quiet NaN, we raise invalid op; this implementation does so.
1234
1235       If we need to do the addition we can do so with normal
1236       precision.  */
1237    if(fs == opOK)
1238      fs = addOrSubtract(addend, rounding_mode, false);
1239  }
1240
1241  return fs;
1242}
1243
1244/* Comparison requires normalized numbers.  */
1245APFloat::cmpResult
1246APFloat::compare(const APFloat &rhs) const
1247{
1248  cmpResult result;
1249
1250  assert(semantics == rhs.semantics);
1251
1252  switch(convolve(category, rhs.category)) {
1253  default:
1254    assert(0);
1255
1256  case convolve(fcNaN, fcZero):
1257  case convolve(fcNaN, fcNormal):
1258  case convolve(fcNaN, fcInfinity):
1259  case convolve(fcNaN, fcNaN):
1260  case convolve(fcZero, fcNaN):
1261  case convolve(fcNormal, fcNaN):
1262  case convolve(fcInfinity, fcNaN):
1263    return cmpUnordered;
1264
1265  case convolve(fcInfinity, fcNormal):
1266  case convolve(fcInfinity, fcZero):
1267  case convolve(fcNormal, fcZero):
1268    if(sign)
1269      return cmpLessThan;
1270    else
1271      return cmpGreaterThan;
1272
1273  case convolve(fcNormal, fcInfinity):
1274  case convolve(fcZero, fcInfinity):
1275  case convolve(fcZero, fcNormal):
1276    if(rhs.sign)
1277      return cmpGreaterThan;
1278    else
1279      return cmpLessThan;
1280
1281  case convolve(fcInfinity, fcInfinity):
1282    if(sign == rhs.sign)
1283      return cmpEqual;
1284    else if(sign)
1285      return cmpLessThan;
1286    else
1287      return cmpGreaterThan;
1288
1289  case convolve(fcZero, fcZero):
1290    return cmpEqual;
1291
1292  case convolve(fcNormal, fcNormal):
1293    break;
1294  }
1295
1296  /* Two normal numbers.  Do they have the same sign?  */
1297  if(sign != rhs.sign) {
1298    if(sign)
1299      result = cmpLessThan;
1300    else
1301      result = cmpGreaterThan;
1302  } else {
1303    /* Compare absolute values; invert result if negative.  */
1304    result = compareAbsoluteValue(rhs);
1305
1306    if(sign) {
1307      if(result == cmpLessThan)
1308	result = cmpGreaterThan;
1309      else if(result == cmpGreaterThan)
1310	result = cmpLessThan;
1311    }
1312  }
1313
1314  return result;
1315}
1316
1317APFloat::opStatus
1318APFloat::convert(const fltSemantics &toSemantics,
1319		 roundingMode rounding_mode)
1320{
1321  lostFraction lostFraction;
1322  unsigned int newPartCount, oldPartCount;
1323  opStatus fs;
1324
1325  lostFraction = lfExactlyZero;
1326  newPartCount = partCountForBits(toSemantics.precision + 1);
1327  oldPartCount = partCount();
1328
1329  /* Handle storage complications.  If our new form is wider,
1330     re-allocate our bit pattern into wider storage.  If it is
1331     narrower, we ignore the excess parts, but if narrowing to a
1332     single part we need to free the old storage.
1333     Be careful not to reference significandParts for zeroes
1334     and infinities, since it aborts.  */
1335  if (newPartCount > oldPartCount) {
1336    integerPart *newParts;
1337    newParts = new integerPart[newPartCount];
1338    APInt::tcSet(newParts, 0, newPartCount);
1339    if (category==fcNormal || category==fcNaN)
1340      APInt::tcAssign(newParts, significandParts(), oldPartCount);
1341    freeSignificand();
1342    significand.parts = newParts;
1343  } else if (newPartCount < oldPartCount) {
1344    /* Capture any lost fraction through truncation of parts so we get
1345       correct rounding whilst normalizing.  */
1346    if (category==fcNormal)
1347      lostFraction = lostFractionThroughTruncation
1348        (significandParts(), oldPartCount, toSemantics.precision);
1349    if (newPartCount == 1) {
1350        integerPart newPart = 0;
1351        if (category==fcNormal || category==fcNaN)
1352          newPart = significandParts()[0];
1353        freeSignificand();
1354        significand.part = newPart;
1355    }
1356  }
1357
1358  if(category == fcNormal) {
1359    /* Re-interpret our bit-pattern.  */
1360    exponent += toSemantics.precision - semantics->precision;
1361    semantics = &toSemantics;
1362    fs = normalize(rounding_mode, lostFraction);
1363  } else if (category == fcNaN) {
1364    int shift = toSemantics.precision - semantics->precision;
1365    // No normalization here, just truncate
1366    if (shift>0)
1367      APInt::tcShiftLeft(significandParts(), newPartCount, shift);
1368    else if (shift < 0)
1369      APInt::tcShiftRight(significandParts(), newPartCount, -shift);
1370    // gcc forces the Quiet bit on, which means (float)(double)(float_sNan)
1371    // does not give you back the same bits.  This is dubious, and we
1372    // don't currently do it.  You're really supposed to get
1373    // an invalid operation signal at runtime, but nobody does that.
1374    semantics = &toSemantics;
1375    fs = opOK;
1376  } else {
1377    semantics = &toSemantics;
1378    fs = opOK;
1379  }
1380
1381  return fs;
1382}
1383
1384/* Convert a floating point number to an integer according to the
1385   rounding mode.  If the rounded integer value is out of range this
1386   returns an invalid operation exception.  If the rounded value is in
1387   range but the floating point number is not the exact integer, the C
1388   standard doesn't require an inexact exception to be raised.  IEEE
1389   854 does require it so we do that.
1390
1391   Note that for conversions to integer type the C standard requires
1392   round-to-zero to always be used.  */
1393APFloat::opStatus
1394APFloat::convertToInteger(integerPart *parts, unsigned int width,
1395			  bool isSigned,
1396			  roundingMode rounding_mode) const
1397{
1398  lostFraction lost_fraction;
1399  unsigned int msb, partsCount;
1400  int bits;
1401
1402  /* Handle the three special cases first.  */
1403  if(category == fcInfinity || category == fcNaN)
1404    return opInvalidOp;
1405
1406  partsCount = partCountForBits(width);
1407
1408  if(category == fcZero) {
1409    APInt::tcSet(parts, 0, partsCount);
1410    return opOK;
1411  }
1412
1413  /* Shift the bit pattern so the fraction is lost.  */
1414  APFloat tmp(*this);
1415
1416  bits = (int) semantics->precision - 1 - exponent;
1417
1418  if(bits > 0) {
1419    lost_fraction = tmp.shiftSignificandRight(bits);
1420  } else {
1421    tmp.shiftSignificandLeft(-bits);
1422    lost_fraction = lfExactlyZero;
1423  }
1424
1425  if(lost_fraction != lfExactlyZero
1426     && tmp.roundAwayFromZero(rounding_mode, lost_fraction))
1427    tmp.incrementSignificand();
1428
1429  msb = tmp.significandMSB();
1430
1431  /* Negative numbers cannot be represented as unsigned.  */
1432  if(!isSigned && tmp.sign && msb != -1U)
1433    return opInvalidOp;
1434
1435  /* It takes exponent + 1 bits to represent the truncated floating
1436     point number without its sign.  We lose a bit for the sign, but
1437     the maximally negative integer is a special case.  */
1438  if(msb + 1 > width)		/* !! Not same as msb >= width !! */
1439    return opInvalidOp;
1440
1441  if(isSigned && msb + 1 == width
1442     && (!tmp.sign || tmp.significandLSB() != msb))
1443    return opInvalidOp;
1444
1445  APInt::tcAssign(parts, tmp.significandParts(), partsCount);
1446
1447  if(tmp.sign)
1448    APInt::tcNegate(parts, partsCount);
1449
1450  if(lost_fraction == lfExactlyZero)
1451    return opOK;
1452  else
1453    return opInexact;
1454}
1455
1456APFloat::opStatus
1457APFloat::convertFromUnsignedInteger(integerPart *parts,
1458				    unsigned int partCount,
1459				    roundingMode rounding_mode)
1460{
1461  unsigned int msb, precision;
1462  lostFraction lost_fraction;
1463
1464  msb = APInt::tcMSB(parts, partCount) + 1;
1465  precision = semantics->precision;
1466
1467  category = fcNormal;
1468  exponent = precision - 1;
1469
1470  if(msb > precision) {
1471    exponent += (msb - precision);
1472    lost_fraction = shiftRight(parts, partCount, msb - precision);
1473    msb = precision;
1474  } else
1475    lost_fraction = lfExactlyZero;
1476
1477  /* Copy the bit image.  */
1478  zeroSignificand();
1479  APInt::tcAssign(significandParts(), parts, partCountForBits(msb));
1480
1481  return normalize(rounding_mode, lost_fraction);
1482}
1483
1484APFloat::opStatus
1485APFloat::convertFromInteger(const integerPart *parts, unsigned int width,
1486			    bool isSigned, roundingMode rounding_mode)
1487{
1488  unsigned int partCount = partCountForBits(width);
1489  opStatus status;
1490  APInt api = APInt(width, partCount, parts);
1491  integerPart *copy = new integerPart[partCount];
1492
1493  sign = false;
1494  if(isSigned) {
1495    if (APInt::tcExtractBit(parts, width - 1)) {
1496      sign = true;
1497      if (width < partCount * integerPartWidth)
1498        api = api.sext(partCount * integerPartWidth);
1499    }
1500    else if (width < partCount * integerPartWidth)
1501      api = api.zext(partCount * integerPartWidth);
1502  } else {
1503    if (width < partCount * integerPartWidth)
1504      api = api.zext(partCount * integerPartWidth);
1505  }
1506
1507  APInt::tcAssign(copy, api.getRawData(), partCount);
1508  status = convertFromUnsignedInteger(copy, partCount, rounding_mode);
1509  return status;
1510}
1511
1512APFloat::opStatus
1513APFloat::convertFromHexadecimalString(const char *p,
1514				      roundingMode rounding_mode)
1515{
1516  lostFraction lost_fraction;
1517  integerPart *significand;
1518  unsigned int bitPos, partsCount;
1519  const char *dot, *firstSignificantDigit;
1520
1521  zeroSignificand();
1522  exponent = 0;
1523  category = fcNormal;
1524
1525  significand = significandParts();
1526  partsCount = partCount();
1527  bitPos = partsCount * integerPartWidth;
1528
1529  /* Skip leading zeroes and any(hexa)decimal point.  */
1530  p = skipLeadingZeroesAndAnyDot(p, &dot);
1531  firstSignificantDigit = p;
1532
1533  for(;;) {
1534    integerPart hex_value;
1535
1536    if(*p == '.') {
1537      assert(dot == 0);
1538      dot = p++;
1539    }
1540
1541    hex_value = hexDigitValue(*p);
1542    if(hex_value == -1U) {
1543      lost_fraction = lfExactlyZero;
1544      break;
1545    }
1546
1547    p++;
1548
1549    /* Store the number whilst 4-bit nibbles remain.  */
1550    if(bitPos) {
1551      bitPos -= 4;
1552      hex_value <<= bitPos % integerPartWidth;
1553      significand[bitPos / integerPartWidth] |= hex_value;
1554    } else {
1555      lost_fraction = trailingHexadecimalFraction(p, hex_value);
1556      while(hexDigitValue(*p) != -1U)
1557	p++;
1558      break;
1559    }
1560  }
1561
1562  /* Hex floats require an exponent but not a hexadecimal point.  */
1563  assert(*p == 'p' || *p == 'P');
1564
1565  /* Ignore the exponent if we are zero.  */
1566  if(p != firstSignificantDigit) {
1567    int expAdjustment;
1568
1569    /* Implicit hexadecimal point?  */
1570    if(!dot)
1571      dot = p;
1572
1573    /* Calculate the exponent adjustment implicit in the number of
1574       significant digits.  */
1575    expAdjustment = dot - firstSignificantDigit;
1576    if(expAdjustment < 0)
1577      expAdjustment++;
1578    expAdjustment = expAdjustment * 4 - 1;
1579
1580    /* Adjust for writing the significand starting at the most
1581       significant nibble.  */
1582    expAdjustment += semantics->precision;
1583    expAdjustment -= partsCount * integerPartWidth;
1584
1585    /* Adjust for the given exponent.  */
1586    exponent = totalExponent(p, expAdjustment);
1587  }
1588
1589  return normalize(rounding_mode, lost_fraction);
1590}
1591
1592APFloat::opStatus
1593APFloat::convertFromString(const char *p, roundingMode rounding_mode) {
1594  /* Handle a leading minus sign.  */
1595  if(*p == '-')
1596    sign = 1, p++;
1597  else
1598    sign = 0;
1599
1600  if(p[0] == '0' && (p[1] == 'x' || p[1] == 'X'))
1601    return convertFromHexadecimalString(p + 2, rounding_mode);
1602
1603  assert(0 && "Decimal to binary conversions not yet implemented");
1604  abort();
1605}
1606
1607// For good performance it is desirable for different APFloats
1608// to produce different integers.
1609uint32_t
1610APFloat::getHashValue() const {
1611  if (category==fcZero) return sign<<8 | semantics->precision ;
1612  else if (category==fcInfinity) return sign<<9 | semantics->precision;
1613  else if (category==fcNaN) return 1<<10 | semantics->precision;
1614  else {
1615    uint32_t hash = sign<<11 | semantics->precision | exponent<<12;
1616    const integerPart* p = significandParts();
1617    for (int i=partCount(); i>0; i--, p++)
1618      hash ^= ((uint32_t)*p) ^ (*p)>>32;
1619    return hash;
1620  }
1621}
1622
1623// Conversion from APFloat to/from host float/double.  It may eventually be
1624// possible to eliminate these and have everybody deal with APFloats, but that
1625// will take a while.  This approach will not easily extend to long double.
1626// Current implementation requires integerPartWidth==64, which is correct at
1627// the moment but could be made more general.
1628
1629// Denormals have exponent minExponent in APFloat, but minExponent-1 in
1630// the actual IEEE respresentations.  We compensate for that here.
1631
1632APInt
1633APFloat::convertF80LongDoubleAPFloatToAPInt() const {
1634  assert(semantics == (const llvm::fltSemantics* const)&x87DoubleExtended);
1635  assert (partCount()==2);
1636
1637  uint64_t myexponent, mysignificand;
1638
1639  if (category==fcNormal) {
1640    myexponent = exponent+16383; //bias
1641    mysignificand = significandParts()[0];
1642    if (myexponent==1 && !(mysignificand & 0x8000000000000000ULL))
1643      myexponent = 0;   // denormal
1644  } else if (category==fcZero) {
1645    myexponent = 0;
1646    mysignificand = 0;
1647  } else if (category==fcInfinity) {
1648    myexponent = 0x7fff;
1649    mysignificand = 0x8000000000000000ULL;
1650  } else if (category==fcNaN) {
1651    myexponent = 0x7fff;
1652    mysignificand = significandParts()[0];
1653  } else
1654    assert(0);
1655
1656  uint64_t words[2];
1657  words[0] =  (((uint64_t)sign & 1) << 63) |
1658              ((myexponent & 0x7fff) <<  48) |
1659              ((mysignificand >>16) & 0xffffffffffffLL);
1660  words[1] = mysignificand & 0xffff;
1661  APInt api(80, 2, words);
1662  return api;
1663}
1664
1665APInt
1666APFloat::convertDoubleAPFloatToAPInt() const {
1667  assert(semantics == (const llvm::fltSemantics*)&IEEEdouble);
1668  assert (partCount()==1);
1669
1670  uint64_t myexponent, mysignificand;
1671
1672  if (category==fcNormal) {
1673    myexponent = exponent+1023; //bias
1674    mysignificand = *significandParts();
1675    if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
1676      myexponent = 0;   // denormal
1677  } else if (category==fcZero) {
1678    myexponent = 0;
1679    mysignificand = 0;
1680  } else if (category==fcInfinity) {
1681    myexponent = 0x7ff;
1682    mysignificand = 0;
1683  } else if (category==fcNaN) {
1684    myexponent = 0x7ff;
1685    mysignificand = *significandParts();
1686  } else
1687    assert(0);
1688
1689  APInt api(64, (((((uint64_t)sign & 1) << 63) |
1690                 ((myexponent & 0x7ff) <<  52) |
1691                 (mysignificand & 0xfffffffffffffLL))));
1692  return api;
1693}
1694
1695APInt
1696APFloat::convertFloatAPFloatToAPInt() const {
1697  assert(semantics == (const llvm::fltSemantics*)&IEEEsingle);
1698  assert (partCount()==1);
1699
1700  uint32_t myexponent, mysignificand;
1701
1702  if (category==fcNormal) {
1703    myexponent = exponent+127; //bias
1704    mysignificand = *significandParts();
1705    if (myexponent == 1 && !(mysignificand & 0x400000))
1706      myexponent = 0;   // denormal
1707  } else if (category==fcZero) {
1708    myexponent = 0;
1709    mysignificand = 0;
1710  } else if (category==fcInfinity) {
1711    myexponent = 0xff;
1712    mysignificand = 0;
1713  } else if (category==fcNaN) {
1714    myexponent = 0xff;
1715    mysignificand = *significandParts();
1716  } else
1717    assert(0);
1718
1719  APInt api(32, (((sign&1) << 31) | ((myexponent&0xff) << 23) |
1720                 (mysignificand & 0x7fffff)));
1721  return api;
1722}
1723
1724APInt
1725APFloat::convertToAPInt() const {
1726  if (semantics == (const llvm::fltSemantics* const)&IEEEsingle)
1727    return convertFloatAPFloatToAPInt();
1728  else if (semantics == (const llvm::fltSemantics* const)&IEEEdouble)
1729    return convertDoubleAPFloatToAPInt();
1730  else if (semantics == (const llvm::fltSemantics* const)&x87DoubleExtended)
1731    return convertF80LongDoubleAPFloatToAPInt();
1732  else
1733    assert(0);
1734}
1735
1736float
1737APFloat::convertToFloat() const {
1738  assert(semantics == (const llvm::fltSemantics* const)&IEEEsingle);
1739  APInt api = convertToAPInt();
1740  return api.bitsToFloat();
1741}
1742
1743double
1744APFloat::convertToDouble() const {
1745  assert(semantics == (const llvm::fltSemantics* const)&IEEEdouble);
1746  APInt api = convertToAPInt();
1747  return api.bitsToDouble();
1748}
1749
1750/// Integer bit is explicit in this format.  Current Intel book does not
1751/// define meaning of:
1752///  exponent = all 1's, integer bit not set.
1753///  exponent = 0, integer bit set. (formerly "psuedodenormals")
1754///  exponent!=0 nor all 1's, integer bit not set. (formerly "unnormals")
1755void
1756APFloat::initFromF80LongDoubleAPInt(const APInt &api) {
1757  assert(api.getBitWidth()==80);
1758  uint64_t i1 = api.getRawData()[0];
1759  uint64_t i2 = api.getRawData()[1];
1760  uint64_t myexponent = (i1 >> 48) & 0x7fff;
1761  uint64_t mysignificand = ((i1 << 16) &  0xffffffffffff0000ULL) |
1762                          (i2 & 0xffff);
1763
1764  initialize(&APFloat::x87DoubleExtended);
1765  assert(partCount()==2);
1766
1767  sign = i1>>63;
1768  if (myexponent==0 && mysignificand==0) {
1769    // exponent, significand meaningless
1770    category = fcZero;
1771  } else if (myexponent==0x7fff && mysignificand==0x8000000000000000ULL) {
1772    // exponent, significand meaningless
1773    category = fcInfinity;
1774  } else if (myexponent==0x7fff && mysignificand!=0x8000000000000000ULL) {
1775    // exponent meaningless
1776    category = fcNaN;
1777    significandParts()[0] = mysignificand;
1778    significandParts()[1] = 0;
1779  } else {
1780    category = fcNormal;
1781    exponent = myexponent - 16383;
1782    significandParts()[0] = mysignificand;
1783    significandParts()[1] = 0;
1784    if (myexponent==0)          // denormal
1785      exponent = -16382;
1786 }
1787}
1788
1789void
1790APFloat::initFromDoubleAPInt(const APInt &api) {
1791  assert(api.getBitWidth()==64);
1792  uint64_t i = *api.getRawData();
1793  uint64_t myexponent = (i >> 52) & 0x7ff;
1794  uint64_t mysignificand = i & 0xfffffffffffffLL;
1795
1796  initialize(&APFloat::IEEEdouble);
1797  assert(partCount()==1);
1798
1799  sign = i>>63;
1800  if (myexponent==0 && mysignificand==0) {
1801    // exponent, significand meaningless
1802    category = fcZero;
1803  } else if (myexponent==0x7ff && mysignificand==0) {
1804    // exponent, significand meaningless
1805    category = fcInfinity;
1806  } else if (myexponent==0x7ff && mysignificand!=0) {
1807    // exponent meaningless
1808    category = fcNaN;
1809    *significandParts() = mysignificand;
1810  } else {
1811    category = fcNormal;
1812    exponent = myexponent - 1023;
1813    *significandParts() = mysignificand;
1814    if (myexponent==0)          // denormal
1815      exponent = -1022;
1816    else
1817      *significandParts() |= 0x10000000000000LL;  // integer bit
1818 }
1819}
1820
1821void
1822APFloat::initFromFloatAPInt(const APInt & api) {
1823  assert(api.getBitWidth()==32);
1824  uint32_t i = (uint32_t)*api.getRawData();
1825  uint32_t myexponent = (i >> 23) & 0xff;
1826  uint32_t mysignificand = i & 0x7fffff;
1827
1828  initialize(&APFloat::IEEEsingle);
1829  assert(partCount()==1);
1830
1831  sign = i >> 31;
1832  if (myexponent==0 && mysignificand==0) {
1833    // exponent, significand meaningless
1834    category = fcZero;
1835  } else if (myexponent==0xff && mysignificand==0) {
1836    // exponent, significand meaningless
1837    category = fcInfinity;
1838  } else if (myexponent==0xff && mysignificand!=0) {
1839    // sign, exponent, significand meaningless
1840    category = fcNaN;
1841    *significandParts() = mysignificand;
1842  } else {
1843    category = fcNormal;
1844    exponent = myexponent - 127;  //bias
1845    *significandParts() = mysignificand;
1846    if (myexponent==0)    // denormal
1847      exponent = -126;
1848    else
1849      *significandParts() |= 0x800000; // integer bit
1850  }
1851}
1852
1853/// Treat api as containing the bits of a floating point number.  Currently
1854/// we infer the floating point type from the size of the APInt.  FIXME: This
1855/// breaks when we get to PPC128 and IEEE128 (but both cannot exist in the
1856/// same compile...)
1857void
1858APFloat::initFromAPInt(const APInt& api) {
1859  if (api.getBitWidth() == 32)
1860    return initFromFloatAPInt(api);
1861  else if (api.getBitWidth()==64)
1862    return initFromDoubleAPInt(api);
1863  else if (api.getBitWidth()==80)
1864    return initFromF80LongDoubleAPInt(api);
1865  else
1866    assert(0);
1867}
1868
1869APFloat::APFloat(const APInt& api) {
1870  initFromAPInt(api);
1871}
1872
1873APFloat::APFloat(float f) {
1874  APInt api = APInt(32, 0);
1875  initFromAPInt(api.floatToBits(f));
1876}
1877
1878APFloat::APFloat(double d) {
1879  APInt api = APInt(64, 0);
1880  initFromAPInt(api.doubleToBits(d));
1881}
1882
1883