APFloat.cpp revision 3f6eb7419de437436265831fce92f62498556e08 
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 +
343                          semantics->implicitIntegerBit ? 1 : 0);
344}
345
346unsigned int
347APFloat::semanticsPrecision(const fltSemantics &semantics)
348{
349  return semantics.precision;
350}
351
352const integerPart *
353APFloat::significandParts() const
354{
355  return const_cast<APFloat *>(this)->significandParts();
356}
357
358integerPart *
359APFloat::significandParts()
360{
361  assert(category == fcNormal || category == fcNaN);
362
363  if(partCount() > 1)
364    return significand.parts;
365  else
366    return &significand.part;
367}
368
369/* Combine the effect of two lost fractions.  */
370lostFraction
371APFloat::combineLostFractions(lostFraction moreSignificant,
372			      lostFraction lessSignificant)
373{
374  if(lessSignificant != lfExactlyZero) {
375    if(moreSignificant == lfExactlyZero)
376      moreSignificant = lfLessThanHalf;
377    else if(moreSignificant == lfExactlyHalf)
378      moreSignificant = lfMoreThanHalf;
379  }
380
381  return moreSignificant;
382}
383
384void
385APFloat::zeroSignificand()
386{
387  category = fcNormal;
388  APInt::tcSet(significandParts(), 0, partCount());
389}
390
391/* Increment an fcNormal floating point number's significand.  */
392void
393APFloat::incrementSignificand()
394{
395  integerPart carry;
396
397  carry = APInt::tcIncrement(significandParts(), partCount());
398
399  /* Our callers should never cause us to overflow.  */
400  assert(carry == 0);
401}
402
403/* Add the significand of the RHS.  Returns the carry flag.  */
404integerPart
405APFloat::addSignificand(const APFloat &rhs)
406{
407  integerPart *parts;
408
409  parts = significandParts();
410
411  assert(semantics == rhs.semantics);
412  assert(exponent == rhs.exponent);
413
414  return APInt::tcAdd(parts, rhs.significandParts(), 0, partCount());
415}
416
417/* Subtract the significand of the RHS with a borrow flag.  Returns
418   the borrow flag.  */
419integerPart
420APFloat::subtractSignificand(const APFloat &rhs, integerPart borrow)
421{
422  integerPart *parts;
423
424  parts = significandParts();
425
426  assert(semantics == rhs.semantics);
427  assert(exponent == rhs.exponent);
428
429  return APInt::tcSubtract(parts, rhs.significandParts(), borrow,
430			   partCount());
431}
432
433/* Multiply the significand of the RHS.  If ADDEND is non-NULL, add it
434   on to the full-precision result of the multiplication.  Returns the
435   lost fraction.  */
436lostFraction
437APFloat::multiplySignificand(const APFloat &rhs, const APFloat *addend)
438{
439  unsigned int omsb;	// One, not zero, based MSB.
440  unsigned int partsCount, newPartsCount, precision;
441  integerPart *lhsSignificand;
442  integerPart scratch[4];
443  integerPart *fullSignificand;
444  lostFraction lost_fraction;
445
446  assert(semantics == rhs.semantics);
447
448  precision = semantics->precision;
449  newPartsCount = partCountForBits(precision * 2);
450
451  if(newPartsCount > 4)
452    fullSignificand = new integerPart[newPartsCount];
453  else
454    fullSignificand = scratch;
455
456  lhsSignificand = significandParts();
457  partsCount = partCount();
458
459  APInt::tcFullMultiply(fullSignificand, lhsSignificand,
460			rhs.significandParts(), partsCount);
461
462  lost_fraction = lfExactlyZero;
463  omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
464  exponent += rhs.exponent;
465
466  if(addend) {
467    Significand savedSignificand = significand;
468    const fltSemantics *savedSemantics = semantics;
469    fltSemantics extendedSemantics;
470    opStatus status;
471    unsigned int extendedPrecision;
472
473    /* Normalize our MSB.  */
474    extendedPrecision = precision + precision - 1;
475    if(omsb != extendedPrecision)
476      {
477	APInt::tcShiftLeft(fullSignificand, newPartsCount,
478			   extendedPrecision - omsb);
479	exponent -= extendedPrecision - omsb;
480      }
481
482    /* Create new semantics.  */
483    extendedSemantics = *semantics;
484    extendedSemantics.precision = extendedPrecision;
485
486    if(newPartsCount == 1)
487      significand.part = fullSignificand[0];
488    else
489      significand.parts = fullSignificand;
490    semantics = &extendedSemantics;
491
492    APFloat extendedAddend(*addend);
493    status = extendedAddend.convert(extendedSemantics, rmTowardZero);
494    assert(status == opOK);
495    lost_fraction = addOrSubtractSignificand(extendedAddend, false);
496
497    /* Restore our state.  */
498    if(newPartsCount == 1)
499      fullSignificand[0] = significand.part;
500    significand = savedSignificand;
501    semantics = savedSemantics;
502
503    omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
504  }
505
506  exponent -= (precision - 1);
507
508  if(omsb > precision) {
509    unsigned int bits, significantParts;
510    lostFraction lf;
511
512    bits = omsb - precision;
513    significantParts = partCountForBits(omsb);
514    lf = shiftRight(fullSignificand, significantParts, bits);
515    lost_fraction = combineLostFractions(lf, lost_fraction);
516    exponent += bits;
517  }
518
519  APInt::tcAssign(lhsSignificand, fullSignificand, partsCount);
520
521  if(newPartsCount > 4)
522    delete [] fullSignificand;
523
524  return lost_fraction;
525}
526
527/* Multiply the significands of LHS and RHS to DST.  */
528lostFraction
529APFloat::divideSignificand(const APFloat &rhs)
530{
531  unsigned int bit, i, partsCount;
532  const integerPart *rhsSignificand;
533  integerPart *lhsSignificand, *dividend, *divisor;
534  integerPart scratch[4];
535  lostFraction lost_fraction;
536
537  assert(semantics == rhs.semantics);
538
539  lhsSignificand = significandParts();
540  rhsSignificand = rhs.significandParts();
541  partsCount = partCount();
542
543  if(partsCount > 2)
544    dividend = new integerPart[partsCount * 2];
545  else
546    dividend = scratch;
547
548  divisor = dividend + partsCount;
549
550  /* Copy the dividend and divisor as they will be modified in-place.  */
551  for(i = 0; i < partsCount; i++) {
552    dividend[i] = lhsSignificand[i];
553    divisor[i] = rhsSignificand[i];
554    lhsSignificand[i] = 0;
555  }
556
557  exponent -= rhs.exponent;
558
559  unsigned int precision = semantics->precision;
560
561  /* Normalize the divisor.  */
562  bit = precision - APInt::tcMSB(divisor, partsCount) - 1;
563  if(bit) {
564    exponent += bit;
565    APInt::tcShiftLeft(divisor, partsCount, bit);
566  }
567
568  /* Normalize the dividend.  */
569  bit = precision - APInt::tcMSB(dividend, partsCount) - 1;
570  if(bit) {
571    exponent -= bit;
572    APInt::tcShiftLeft(dividend, partsCount, bit);
573  }
574
575  if(APInt::tcCompare(dividend, divisor, partsCount) < 0) {
576    exponent--;
577    APInt::tcShiftLeft(dividend, partsCount, 1);
578    assert(APInt::tcCompare(dividend, divisor, partsCount) >= 0);
579  }
580
581  /* Long division.  */
582  for(bit = precision; bit; bit -= 1) {
583    if(APInt::tcCompare(dividend, divisor, partsCount) >= 0) {
584      APInt::tcSubtract(dividend, divisor, 0, partsCount);
585      APInt::tcSetBit(lhsSignificand, bit - 1);
586    }
587
588    APInt::tcShiftLeft(dividend, partsCount, 1);
589  }
590
591  /* Figure out the lost fraction.  */
592  int cmp = APInt::tcCompare(dividend, divisor, partsCount);
593
594  if(cmp > 0)
595    lost_fraction = lfMoreThanHalf;
596  else if(cmp == 0)
597    lost_fraction = lfExactlyHalf;
598  else if(APInt::tcIsZero(dividend, partsCount))
599    lost_fraction = lfExactlyZero;
600  else
601    lost_fraction = lfLessThanHalf;
602
603  if(partsCount > 2)
604    delete [] dividend;
605
606  return lost_fraction;
607}
608
609unsigned int
610APFloat::significandMSB() const
611{
612  return APInt::tcMSB(significandParts(), partCount());
613}
614
615unsigned int
616APFloat::significandLSB() const
617{
618  return APInt::tcLSB(significandParts(), partCount());
619}
620
621/* Note that a zero result is NOT normalized to fcZero.  */
622lostFraction
623APFloat::shiftSignificandRight(unsigned int bits)
624{
625  /* Our exponent should not overflow.  */
626  assert((exponent_t) (exponent + bits) >= exponent);
627
628  exponent += bits;
629
630  return shiftRight(significandParts(), partCount(), bits);
631}
632
633/* Shift the significand left BITS bits, subtract BITS from its exponent.  */
634void
635APFloat::shiftSignificandLeft(unsigned int bits)
636{
637  assert(bits < semantics->precision);
638
639  if(bits) {
640    unsigned int partsCount = partCount();
641
642    APInt::tcShiftLeft(significandParts(), partsCount, bits);
643    exponent -= bits;
644
645    assert(!APInt::tcIsZero(significandParts(), partsCount));
646  }
647}
648
649APFloat::cmpResult
650APFloat::compareAbsoluteValue(const APFloat &rhs) const
651{
652  int compare;
653
654  assert(semantics == rhs.semantics);
655  assert(category == fcNormal);
656  assert(rhs.category == fcNormal);
657
658  compare = exponent - rhs.exponent;
659
660  /* If exponents are equal, do an unsigned bignum comparison of the
661     significands.  */
662  if(compare == 0)
663    compare = APInt::tcCompare(significandParts(), rhs.significandParts(),
664			       partCount());
665
666  if(compare > 0)
667    return cmpGreaterThan;
668  else if(compare < 0)
669    return cmpLessThan;
670  else
671    return cmpEqual;
672}
673
674/* Handle overflow.  Sign is preserved.  We either become infinity or
675   the largest finite number.  */
676APFloat::opStatus
677APFloat::handleOverflow(roundingMode rounding_mode)
678{
679  /* Infinity?  */
680  if(rounding_mode == rmNearestTiesToEven
681     || rounding_mode == rmNearestTiesToAway
682     || (rounding_mode == rmTowardPositive && !sign)
683     || (rounding_mode == rmTowardNegative && sign))
684    {
685      category = fcInfinity;
686      return (opStatus) (opOverflow | opInexact);
687    }
688
689  /* Otherwise we become the largest finite number.  */
690  category = fcNormal;
691  exponent = semantics->maxExponent;
692  APInt::tcSetLeastSignificantBits(significandParts(), partCount(),
693				   semantics->precision);
694
695  return opInexact;
696}
697
698/* This routine must work for fcZero of both signs, and fcNormal
699   numbers.  */
700bool
701APFloat::roundAwayFromZero(roundingMode rounding_mode,
702			   lostFraction lost_fraction)
703{
704  /* NaNs and infinities should not have lost fractions.  */
705  assert(category == fcNormal || category == fcZero);
706
707  /* Our caller has already handled this case.  */
708  assert(lost_fraction != lfExactlyZero);
709
710  switch(rounding_mode) {
711  default:
712    assert(0);
713
714  case rmNearestTiesToAway:
715    return lost_fraction == lfExactlyHalf || lost_fraction == lfMoreThanHalf;
716
717  case rmNearestTiesToEven:
718    if(lost_fraction == lfMoreThanHalf)
719      return true;
720
721    /* Our zeroes don't have a significand to test.  */
722    if(lost_fraction == lfExactlyHalf && category != fcZero)
723      return significandParts()[0] & 1;
724
725    return false;
726
727  case rmTowardZero:
728    return false;
729
730  case rmTowardPositive:
731    return sign == false;
732
733  case rmTowardNegative:
734    return sign == true;
735  }
736}
737
738APFloat::opStatus
739APFloat::normalize(roundingMode rounding_mode,
740		   lostFraction lost_fraction)
741{
742  unsigned int omsb;		/* One, not zero, based MSB.  */
743  int exponentChange;
744
745  if(category != fcNormal)
746    return opOK;
747
748  /* Before rounding normalize the exponent of fcNormal numbers.  */
749  omsb = significandMSB() + 1;
750
751  if(omsb) {
752    /* OMSB is numbered from 1.  We want to place it in the integer
753       bit numbered PRECISON if possible, with a compensating change in
754       the exponent.  */
755    exponentChange = omsb - semantics->precision;
756
757    /* If the resulting exponent is too high, overflow according to
758       the rounding mode.  */
759    if(exponent + exponentChange > semantics->maxExponent)
760      return handleOverflow(rounding_mode);
761
762    /* Subnormal numbers have exponent minExponent, and their MSB
763       is forced based on that.  */
764    if(exponent + exponentChange < semantics->minExponent)
765      exponentChange = semantics->minExponent - exponent;
766
767    /* Shifting left is easy as we don't lose precision.  */
768    if(exponentChange < 0) {
769      assert(lost_fraction == lfExactlyZero);
770
771      shiftSignificandLeft(-exponentChange);
772
773      return opOK;
774    }
775
776    if(exponentChange > 0) {
777      lostFraction lf;
778
779      /* Shift right and capture any new lost fraction.  */
780      lf = shiftSignificandRight(exponentChange);
781
782      lost_fraction = combineLostFractions(lf, lost_fraction);
783
784      /* Keep OMSB up-to-date.  */
785      if(omsb > (unsigned) exponentChange)
786	omsb -= (unsigned) exponentChange;
787      else
788	omsb = 0;
789    }
790  }
791
792  /* Now round the number according to rounding_mode given the lost
793     fraction.  */
794
795  /* As specified in IEEE 754, since we do not trap we do not report
796     underflow for exact results.  */
797  if(lost_fraction == lfExactlyZero) {
798    /* Canonicalize zeroes.  */
799    if(omsb == 0)
800      category = fcZero;
801
802    return opOK;
803  }
804
805  /* Increment the significand if we're rounding away from zero.  */
806  if(roundAwayFromZero(rounding_mode, lost_fraction)) {
807    if(omsb == 0)
808      exponent = semantics->minExponent;
809
810    incrementSignificand();
811    omsb = significandMSB() + 1;
812
813    /* Did the significand increment overflow?  */
814    if(omsb == (unsigned) semantics->precision + 1) {
815      /* Renormalize by incrementing the exponent and shifting our
816	 significand right one.  However if we already have the
817	 maximum exponent we overflow to infinity.  */
818      if(exponent == semantics->maxExponent) {
819	category = fcInfinity;
820
821	return (opStatus) (opOverflow | opInexact);
822      }
823
824      shiftSignificandRight(1);
825
826      return opInexact;
827    }
828  }
829
830  /* The normal case - we were and are not denormal, and any
831     significand increment above didn't overflow.  */
832  if(omsb == semantics->precision)
833    return opInexact;
834
835  /* We have a non-zero denormal.  */
836  assert(omsb < semantics->precision);
837  assert(exponent == semantics->minExponent);
838
839  /* Canonicalize zeroes.  */
840  if(omsb == 0)
841    category = fcZero;
842
843  /* The fcZero case is a denormal that underflowed to zero.  */
844  return (opStatus) (opUnderflow | opInexact);
845}
846
847APFloat::opStatus
848APFloat::addOrSubtractSpecials(const APFloat &rhs, bool subtract)
849{
850  switch(convolve(category, rhs.category)) {
851  default:
852    assert(0);
853
854  case convolve(fcNaN, fcZero):
855  case convolve(fcNaN, fcNormal):
856  case convolve(fcNaN, fcInfinity):
857  case convolve(fcNaN, fcNaN):
858  case convolve(fcNormal, fcZero):
859  case convolve(fcInfinity, fcNormal):
860  case convolve(fcInfinity, fcZero):
861    return opOK;
862
863  case convolve(fcZero, fcNaN):
864  case convolve(fcNormal, fcNaN):
865  case convolve(fcInfinity, fcNaN):
866    category = fcNaN;
867    copySignificand(rhs);
868    return opOK;
869
870  case convolve(fcNormal, fcInfinity):
871  case convolve(fcZero, fcInfinity):
872    category = fcInfinity;
873    sign = rhs.sign ^ subtract;
874    return opOK;
875
876  case convolve(fcZero, fcNormal):
877    assign(rhs);
878    sign = rhs.sign ^ subtract;
879    return opOK;
880
881  case convolve(fcZero, fcZero):
882    /* Sign depends on rounding mode; handled by caller.  */
883    return opOK;
884
885  case convolve(fcInfinity, fcInfinity):
886    /* Differently signed infinities can only be validly
887       subtracted.  */
888    if(sign ^ rhs.sign != subtract) {
889      category = fcNaN;
890      // Arbitrary but deterministic value for significand
891      APInt::tcSet(significandParts(), ~0U, partCount());
892      return opInvalidOp;
893    }
894
895    return opOK;
896
897  case convolve(fcNormal, fcNormal):
898    return opDivByZero;
899  }
900}
901
902/* Add or subtract two normal numbers.  */
903lostFraction
904APFloat::addOrSubtractSignificand(const APFloat &rhs, bool subtract)
905{
906  integerPart carry;
907  lostFraction lost_fraction;
908  int bits;
909
910  /* Determine if the operation on the absolute values is effectively
911     an addition or subtraction.  */
912  subtract ^= (sign ^ rhs.sign);
913
914  /* Are we bigger exponent-wise than the RHS?  */
915  bits = exponent - rhs.exponent;
916
917  /* Subtraction is more subtle than one might naively expect.  */
918  if(subtract) {
919    APFloat temp_rhs(rhs);
920    bool reverse;
921
922    if (bits == 0) {
923      reverse = compareAbsoluteValue(temp_rhs) == cmpLessThan;
924      lost_fraction = lfExactlyZero;
925    } else if (bits > 0) {
926      lost_fraction = temp_rhs.shiftSignificandRight(bits - 1);
927      shiftSignificandLeft(1);
928      reverse = false;
929    } else {
930      lost_fraction = shiftSignificandRight(-bits - 1);
931      temp_rhs.shiftSignificandLeft(1);
932      reverse = true;
933    }
934
935    if (reverse) {
936      carry = temp_rhs.subtractSignificand
937	(*this, lost_fraction != lfExactlyZero);
938      copySignificand(temp_rhs);
939      sign = !sign;
940    } else {
941      carry = subtractSignificand
942	(temp_rhs, lost_fraction != lfExactlyZero);
943    }
944
945    /* Invert the lost fraction - it was on the RHS and
946       subtracted.  */
947    if(lost_fraction == lfLessThanHalf)
948      lost_fraction = lfMoreThanHalf;
949    else if(lost_fraction == lfMoreThanHalf)
950      lost_fraction = lfLessThanHalf;
951
952    /* The code above is intended to ensure that no borrow is
953       necessary.  */
954    assert(!carry);
955  } else {
956    if(bits > 0) {
957      APFloat temp_rhs(rhs);
958
959      lost_fraction = temp_rhs.shiftSignificandRight(bits);
960      carry = addSignificand(temp_rhs);
961    } else {
962      lost_fraction = shiftSignificandRight(-bits);
963      carry = addSignificand(rhs);
964    }
965
966    /* We have a guard bit; generating a carry cannot happen.  */
967    assert(!carry);
968  }
969
970  return lost_fraction;
971}
972
973APFloat::opStatus
974APFloat::multiplySpecials(const APFloat &rhs)
975{
976  switch(convolve(category, rhs.category)) {
977  default:
978    assert(0);
979
980  case convolve(fcNaN, fcZero):
981  case convolve(fcNaN, fcNormal):
982  case convolve(fcNaN, fcInfinity):
983  case convolve(fcNaN, fcNaN):
984    return opOK;
985
986  case convolve(fcZero, fcNaN):
987  case convolve(fcNormal, fcNaN):
988  case convolve(fcInfinity, fcNaN):
989    category = fcNaN;
990    copySignificand(rhs);
991    return opOK;
992
993  case convolve(fcNormal, fcInfinity):
994  case convolve(fcInfinity, fcNormal):
995  case convolve(fcInfinity, fcInfinity):
996    category = fcInfinity;
997    return opOK;
998
999  case convolve(fcZero, fcNormal):
1000  case convolve(fcNormal, fcZero):
1001  case convolve(fcZero, fcZero):
1002    category = fcZero;
1003    return opOK;
1004
1005  case convolve(fcZero, fcInfinity):
1006  case convolve(fcInfinity, fcZero):
1007    category = fcNaN;
1008    // Arbitrary but deterministic value for significand
1009    APInt::tcSet(significandParts(), ~0U, partCount());
1010    return opInvalidOp;
1011
1012  case convolve(fcNormal, fcNormal):
1013    return opOK;
1014  }
1015}
1016
1017APFloat::opStatus
1018APFloat::divideSpecials(const APFloat &rhs)
1019{
1020  switch(convolve(category, rhs.category)) {
1021  default:
1022    assert(0);
1023
1024  case convolve(fcNaN, fcZero):
1025  case convolve(fcNaN, fcNormal):
1026  case convolve(fcNaN, fcInfinity):
1027  case convolve(fcNaN, fcNaN):
1028  case convolve(fcInfinity, fcZero):
1029  case convolve(fcInfinity, fcNormal):
1030  case convolve(fcZero, fcInfinity):
1031  case convolve(fcZero, fcNormal):
1032    return opOK;
1033
1034  case convolve(fcZero, fcNaN):
1035  case convolve(fcNormal, fcNaN):
1036  case convolve(fcInfinity, fcNaN):
1037    category = fcNaN;
1038    copySignificand(rhs);
1039    return opOK;
1040
1041  case convolve(fcNormal, fcInfinity):
1042    category = fcZero;
1043    return opOK;
1044
1045  case convolve(fcNormal, fcZero):
1046    category = fcInfinity;
1047    return opDivByZero;
1048
1049  case convolve(fcInfinity, fcInfinity):
1050  case convolve(fcZero, fcZero):
1051    category = fcNaN;
1052    // Arbitrary but deterministic value for significand
1053    APInt::tcSet(significandParts(), ~0U, partCount());
1054    return opInvalidOp;
1055
1056  case convolve(fcNormal, fcNormal):
1057    return opOK;
1058  }
1059}
1060
1061/* Change sign.  */
1062void
1063APFloat::changeSign()
1064{
1065  /* Look mummy, this one's easy.  */
1066  sign = !sign;
1067}
1068
1069void
1070APFloat::clearSign()
1071{
1072  /* So is this one. */
1073  sign = 0;
1074}
1075
1076void
1077APFloat::copySign(const APFloat &rhs)
1078{
1079  /* And this one. */
1080  sign = rhs.sign;
1081}
1082
1083/* Normalized addition or subtraction.  */
1084APFloat::opStatus
1085APFloat::addOrSubtract(const APFloat &rhs, roundingMode rounding_mode,
1086		       bool subtract)
1087{
1088  opStatus fs;
1089
1090  fs = addOrSubtractSpecials(rhs, subtract);
1091
1092  /* This return code means it was not a simple case.  */
1093  if(fs == opDivByZero) {
1094    lostFraction lost_fraction;
1095
1096    lost_fraction = addOrSubtractSignificand(rhs, subtract);
1097    fs = normalize(rounding_mode, lost_fraction);
1098
1099    /* Can only be zero if we lost no fraction.  */
1100    assert(category != fcZero || lost_fraction == lfExactlyZero);
1101  }
1102
1103  /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1104     positive zero unless rounding to minus infinity, except that
1105     adding two like-signed zeroes gives that zero.  */
1106  if(category == fcZero) {
1107    if(rhs.category != fcZero || (sign == rhs.sign) == subtract)
1108      sign = (rounding_mode == rmTowardNegative);
1109  }
1110
1111  return fs;
1112}
1113
1114/* Normalized addition.  */
1115APFloat::opStatus
1116APFloat::add(const APFloat &rhs, roundingMode rounding_mode)
1117{
1118  return addOrSubtract(rhs, rounding_mode, false);
1119}
1120
1121/* Normalized subtraction.  */
1122APFloat::opStatus
1123APFloat::subtract(const APFloat &rhs, roundingMode rounding_mode)
1124{
1125  return addOrSubtract(rhs, rounding_mode, true);
1126}
1127
1128/* Normalized multiply.  */
1129APFloat::opStatus
1130APFloat::multiply(const APFloat &rhs, roundingMode rounding_mode)
1131{
1132  opStatus fs;
1133
1134  sign ^= rhs.sign;
1135  fs = multiplySpecials(rhs);
1136
1137  if(category == fcNormal) {
1138    lostFraction lost_fraction = multiplySignificand(rhs, 0);
1139    fs = normalize(rounding_mode, lost_fraction);
1140    if(lost_fraction != lfExactlyZero)
1141      fs = (opStatus) (fs | opInexact);
1142  }
1143
1144  return fs;
1145}
1146
1147/* Normalized divide.  */
1148APFloat::opStatus
1149APFloat::divide(const APFloat &rhs, roundingMode rounding_mode)
1150{
1151  opStatus fs;
1152
1153  sign ^= rhs.sign;
1154  fs = divideSpecials(rhs);
1155
1156  if(category == fcNormal) {
1157    lostFraction lost_fraction = divideSignificand(rhs);
1158    fs = normalize(rounding_mode, lost_fraction);
1159    if(lost_fraction != lfExactlyZero)
1160      fs = (opStatus) (fs | opInexact);
1161  }
1162
1163  return fs;
1164}
1165
1166/* Normalized remainder. */
1167APFloat::opStatus
1168APFloat::mod(const APFloat &rhs, roundingMode rounding_mode)
1169{
1170  opStatus fs;
1171  APFloat V = *this;
1172  unsigned int origSign = sign;
1173  fs = V.divide(rhs, rmNearestTiesToEven);
1174  if (fs == opDivByZero)
1175    return fs;
1176
1177  int parts = partCount();
1178  integerPart *x = new integerPart[parts];
1179  fs = V.convertToInteger(x, parts * integerPartWidth, true,
1180                          rmNearestTiesToEven);
1181  if (fs==opInvalidOp)
1182    return fs;
1183
1184  fs = V.convertFromInteger(x, parts, true, 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  unsigned int newPartCount;
1322  opStatus fs;
1323
1324  newPartCount = partCountForBits(toSemantics.precision + 1);
1325
1326  /* If our new form is wider, re-allocate our bit pattern into wider
1327     storage.  */
1328  if(newPartCount > partCount()) {
1329    integerPart *newParts;
1330
1331    newParts = new integerPart[newPartCount];
1332    APInt::tcSet(newParts, 0, newPartCount);
1333    APInt::tcAssign(newParts, significandParts(), partCount());
1334    freeSignificand();
1335    significand.parts = newParts;
1336  }
1337
1338  if(category == fcNormal) {
1339    /* Re-interpret our bit-pattern.  */
1340    exponent += toSemantics.precision - semantics->precision;
1341    semantics = &toSemantics;
1342    fs = normalize(rounding_mode, lfExactlyZero);
1343  } else {
1344    semantics = &toSemantics;
1345    fs = opOK;
1346  }
1347
1348  return fs;
1349}
1350
1351/* Convert a floating point number to an integer according to the
1352   rounding mode.  If the rounded integer value is out of range this
1353   returns an invalid operation exception.  If the rounded value is in
1354   range but the floating point number is not the exact integer, the C
1355   standard doesn't require an inexact exception to be raised.  IEEE
1356   854 does require it so we do that.
1357
1358   Note that for conversions to integer type the C standard requires
1359   round-to-zero to always be used.  */
1360APFloat::opStatus
1361APFloat::convertToInteger(integerPart *parts, unsigned int width,
1362			  bool isSigned,
1363			  roundingMode rounding_mode) const
1364{
1365  lostFraction lost_fraction;
1366  unsigned int msb, partsCount;
1367  int bits;
1368
1369  /* Handle the three special cases first.  */
1370  if(category == fcInfinity || category == fcNaN)
1371    return opInvalidOp;
1372
1373  partsCount = partCountForBits(width);
1374
1375  if(category == fcZero) {
1376    APInt::tcSet(parts, 0, partsCount);
1377    return opOK;
1378  }
1379
1380  /* Shift the bit pattern so the fraction is lost.  */
1381  APFloat tmp(*this);
1382
1383  bits = (int) semantics->precision - 1 - exponent;
1384
1385  if(bits > 0) {
1386    lost_fraction = tmp.shiftSignificandRight(bits);
1387  } else {
1388    tmp.shiftSignificandLeft(-bits);
1389    lost_fraction = lfExactlyZero;
1390  }
1391
1392  if(lost_fraction != lfExactlyZero
1393     && tmp.roundAwayFromZero(rounding_mode, lost_fraction))
1394    tmp.incrementSignificand();
1395
1396  msb = tmp.significandMSB();
1397
1398  /* Negative numbers cannot be represented as unsigned.  */
1399  if(!isSigned && tmp.sign && msb != -1U)
1400    return opInvalidOp;
1401
1402  /* It takes exponent + 1 bits to represent the truncated floating
1403     point number without its sign.  We lose a bit for the sign, but
1404     the maximally negative integer is a special case.  */
1405  if(msb + 1 > width)		/* !! Not same as msb >= width !! */
1406    return opInvalidOp;
1407
1408  if(isSigned && msb + 1 == width
1409     && (!tmp.sign || tmp.significandLSB() != msb))
1410    return opInvalidOp;
1411
1412  APInt::tcAssign(parts, tmp.significandParts(), partsCount);
1413
1414  if(tmp.sign)
1415    APInt::tcNegate(parts, partsCount);
1416
1417  if(lost_fraction == lfExactlyZero)
1418    return opOK;
1419  else
1420    return opInexact;
1421}
1422
1423APFloat::opStatus
1424APFloat::convertFromUnsignedInteger(integerPart *parts,
1425				    unsigned int partCount,
1426				    roundingMode rounding_mode)
1427{
1428  unsigned int msb, precision;
1429  lostFraction lost_fraction;
1430
1431  msb = APInt::tcMSB(parts, partCount) + 1;
1432  precision = semantics->precision;
1433
1434  category = fcNormal;
1435  exponent = precision - 1;
1436
1437  if(msb > precision) {
1438    exponent += (msb - precision);
1439    lost_fraction = shiftRight(parts, partCount, msb - precision);
1440    msb = precision;
1441  } else
1442    lost_fraction = lfExactlyZero;
1443
1444  /* Copy the bit image.  */
1445  zeroSignificand();
1446  APInt::tcAssign(significandParts(), parts, partCountForBits(msb));
1447
1448  return normalize(rounding_mode, lost_fraction);
1449}
1450
1451APFloat::opStatus
1452APFloat::convertFromInteger(const integerPart *parts,
1453			    unsigned int partCount, bool isSigned,
1454			    roundingMode rounding_mode)
1455{
1456  unsigned int width;
1457  opStatus status;
1458  integerPart *copy;
1459
1460  copy = new integerPart[partCount];
1461  APInt::tcAssign(copy, parts, partCount);
1462
1463  width = partCount * integerPartWidth;
1464
1465  sign = false;
1466  if(isSigned && APInt::tcExtractBit(parts, width - 1)) {
1467    sign = true;
1468    APInt::tcNegate(copy, partCount);
1469  }
1470
1471  status = convertFromUnsignedInteger(copy, partCount, rounding_mode);
1472  delete [] copy;
1473
1474  return status;
1475}
1476
1477APFloat::opStatus
1478APFloat::convertFromHexadecimalString(const char *p,
1479				      roundingMode rounding_mode)
1480{
1481  lostFraction lost_fraction;
1482  integerPart *significand;
1483  unsigned int bitPos, partsCount;
1484  const char *dot, *firstSignificantDigit;
1485
1486  zeroSignificand();
1487  exponent = 0;
1488  category = fcNormal;
1489
1490  significand = significandParts();
1491  partsCount = partCount();
1492  bitPos = partsCount * integerPartWidth;
1493
1494  /* Skip leading zeroes and any(hexa)decimal point.  */
1495  p = skipLeadingZeroesAndAnyDot(p, &dot);
1496  firstSignificantDigit = p;
1497
1498  for(;;) {
1499    integerPart hex_value;
1500
1501    if(*p == '.') {
1502      assert(dot == 0);
1503      dot = p++;
1504    }
1505
1506    hex_value = hexDigitValue(*p);
1507    if(hex_value == -1U) {
1508      lost_fraction = lfExactlyZero;
1509      break;
1510    }
1511
1512    p++;
1513
1514    /* Store the number whilst 4-bit nibbles remain.  */
1515    if(bitPos) {
1516      bitPos -= 4;
1517      hex_value <<= bitPos % integerPartWidth;
1518      significand[bitPos / integerPartWidth] |= hex_value;
1519    } else {
1520      lost_fraction = trailingHexadecimalFraction(p, hex_value);
1521      while(hexDigitValue(*p) != -1U)
1522	p++;
1523      break;
1524    }
1525  }
1526
1527  /* Hex floats require an exponent but not a hexadecimal point.  */
1528  assert(*p == 'p' || *p == 'P');
1529
1530  /* Ignore the exponent if we are zero.  */
1531  if(p != firstSignificantDigit) {
1532    int expAdjustment;
1533
1534    /* Implicit hexadecimal point?  */
1535    if(!dot)
1536      dot = p;
1537
1538    /* Calculate the exponent adjustment implicit in the number of
1539       significant digits.  */
1540    expAdjustment = dot - firstSignificantDigit;
1541    if(expAdjustment < 0)
1542      expAdjustment++;
1543    expAdjustment = expAdjustment * 4 - 1;
1544
1545    /* Adjust for writing the significand starting at the most
1546       significant nibble.  */
1547    expAdjustment += semantics->precision;
1548    expAdjustment -= partsCount * integerPartWidth;
1549
1550    /* Adjust for the given exponent.  */
1551    exponent = totalExponent(p, expAdjustment);
1552  }
1553
1554  return normalize(rounding_mode, lost_fraction);
1555}
1556
1557APFloat::opStatus
1558APFloat::convertFromString(const char *p, roundingMode rounding_mode) {
1559  /* Handle a leading minus sign.  */
1560  if(*p == '-')
1561    sign = 1, p++;
1562  else
1563    sign = 0;
1564
1565  if(p[0] == '0' && (p[1] == 'x' || p[1] == 'X'))
1566    return convertFromHexadecimalString(p + 2, rounding_mode);
1567
1568  assert(0 && "Decimal to binary conversions not yet implemented");
1569  abort();
1570}
1571
1572// For good performance it is desirable for different APFloats
1573// to produce different integers.
1574uint32_t
1575APFloat::getHashValue() const {
1576  if (category==fcZero) return sign<<8 | semantics->precision ;
1577  else if (category==fcInfinity) return sign<<9 | semantics->precision;
1578  else if (category==fcNaN) return 1<<10 | semantics->precision;
1579  else {
1580    uint32_t hash = sign<<11 | semantics->precision | exponent<<12;
1581    const integerPart* p = significandParts();
1582    for (int i=partCount(); i>0; i--, p++)
1583      hash ^= ((uint32_t)*p) ^ (*p)>>32;
1584    return hash;
1585  }
1586}
1587
1588// Conversion from APFloat to/from host float/double.  It may eventually be
1589// possible to eliminate these and have everybody deal with APFloats, but that
1590// will take a while.  This approach will not easily extend to long double.
1591// Current implementation requires partCount()==1, which is correct at the
1592// moment but could be made more general.
1593
1594// Denormals have exponent minExponent in APFloat, but minExponent-1 in
1595// the actual IEEE respresentation.  We compensate for that here.
1596
1597APInt
1598APFloat::convertF80LongDoubleAPFloatToAPInt() const {
1599  assert(semantics == (const llvm::fltSemantics* const)&x87DoubleExtended);
1600  assert (partCount()==1);
1601
1602  uint64_t myexponent, mysignificand;
1603
1604  if (category==fcNormal) {
1605    myexponent = exponent+16383; //bias
1606    mysignificand = *significandParts();
1607    if (myexponent==1 && !(mysignificand & 0x8000000000000000ULL))
1608      myexponent = 0;   // denormal
1609  } else if (category==fcZero) {
1610    myexponent = 0;
1611    mysignificand = 0;
1612  } else if (category==fcInfinity) {
1613    myexponent = 0x7fff;
1614    mysignificand = 0x8000000000000000ULL;
1615  } else if (category==fcNaN) {
1616    myexponent = 0x7fff;
1617    mysignificand = *significandParts();
1618  } else
1619    assert(0);
1620
1621  uint64_t words[2];
1622  words[0] =  (((uint64_t)sign & 1) << 63) |
1623              ((myexponent & 0x7fff) <<  48) |
1624              ((mysignificand >>16) & 0xffffffffffffLL);
1625  words[1] = mysignificand & 0xffff;
1626  APInt api(80, 2, words);
1627  return api;
1628}
1629
1630APInt
1631APFloat::convertDoubleAPFloatToAPInt() const {
1632  assert(semantics == (const llvm::fltSemantics* const)&IEEEdouble);
1633  assert (partCount()==1);
1634
1635  uint64_t myexponent, mysignificand;
1636
1637  if (category==fcNormal) {
1638    myexponent = exponent+1023; //bias
1639    mysignificand = *significandParts();
1640    if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
1641      myexponent = 0;   // denormal
1642  } else if (category==fcZero) {
1643    myexponent = 0;
1644    mysignificand = 0;
1645  } else if (category==fcInfinity) {
1646    myexponent = 0x7ff;
1647    mysignificand = 0;
1648  } else if (category==fcNaN) {
1649    myexponent = 0x7ff;
1650    mysignificand = *significandParts();
1651  } else
1652    assert(0);
1653
1654  APInt api(64, (((((uint64_t)sign & 1) << 63) |
1655                 ((myexponent & 0x7ff) <<  52) |
1656                 (mysignificand & 0xfffffffffffffLL))));
1657  return api;
1658}
1659
1660APInt
1661APFloat::convertFloatAPFloatToAPInt() const {
1662  assert(semantics == (const llvm::fltSemantics* const)&IEEEsingle);
1663  assert (partCount()==1);
1664
1665  uint32_t myexponent, mysignificand;
1666
1667  if (category==fcNormal) {
1668    myexponent = exponent+127; //bias
1669    mysignificand = *significandParts();
1670    if (myexponent == 1 && !(mysignificand & 0x400000))
1671      myexponent = 0;   // denormal
1672  } else if (category==fcZero) {
1673    myexponent = 0;
1674    mysignificand = 0;
1675  } else if (category==fcInfinity) {
1676    myexponent = 0xff;
1677    mysignificand = 0;
1678  } else if (category==fcNaN) {
1679    myexponent = 0xff;
1680    mysignificand = *significandParts();
1681  } else
1682    assert(0);
1683
1684  APInt api(32, (((sign&1) << 31) | ((myexponent&0xff) << 23) |
1685                 (mysignificand & 0x7fffff)));
1686  return api;
1687}
1688
1689APInt
1690APFloat::convertToAPInt() const {
1691  if (semantics == (const llvm::fltSemantics* const)&IEEEsingle)
1692    return convertFloatAPFloatToAPInt();
1693  else if (semantics == (const llvm::fltSemantics* const)&IEEEdouble)
1694    return convertDoubleAPFloatToAPInt();
1695  else if (semantics == (const llvm::fltSemantics* const)&x87DoubleExtended)
1696    return convertF80LongDoubleAPFloatToAPInt();
1697  else
1698    assert(0);
1699}
1700
1701float
1702APFloat::convertToFloat() const {
1703  assert(semantics == (const llvm::fltSemantics* const)&IEEEsingle);
1704  APInt api = convertToAPInt();
1705  return api.bitsToFloat();
1706}
1707
1708double
1709APFloat::convertToDouble() const {
1710  assert(semantics == (const llvm::fltSemantics* const)&IEEEdouble);
1711  APInt api = convertToAPInt();
1712  return api.bitsToDouble();
1713}
1714
1715/// Integer bit is explicit in this format.  Current Intel book does not
1716/// define meaning of:
1717///  exponent = all 1's, integer bit not set.
1718///  exponent = 0, integer bit set. (formerly "psuedodenormals")
1719///  exponent!=0 nor all 1's, integer bit not set. (formerly "unnormals")
1720void
1721APFloat::initFromF80LongDoubleAPInt(const APInt &api) {
1722  assert(api.getBitWidth()==80);
1723  uint64_t i1 = api.getRawData()[0];
1724  uint64_t i2 = api.getRawData()[1];
1725  uint64_t myexponent = (i1 >> 48) & 0x7fff;
1726  uint64_t mysignificand = ((i1 << 16) &  0xffffffffffff0000ULL) |
1727                          (i2 & 0xffff);
1728
1729  initialize(&APFloat::x87DoubleExtended);
1730  assert(partCount()==1);
1731
1732  sign = i1>>63;
1733  if (myexponent==0 && mysignificand==0) {
1734    // exponent, significand meaningless
1735    category = fcZero;
1736  } else if (myexponent==0x7fff && mysignificand==0x8000000000000000ULL) {
1737    // exponent, significand meaningless
1738    category = fcInfinity;
1739  } else if (myexponent==0x7fff && mysignificand!=0x8000000000000000ULL) {
1740    // exponent meaningless
1741    category = fcNaN;
1742    *significandParts() = mysignificand;
1743  } else {
1744    category = fcNormal;
1745    exponent = myexponent - 16383;
1746    *significandParts() = mysignificand;
1747    if (myexponent==0)          // denormal
1748      exponent = -16382;
1749 }
1750}
1751
1752void
1753APFloat::initFromDoubleAPInt(const APInt &api) {
1754  assert(api.getBitWidth()==64);
1755  uint64_t i = *api.getRawData();
1756  uint64_t myexponent = (i >> 52) & 0x7ff;
1757  uint64_t mysignificand = i & 0xfffffffffffffLL;
1758
1759  initialize(&APFloat::IEEEdouble);
1760  assert(partCount()==1);
1761
1762  sign = i>>63;
1763  if (myexponent==0 && mysignificand==0) {
1764    // exponent, significand meaningless
1765    category = fcZero;
1766  } else if (myexponent==0x7ff && mysignificand==0) {
1767    // exponent, significand meaningless
1768    category = fcInfinity;
1769  } else if (myexponent==0x7ff && mysignificand!=0) {
1770    // exponent meaningless
1771    category = fcNaN;
1772    *significandParts() = mysignificand;
1773  } else {
1774    category = fcNormal;
1775    exponent = myexponent - 1023;
1776    *significandParts() = mysignificand;
1777    if (myexponent==0)          // denormal
1778      exponent = -1022;
1779    else
1780      *significandParts() |= 0x10000000000000LL;  // integer bit
1781 }
1782}
1783
1784void
1785APFloat::initFromFloatAPInt(const APInt & api) {
1786  assert(api.getBitWidth()==32);
1787  uint32_t i = (uint32_t)*api.getRawData();
1788  uint32_t myexponent = (i >> 23) & 0xff;
1789  uint32_t mysignificand = i & 0x7fffff;
1790
1791  initialize(&APFloat::IEEEsingle);
1792  assert(partCount()==1);
1793
1794  sign = i >> 31;
1795  if (myexponent==0 && mysignificand==0) {
1796    // exponent, significand meaningless
1797    category = fcZero;
1798  } else if (myexponent==0xff && mysignificand==0) {
1799    // exponent, significand meaningless
1800    category = fcInfinity;
1801  } else if (myexponent==0xff && (mysignificand & 0x400000)) {
1802    // sign, exponent, significand meaningless
1803    category = fcNaN;
1804    *significandParts() = mysignificand;
1805  } else {
1806    category = fcNormal;
1807    exponent = myexponent - 127;  //bias
1808    *significandParts() = mysignificand;
1809    if (myexponent==0)    // denormal
1810      exponent = -126;
1811    else
1812      *significandParts() |= 0x800000; // integer bit
1813  }
1814}
1815
1816/// Treat api as containing the bits of a floating point number.  Currently
1817/// we infer the floating point type from the size of the APInt.  FIXME: This
1818/// breaks when we get to PPC128 and IEEE128 (but both cannot exist in the
1819/// same compile...)
1820void
1821APFloat::initFromAPInt(const APInt& api) {
1822  if (api.getBitWidth() == 32)
1823    return initFromFloatAPInt(api);
1824  else if (api.getBitWidth()==64)
1825    return initFromDoubleAPInt(api);
1826  else if (api.getBitWidth()==80)
1827    return initFromF80LongDoubleAPInt(api);
1828  else
1829    assert(0);
1830}
1831
1832APFloat::APFloat(const APInt& api) {
1833  initFromAPInt(api);
1834}
1835
1836APFloat::APFloat(float f) {
1837  APInt api = APInt(32, 0);
1838  initFromAPInt(api.floatToBits(f));
1839}
1840
1841APFloat::APFloat(double d) {
1842  APInt api = APInt(64, 0);
1843  initFromAPInt(api.doubleToBits(d));
1844}
1845
1846