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