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/APSInt.h"
17#include "llvm/ADT/FoldingSet.h"
18#include "llvm/ADT/Hashing.h"
19#include "llvm/ADT/StringExtras.h"
20#include "llvm/ADT/StringRef.h"
21#include "llvm/Support/ErrorHandling.h"
22#include "llvm/Support/MathExtras.h"
23#include <cstring>
24#include <limits.h>
25
26using namespace llvm;
27
28/// A macro used to combine two fcCategory enums into one key which can be used
29/// in a switch statement to classify how the interaction of two APFloat's
30/// categories affects an operation.
31///
32/// TODO: If clang source code is ever allowed to use constexpr in its own
33/// codebase, change this into a static inline function.
34#define PackCategoriesIntoKey(_lhs, _rhs) ((_lhs) * 4 + (_rhs))
35
36/* Assumed in hexadecimal significand parsing, and conversion to
37   hexadecimal strings.  */
38#define COMPILE_TIME_ASSERT(cond) extern int CTAssert[(cond) ? 1 : -1]
39COMPILE_TIME_ASSERT(integerPartWidth % 4 == 0);
40
41namespace llvm {
42
43  /* Represents floating point arithmetic semantics.  */
44  struct fltSemantics {
45    /* The largest E such that 2^E is representable; this matches the
46       definition of IEEE 754.  */
47    APFloat::ExponentType maxExponent;
48
49    /* The smallest E such that 2^E is a normalized number; this
50       matches the definition of IEEE 754.  */
51    APFloat::ExponentType minExponent;
52
53    /* Number of bits in the significand.  This includes the integer
54       bit.  */
55    unsigned int precision;
56  };
57
58  const fltSemantics APFloat::IEEEhalf = { 15, -14, 11 };
59  const fltSemantics APFloat::IEEEsingle = { 127, -126, 24 };
60  const fltSemantics APFloat::IEEEdouble = { 1023, -1022, 53 };
61  const fltSemantics APFloat::IEEEquad = { 16383, -16382, 113 };
62  const fltSemantics APFloat::x87DoubleExtended = { 16383, -16382, 64 };
63  const fltSemantics APFloat::Bogus = { 0, 0, 0 };
64
65  /* The PowerPC format consists of two doubles.  It does not map cleanly
66     onto the usual format above.  It is approximated using twice the
67     mantissa bits.  Note that for exponents near the double minimum,
68     we no longer can represent the full 106 mantissa bits, so those
69     will be treated as denormal numbers.
70
71     FIXME: While this approximation is equivalent to what GCC uses for
72     compile-time arithmetic on PPC double-double numbers, it is not able
73     to represent all possible values held by a PPC double-double number,
74     for example: (long double) 1.0 + (long double) 0x1p-106
75     Should this be replaced by a full emulation of PPC double-double?  */
76  const fltSemantics APFloat::PPCDoubleDouble = { 1023, -1022 + 53, 53 + 53 };
77
78  /* A tight upper bound on number of parts required to hold the value
79     pow(5, power) is
80
81       power * 815 / (351 * integerPartWidth) + 1
82
83     However, whilst the result may require only this many parts,
84     because we are multiplying two values to get it, the
85     multiplication may require an extra part with the excess part
86     being zero (consider the trivial case of 1 * 1, tcFullMultiply
87     requires two parts to hold the single-part result).  So we add an
88     extra one to guarantee enough space whilst multiplying.  */
89  const unsigned int maxExponent = 16383;
90  const unsigned int maxPrecision = 113;
91  const unsigned int maxPowerOfFiveExponent = maxExponent + maxPrecision - 1;
92  const unsigned int maxPowerOfFiveParts = 2 + ((maxPowerOfFiveExponent * 815)
93                                                / (351 * integerPartWidth));
94}
95
96/* A bunch of private, handy routines.  */
97
98static inline unsigned int
99partCountForBits(unsigned int bits)
100{
101  return ((bits) + integerPartWidth - 1) / integerPartWidth;
102}
103
104/* Returns 0U-9U.  Return values >= 10U are not digits.  */
105static inline unsigned int
106decDigitValue(unsigned int c)
107{
108  return c - '0';
109}
110
111/* Return the value of a decimal exponent of the form
112   [+-]ddddddd.
113
114   If the exponent overflows, returns a large exponent with the
115   appropriate sign.  */
116static int
117readExponent(StringRef::iterator begin, StringRef::iterator end)
118{
119  bool isNegative;
120  unsigned int absExponent;
121  const unsigned int overlargeExponent = 24000;  /* FIXME.  */
122  StringRef::iterator p = begin;
123
124  assert(p != end && "Exponent has no digits");
125
126  isNegative = (*p == '-');
127  if (*p == '-' || *p == '+') {
128    p++;
129    assert(p != end && "Exponent has no digits");
130  }
131
132  absExponent = decDigitValue(*p++);
133  assert(absExponent < 10U && "Invalid character in exponent");
134
135  for (; p != end; ++p) {
136    unsigned int value;
137
138    value = decDigitValue(*p);
139    assert(value < 10U && "Invalid character in exponent");
140
141    value += absExponent * 10;
142    if (absExponent >= overlargeExponent) {
143      absExponent = overlargeExponent;
144      p = end;  /* outwit assert below */
145      break;
146    }
147    absExponent = value;
148  }
149
150  assert(p == end && "Invalid exponent in exponent");
151
152  if (isNegative)
153    return -(int) absExponent;
154  else
155    return (int) absExponent;
156}
157
158/* This is ugly and needs cleaning up, but I don't immediately see
159   how whilst remaining safe.  */
160static int
161totalExponent(StringRef::iterator p, StringRef::iterator end,
162              int exponentAdjustment)
163{
164  int unsignedExponent;
165  bool negative, overflow;
166  int exponent = 0;
167
168  assert(p != end && "Exponent has no digits");
169
170  negative = *p == '-';
171  if (*p == '-' || *p == '+') {
172    p++;
173    assert(p != end && "Exponent has no digits");
174  }
175
176  unsignedExponent = 0;
177  overflow = false;
178  for (; p != end; ++p) {
179    unsigned int value;
180
181    value = decDigitValue(*p);
182    assert(value < 10U && "Invalid character in exponent");
183
184    unsignedExponent = unsignedExponent * 10 + value;
185    if (unsignedExponent > 32767) {
186      overflow = true;
187      break;
188    }
189  }
190
191  if (exponentAdjustment > 32767 || exponentAdjustment < -32768)
192    overflow = true;
193
194  if (!overflow) {
195    exponent = unsignedExponent;
196    if (negative)
197      exponent = -exponent;
198    exponent += exponentAdjustment;
199    if (exponent > 32767 || exponent < -32768)
200      overflow = true;
201  }
202
203  if (overflow)
204    exponent = negative ? -32768: 32767;
205
206  return exponent;
207}
208
209static StringRef::iterator
210skipLeadingZeroesAndAnyDot(StringRef::iterator begin, StringRef::iterator end,
211                           StringRef::iterator *dot)
212{
213  StringRef::iterator p = begin;
214  *dot = end;
215  while (*p == '0' && p != end)
216    p++;
217
218  if (*p == '.') {
219    *dot = p++;
220
221    assert(end - begin != 1 && "Significand has no digits");
222
223    while (*p == '0' && p != end)
224      p++;
225  }
226
227  return p;
228}
229
230/* Given a normal decimal floating point number of the form
231
232     dddd.dddd[eE][+-]ddd
233
234   where the decimal point and exponent are optional, fill out the
235   structure D.  Exponent is appropriate if the significand is
236   treated as an integer, and normalizedExponent if the significand
237   is taken to have the decimal point after a single leading
238   non-zero digit.
239
240   If the value is zero, V->firstSigDigit points to a non-digit, and
241   the return exponent is zero.
242*/
243struct decimalInfo {
244  const char *firstSigDigit;
245  const char *lastSigDigit;
246  int exponent;
247  int normalizedExponent;
248};
249
250static void
251interpretDecimal(StringRef::iterator begin, StringRef::iterator end,
252                 decimalInfo *D)
253{
254  StringRef::iterator dot = end;
255  StringRef::iterator p = skipLeadingZeroesAndAnyDot (begin, end, &dot);
256
257  D->firstSigDigit = p;
258  D->exponent = 0;
259  D->normalizedExponent = 0;
260
261  for (; p != end; ++p) {
262    if (*p == '.') {
263      assert(dot == end && "String contains multiple dots");
264      dot = p++;
265      if (p == end)
266        break;
267    }
268    if (decDigitValue(*p) >= 10U)
269      break;
270  }
271
272  if (p != end) {
273    assert((*p == 'e' || *p == 'E') && "Invalid character in significand");
274    assert(p != begin && "Significand has no digits");
275    assert((dot == end || p - begin != 1) && "Significand has no digits");
276
277    /* p points to the first non-digit in the string */
278    D->exponent = readExponent(p + 1, end);
279
280    /* Implied decimal point?  */
281    if (dot == end)
282      dot = p;
283  }
284
285  /* If number is all zeroes accept any exponent.  */
286  if (p != D->firstSigDigit) {
287    /* Drop insignificant trailing zeroes.  */
288    if (p != begin) {
289      do
290        do
291          p--;
292        while (p != begin && *p == '0');
293      while (p != begin && *p == '.');
294    }
295
296    /* Adjust the exponents for any decimal point.  */
297    D->exponent += static_cast<APFloat::ExponentType>((dot - p) - (dot > p));
298    D->normalizedExponent = (D->exponent +
299              static_cast<APFloat::ExponentType>((p - D->firstSigDigit)
300                                      - (dot > D->firstSigDigit && dot < p)));
301  }
302
303  D->lastSigDigit = p;
304}
305
306/* Return the trailing fraction of a hexadecimal number.
307   DIGITVALUE is the first hex digit of the fraction, P points to
308   the next digit.  */
309static lostFraction
310trailingHexadecimalFraction(StringRef::iterator p, StringRef::iterator end,
311                            unsigned int digitValue)
312{
313  unsigned int hexDigit;
314
315  /* If the first trailing digit isn't 0 or 8 we can work out the
316     fraction immediately.  */
317  if (digitValue > 8)
318    return lfMoreThanHalf;
319  else if (digitValue < 8 && digitValue > 0)
320    return lfLessThanHalf;
321
322  // Otherwise we need to find the first non-zero digit.
323  while (p != end && (*p == '0' || *p == '.'))
324    p++;
325
326  assert(p != end && "Invalid trailing hexadecimal fraction!");
327
328  hexDigit = hexDigitValue(*p);
329
330  /* If we ran off the end it is exactly zero or one-half, otherwise
331     a little more.  */
332  if (hexDigit == -1U)
333    return digitValue == 0 ? lfExactlyZero: lfExactlyHalf;
334  else
335    return digitValue == 0 ? lfLessThanHalf: lfMoreThanHalf;
336}
337
338/* Return the fraction lost were a bignum truncated losing the least
339   significant BITS bits.  */
340static lostFraction
341lostFractionThroughTruncation(const integerPart *parts,
342                              unsigned int partCount,
343                              unsigned int bits)
344{
345  unsigned int lsb;
346
347  lsb = APInt::tcLSB(parts, partCount);
348
349  /* Note this is guaranteed true if bits == 0, or LSB == -1U.  */
350  if (bits <= lsb)
351    return lfExactlyZero;
352  if (bits == lsb + 1)
353    return lfExactlyHalf;
354  if (bits <= partCount * integerPartWidth &&
355      APInt::tcExtractBit(parts, bits - 1))
356    return lfMoreThanHalf;
357
358  return lfLessThanHalf;
359}
360
361/* Shift DST right BITS bits noting lost fraction.  */
362static lostFraction
363shiftRight(integerPart *dst, unsigned int parts, unsigned int bits)
364{
365  lostFraction lost_fraction;
366
367  lost_fraction = lostFractionThroughTruncation(dst, parts, bits);
368
369  APInt::tcShiftRight(dst, parts, bits);
370
371  return lost_fraction;
372}
373
374/* Combine the effect of two lost fractions.  */
375static lostFraction
376combineLostFractions(lostFraction moreSignificant,
377                     lostFraction lessSignificant)
378{
379  if (lessSignificant != lfExactlyZero) {
380    if (moreSignificant == lfExactlyZero)
381      moreSignificant = lfLessThanHalf;
382    else if (moreSignificant == lfExactlyHalf)
383      moreSignificant = lfMoreThanHalf;
384  }
385
386  return moreSignificant;
387}
388
389/* The error from the true value, in half-ulps, on multiplying two
390   floating point numbers, which differ from the value they
391   approximate by at most HUE1 and HUE2 half-ulps, is strictly less
392   than the returned value.
393
394   See "How to Read Floating Point Numbers Accurately" by William D
395   Clinger.  */
396static unsigned int
397HUerrBound(bool inexactMultiply, unsigned int HUerr1, unsigned int HUerr2)
398{
399  assert(HUerr1 < 2 || HUerr2 < 2 || (HUerr1 + HUerr2 < 8));
400
401  if (HUerr1 + HUerr2 == 0)
402    return inexactMultiply * 2;  /* <= inexactMultiply half-ulps.  */
403  else
404    return inexactMultiply + 2 * (HUerr1 + HUerr2);
405}
406
407/* The number of ulps from the boundary (zero, or half if ISNEAREST)
408   when the least significant BITS are truncated.  BITS cannot be
409   zero.  */
410static integerPart
411ulpsFromBoundary(const integerPart *parts, unsigned int bits, bool isNearest)
412{
413  unsigned int count, partBits;
414  integerPart part, boundary;
415
416  assert(bits != 0);
417
418  bits--;
419  count = bits / integerPartWidth;
420  partBits = bits % integerPartWidth + 1;
421
422  part = parts[count] & (~(integerPart) 0 >> (integerPartWidth - partBits));
423
424  if (isNearest)
425    boundary = (integerPart) 1 << (partBits - 1);
426  else
427    boundary = 0;
428
429  if (count == 0) {
430    if (part - boundary <= boundary - part)
431      return part - boundary;
432    else
433      return boundary - part;
434  }
435
436  if (part == boundary) {
437    while (--count)
438      if (parts[count])
439        return ~(integerPart) 0; /* A lot.  */
440
441    return parts[0];
442  } else if (part == boundary - 1) {
443    while (--count)
444      if (~parts[count])
445        return ~(integerPart) 0; /* A lot.  */
446
447    return -parts[0];
448  }
449
450  return ~(integerPart) 0; /* A lot.  */
451}
452
453/* Place pow(5, power) in DST, and return the number of parts used.
454   DST must be at least one part larger than size of the answer.  */
455static unsigned int
456powerOf5(integerPart *dst, unsigned int power)
457{
458  static const integerPart firstEightPowers[] = { 1, 5, 25, 125, 625, 3125,
459                                                  15625, 78125 };
460  integerPart pow5s[maxPowerOfFiveParts * 2 + 5];
461  pow5s[0] = 78125 * 5;
462
463  unsigned int partsCount[16] = { 1 };
464  integerPart scratch[maxPowerOfFiveParts], *p1, *p2, *pow5;
465  unsigned int result;
466  assert(power <= maxExponent);
467
468  p1 = dst;
469  p2 = scratch;
470
471  *p1 = firstEightPowers[power & 7];
472  power >>= 3;
473
474  result = 1;
475  pow5 = pow5s;
476
477  for (unsigned int n = 0; power; power >>= 1, n++) {
478    unsigned int pc;
479
480    pc = partsCount[n];
481
482    /* Calculate pow(5,pow(2,n+3)) if we haven't yet.  */
483    if (pc == 0) {
484      pc = partsCount[n - 1];
485      APInt::tcFullMultiply(pow5, pow5 - pc, pow5 - pc, pc, pc);
486      pc *= 2;
487      if (pow5[pc - 1] == 0)
488        pc--;
489      partsCount[n] = pc;
490    }
491
492    if (power & 1) {
493      integerPart *tmp;
494
495      APInt::tcFullMultiply(p2, p1, pow5, result, pc);
496      result += pc;
497      if (p2[result - 1] == 0)
498        result--;
499
500      /* Now result is in p1 with partsCount parts and p2 is scratch
501         space.  */
502      tmp = p1, p1 = p2, p2 = tmp;
503    }
504
505    pow5 += pc;
506  }
507
508  if (p1 != dst)
509    APInt::tcAssign(dst, p1, result);
510
511  return result;
512}
513
514/* Zero at the end to avoid modular arithmetic when adding one; used
515   when rounding up during hexadecimal output.  */
516static const char hexDigitsLower[] = "0123456789abcdef0";
517static const char hexDigitsUpper[] = "0123456789ABCDEF0";
518static const char infinityL[] = "infinity";
519static const char infinityU[] = "INFINITY";
520static const char NaNL[] = "nan";
521static const char NaNU[] = "NAN";
522
523/* Write out an integerPart in hexadecimal, starting with the most
524   significant nibble.  Write out exactly COUNT hexdigits, return
525   COUNT.  */
526static unsigned int
527partAsHex (char *dst, integerPart part, unsigned int count,
528           const char *hexDigitChars)
529{
530  unsigned int result = count;
531
532  assert(count != 0 && count <= integerPartWidth / 4);
533
534  part >>= (integerPartWidth - 4 * count);
535  while (count--) {
536    dst[count] = hexDigitChars[part & 0xf];
537    part >>= 4;
538  }
539
540  return result;
541}
542
543/* Write out an unsigned decimal integer.  */
544static char *
545writeUnsignedDecimal (char *dst, unsigned int n)
546{
547  char buff[40], *p;
548
549  p = buff;
550  do
551    *p++ = '0' + n % 10;
552  while (n /= 10);
553
554  do
555    *dst++ = *--p;
556  while (p != buff);
557
558  return dst;
559}
560
561/* Write out a signed decimal integer.  */
562static char *
563writeSignedDecimal (char *dst, int value)
564{
565  if (value < 0) {
566    *dst++ = '-';
567    dst = writeUnsignedDecimal(dst, -(unsigned) value);
568  } else
569    dst = writeUnsignedDecimal(dst, value);
570
571  return dst;
572}
573
574/* Constructors.  */
575void
576APFloat::initialize(const fltSemantics *ourSemantics)
577{
578  unsigned int count;
579
580  semantics = ourSemantics;
581  count = partCount();
582  if (count > 1)
583    significand.parts = new integerPart[count];
584}
585
586void
587APFloat::freeSignificand()
588{
589  if (needsCleanup())
590    delete [] significand.parts;
591}
592
593void
594APFloat::assign(const APFloat &rhs)
595{
596  assert(semantics == rhs.semantics);
597
598  sign = rhs.sign;
599  category = rhs.category;
600  exponent = rhs.exponent;
601  if (isFiniteNonZero() || category == fcNaN)
602    copySignificand(rhs);
603}
604
605void
606APFloat::copySignificand(const APFloat &rhs)
607{
608  assert(isFiniteNonZero() || category == fcNaN);
609  assert(rhs.partCount() >= partCount());
610
611  APInt::tcAssign(significandParts(), rhs.significandParts(),
612                  partCount());
613}
614
615/* Make this number a NaN, with an arbitrary but deterministic value
616   for the significand.  If double or longer, this is a signalling NaN,
617   which may not be ideal.  If float, this is QNaN(0).  */
618void APFloat::makeNaN(bool SNaN, bool Negative, const APInt *fill)
619{
620  category = fcNaN;
621  sign = Negative;
622
623  integerPart *significand = significandParts();
624  unsigned numParts = partCount();
625
626  // Set the significand bits to the fill.
627  if (!fill || fill->getNumWords() < numParts)
628    APInt::tcSet(significand, 0, numParts);
629  if (fill) {
630    APInt::tcAssign(significand, fill->getRawData(),
631                    std::min(fill->getNumWords(), numParts));
632
633    // Zero out the excess bits of the significand.
634    unsigned bitsToPreserve = semantics->precision - 1;
635    unsigned part = bitsToPreserve / 64;
636    bitsToPreserve %= 64;
637    significand[part] &= ((1ULL << bitsToPreserve) - 1);
638    for (part++; part != numParts; ++part)
639      significand[part] = 0;
640  }
641
642  unsigned QNaNBit = semantics->precision - 2;
643
644  if (SNaN) {
645    // We always have to clear the QNaN bit to make it an SNaN.
646    APInt::tcClearBit(significand, QNaNBit);
647
648    // If there are no bits set in the payload, we have to set
649    // *something* to make it a NaN instead of an infinity;
650    // conventionally, this is the next bit down from the QNaN bit.
651    if (APInt::tcIsZero(significand, numParts))
652      APInt::tcSetBit(significand, QNaNBit - 1);
653  } else {
654    // We always have to set the QNaN bit to make it a QNaN.
655    APInt::tcSetBit(significand, QNaNBit);
656  }
657
658  // For x87 extended precision, we want to make a NaN, not a
659  // pseudo-NaN.  Maybe we should expose the ability to make
660  // pseudo-NaNs?
661  if (semantics == &APFloat::x87DoubleExtended)
662    APInt::tcSetBit(significand, QNaNBit + 1);
663}
664
665APFloat APFloat::makeNaN(const fltSemantics &Sem, bool SNaN, bool Negative,
666                         const APInt *fill) {
667  APFloat value(Sem, uninitialized);
668  value.makeNaN(SNaN, Negative, fill);
669  return value;
670}
671
672APFloat &
673APFloat::operator=(const APFloat &rhs)
674{
675  if (this != &rhs) {
676    if (semantics != rhs.semantics) {
677      freeSignificand();
678      initialize(rhs.semantics);
679    }
680    assign(rhs);
681  }
682
683  return *this;
684}
685
686bool
687APFloat::isDenormal() const {
688  return isFiniteNonZero() && (exponent == semantics->minExponent) &&
689         (APInt::tcExtractBit(significandParts(),
690                              semantics->precision - 1) == 0);
691}
692
693bool
694APFloat::isSmallest() const {
695  // The smallest number by magnitude in our format will be the smallest
696  // denormal, i.e. the floating point number with exponent being minimum
697  // exponent and significand bitwise equal to 1 (i.e. with MSB equal to 0).
698  return isFiniteNonZero() && exponent == semantics->minExponent &&
699    significandMSB() == 0;
700}
701
702bool APFloat::isSignificandAllOnes() const {
703  // Test if the significand excluding the integral bit is all ones. This allows
704  // us to test for binade boundaries.
705  const integerPart *Parts = significandParts();
706  const unsigned PartCount = partCount();
707  for (unsigned i = 0; i < PartCount - 1; i++)
708    if (~Parts[i])
709      return false;
710
711  // Set the unused high bits to all ones when we compare.
712  const unsigned NumHighBits =
713    PartCount*integerPartWidth - semantics->precision + 1;
714  assert(NumHighBits <= integerPartWidth && "Can not have more high bits to "
715         "fill than integerPartWidth");
716  const integerPart HighBitFill =
717    ~integerPart(0) << (integerPartWidth - NumHighBits);
718  if (~(Parts[PartCount - 1] | HighBitFill))
719    return false;
720
721  return true;
722}
723
724bool APFloat::isSignificandAllZeros() const {
725  // Test if the significand excluding the integral bit is all zeros. This
726  // allows us to test for binade boundaries.
727  const integerPart *Parts = significandParts();
728  const unsigned PartCount = partCount();
729
730  for (unsigned i = 0; i < PartCount - 1; i++)
731    if (Parts[i])
732      return false;
733
734  const unsigned NumHighBits =
735    PartCount*integerPartWidth - semantics->precision + 1;
736  assert(NumHighBits <= integerPartWidth && "Can not have more high bits to "
737         "clear than integerPartWidth");
738  const integerPart HighBitMask = ~integerPart(0) >> NumHighBits;
739
740  if (Parts[PartCount - 1] & HighBitMask)
741    return false;
742
743  return true;
744}
745
746bool
747APFloat::isLargest() const {
748  // The largest number by magnitude in our format will be the floating point
749  // number with maximum exponent and with significand that is all ones.
750  return isFiniteNonZero() && exponent == semantics->maxExponent
751    && isSignificandAllOnes();
752}
753
754bool
755APFloat::bitwiseIsEqual(const APFloat &rhs) const {
756  if (this == &rhs)
757    return true;
758  if (semantics != rhs.semantics ||
759      category != rhs.category ||
760      sign != rhs.sign)
761    return false;
762  if (category==fcZero || category==fcInfinity)
763    return true;
764  else if (isFiniteNonZero() && exponent!=rhs.exponent)
765    return false;
766  else {
767    int i= partCount();
768    const integerPart* p=significandParts();
769    const integerPart* q=rhs.significandParts();
770    for (; i>0; i--, p++, q++) {
771      if (*p != *q)
772        return false;
773    }
774    return true;
775  }
776}
777
778APFloat::APFloat(const fltSemantics &ourSemantics, integerPart value) {
779  initialize(&ourSemantics);
780  sign = 0;
781  category = fcNormal;
782  zeroSignificand();
783  exponent = ourSemantics.precision - 1;
784  significandParts()[0] = value;
785  normalize(rmNearestTiesToEven, lfExactlyZero);
786}
787
788APFloat::APFloat(const fltSemantics &ourSemantics) {
789  initialize(&ourSemantics);
790  category = fcZero;
791  sign = false;
792}
793
794APFloat::APFloat(const fltSemantics &ourSemantics, uninitializedTag tag) {
795  // Allocates storage if necessary but does not initialize it.
796  initialize(&ourSemantics);
797}
798
799APFloat::APFloat(const fltSemantics &ourSemantics, StringRef text) {
800  initialize(&ourSemantics);
801  convertFromString(text, rmNearestTiesToEven);
802}
803
804APFloat::APFloat(const APFloat &rhs) {
805  initialize(rhs.semantics);
806  assign(rhs);
807}
808
809APFloat::~APFloat()
810{
811  freeSignificand();
812}
813
814// Profile - This method 'profiles' an APFloat for use with FoldingSet.
815void APFloat::Profile(FoldingSetNodeID& ID) const {
816  ID.Add(bitcastToAPInt());
817}
818
819unsigned int
820APFloat::partCount() const
821{
822  return partCountForBits(semantics->precision + 1);
823}
824
825unsigned int
826APFloat::semanticsPrecision(const fltSemantics &semantics)
827{
828  return semantics.precision;
829}
830
831const integerPart *
832APFloat::significandParts() const
833{
834  return const_cast<APFloat *>(this)->significandParts();
835}
836
837integerPart *
838APFloat::significandParts()
839{
840  if (partCount() > 1)
841    return significand.parts;
842  else
843    return &significand.part;
844}
845
846void
847APFloat::zeroSignificand()
848{
849  APInt::tcSet(significandParts(), 0, partCount());
850}
851
852/* Increment an fcNormal floating point number's significand.  */
853void
854APFloat::incrementSignificand()
855{
856  integerPart carry;
857
858  carry = APInt::tcIncrement(significandParts(), partCount());
859
860  /* Our callers should never cause us to overflow.  */
861  assert(carry == 0);
862  (void)carry;
863}
864
865/* Add the significand of the RHS.  Returns the carry flag.  */
866integerPart
867APFloat::addSignificand(const APFloat &rhs)
868{
869  integerPart *parts;
870
871  parts = significandParts();
872
873  assert(semantics == rhs.semantics);
874  assert(exponent == rhs.exponent);
875
876  return APInt::tcAdd(parts, rhs.significandParts(), 0, partCount());
877}
878
879/* Subtract the significand of the RHS with a borrow flag.  Returns
880   the borrow flag.  */
881integerPart
882APFloat::subtractSignificand(const APFloat &rhs, integerPart borrow)
883{
884  integerPart *parts;
885
886  parts = significandParts();
887
888  assert(semantics == rhs.semantics);
889  assert(exponent == rhs.exponent);
890
891  return APInt::tcSubtract(parts, rhs.significandParts(), borrow,
892                           partCount());
893}
894
895/* Multiply the significand of the RHS.  If ADDEND is non-NULL, add it
896   on to the full-precision result of the multiplication.  Returns the
897   lost fraction.  */
898lostFraction
899APFloat::multiplySignificand(const APFloat &rhs, const APFloat *addend)
900{
901  unsigned int omsb;        // One, not zero, based MSB.
902  unsigned int partsCount, newPartsCount, precision;
903  integerPart *lhsSignificand;
904  integerPart scratch[4];
905  integerPart *fullSignificand;
906  lostFraction lost_fraction;
907  bool ignored;
908
909  assert(semantics == rhs.semantics);
910
911  precision = semantics->precision;
912  newPartsCount = partCountForBits(precision * 2);
913
914  if (newPartsCount > 4)
915    fullSignificand = new integerPart[newPartsCount];
916  else
917    fullSignificand = scratch;
918
919  lhsSignificand = significandParts();
920  partsCount = partCount();
921
922  APInt::tcFullMultiply(fullSignificand, lhsSignificand,
923                        rhs.significandParts(), partsCount, partsCount);
924
925  lost_fraction = lfExactlyZero;
926  omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
927  exponent += rhs.exponent;
928
929  // Assume the operands involved in the multiplication are single-precision
930  // FP, and the two multiplicants are:
931  //   *this = a23 . a22 ... a0 * 2^e1
932  //     rhs = b23 . b22 ... b0 * 2^e2
933  // the result of multiplication is:
934  //   *this = c47 c46 . c45 ... c0 * 2^(e1+e2)
935  // Note that there are two significant bits at the left-hand side of the
936  // radix point. Move the radix point toward left by one bit, and adjust
937  // exponent accordingly.
938  exponent += 1;
939
940  if (addend) {
941    // The intermediate result of the multiplication has "2 * precision"
942    // signicant bit; adjust the addend to be consistent with mul result.
943    //
944    Significand savedSignificand = significand;
945    const fltSemantics *savedSemantics = semantics;
946    fltSemantics extendedSemantics;
947    opStatus status;
948    unsigned int extendedPrecision;
949
950    /* Normalize our MSB.  */
951    extendedPrecision = 2 * precision;
952    if (omsb != extendedPrecision) {
953      assert(extendedPrecision > omsb);
954      APInt::tcShiftLeft(fullSignificand, newPartsCount,
955                         extendedPrecision - omsb);
956      exponent -= extendedPrecision - omsb;
957    }
958
959    /* Create new semantics.  */
960    extendedSemantics = *semantics;
961    extendedSemantics.precision = extendedPrecision;
962
963    if (newPartsCount == 1)
964      significand.part = fullSignificand[0];
965    else
966      significand.parts = fullSignificand;
967    semantics = &extendedSemantics;
968
969    APFloat extendedAddend(*addend);
970    status = extendedAddend.convert(extendedSemantics, rmTowardZero, &ignored);
971    assert(status == opOK);
972    (void)status;
973    lost_fraction = addOrSubtractSignificand(extendedAddend, false);
974
975    /* Restore our state.  */
976    if (newPartsCount == 1)
977      fullSignificand[0] = significand.part;
978    significand = savedSignificand;
979    semantics = savedSemantics;
980
981    omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
982  }
983
984  // Convert the result having "2 * precision" significant-bits back to the one
985  // having "precision" significant-bits. First, move the radix point from
986  // poision "2*precision - 1" to "precision - 1". The exponent need to be
987  // adjusted by "2*precision - 1" - "precision - 1" = "precision".
988  exponent -= precision;
989
990  // In case MSB resides at the left-hand side of radix point, shift the
991  // mantissa right by some amount to make sure the MSB reside right before
992  // the radix point (i.e. "MSB . rest-significant-bits").
993  //
994  // Note that the result is not normalized when "omsb < precision". So, the
995  // caller needs to call APFloat::normalize() if normalized value is expected.
996  if (omsb > precision) {
997    unsigned int bits, significantParts;
998    lostFraction lf;
999
1000    bits = omsb - precision;
1001    significantParts = partCountForBits(omsb);
1002    lf = shiftRight(fullSignificand, significantParts, bits);
1003    lost_fraction = combineLostFractions(lf, lost_fraction);
1004    exponent += bits;
1005  }
1006
1007  APInt::tcAssign(lhsSignificand, fullSignificand, partsCount);
1008
1009  if (newPartsCount > 4)
1010    delete [] fullSignificand;
1011
1012  return lost_fraction;
1013}
1014
1015/* Multiply the significands of LHS and RHS to DST.  */
1016lostFraction
1017APFloat::divideSignificand(const APFloat &rhs)
1018{
1019  unsigned int bit, i, partsCount;
1020  const integerPart *rhsSignificand;
1021  integerPart *lhsSignificand, *dividend, *divisor;
1022  integerPart scratch[4];
1023  lostFraction lost_fraction;
1024
1025  assert(semantics == rhs.semantics);
1026
1027  lhsSignificand = significandParts();
1028  rhsSignificand = rhs.significandParts();
1029  partsCount = partCount();
1030
1031  if (partsCount > 2)
1032    dividend = new integerPart[partsCount * 2];
1033  else
1034    dividend = scratch;
1035
1036  divisor = dividend + partsCount;
1037
1038  /* Copy the dividend and divisor as they will be modified in-place.  */
1039  for (i = 0; i < partsCount; i++) {
1040    dividend[i] = lhsSignificand[i];
1041    divisor[i] = rhsSignificand[i];
1042    lhsSignificand[i] = 0;
1043  }
1044
1045  exponent -= rhs.exponent;
1046
1047  unsigned int precision = semantics->precision;
1048
1049  /* Normalize the divisor.  */
1050  bit = precision - APInt::tcMSB(divisor, partsCount) - 1;
1051  if (bit) {
1052    exponent += bit;
1053    APInt::tcShiftLeft(divisor, partsCount, bit);
1054  }
1055
1056  /* Normalize the dividend.  */
1057  bit = precision - APInt::tcMSB(dividend, partsCount) - 1;
1058  if (bit) {
1059    exponent -= bit;
1060    APInt::tcShiftLeft(dividend, partsCount, bit);
1061  }
1062
1063  /* Ensure the dividend >= divisor initially for the loop below.
1064     Incidentally, this means that the division loop below is
1065     guaranteed to set the integer bit to one.  */
1066  if (APInt::tcCompare(dividend, divisor, partsCount) < 0) {
1067    exponent--;
1068    APInt::tcShiftLeft(dividend, partsCount, 1);
1069    assert(APInt::tcCompare(dividend, divisor, partsCount) >= 0);
1070  }
1071
1072  /* Long division.  */
1073  for (bit = precision; bit; bit -= 1) {
1074    if (APInt::tcCompare(dividend, divisor, partsCount) >= 0) {
1075      APInt::tcSubtract(dividend, divisor, 0, partsCount);
1076      APInt::tcSetBit(lhsSignificand, bit - 1);
1077    }
1078
1079    APInt::tcShiftLeft(dividend, partsCount, 1);
1080  }
1081
1082  /* Figure out the lost fraction.  */
1083  int cmp = APInt::tcCompare(dividend, divisor, partsCount);
1084
1085  if (cmp > 0)
1086    lost_fraction = lfMoreThanHalf;
1087  else if (cmp == 0)
1088    lost_fraction = lfExactlyHalf;
1089  else if (APInt::tcIsZero(dividend, partsCount))
1090    lost_fraction = lfExactlyZero;
1091  else
1092    lost_fraction = lfLessThanHalf;
1093
1094  if (partsCount > 2)
1095    delete [] dividend;
1096
1097  return lost_fraction;
1098}
1099
1100unsigned int
1101APFloat::significandMSB() const
1102{
1103  return APInt::tcMSB(significandParts(), partCount());
1104}
1105
1106unsigned int
1107APFloat::significandLSB() const
1108{
1109  return APInt::tcLSB(significandParts(), partCount());
1110}
1111
1112/* Note that a zero result is NOT normalized to fcZero.  */
1113lostFraction
1114APFloat::shiftSignificandRight(unsigned int bits)
1115{
1116  /* Our exponent should not overflow.  */
1117  assert((ExponentType) (exponent + bits) >= exponent);
1118
1119  exponent += bits;
1120
1121  return shiftRight(significandParts(), partCount(), bits);
1122}
1123
1124/* Shift the significand left BITS bits, subtract BITS from its exponent.  */
1125void
1126APFloat::shiftSignificandLeft(unsigned int bits)
1127{
1128  assert(bits < semantics->precision);
1129
1130  if (bits) {
1131    unsigned int partsCount = partCount();
1132
1133    APInt::tcShiftLeft(significandParts(), partsCount, bits);
1134    exponent -= bits;
1135
1136    assert(!APInt::tcIsZero(significandParts(), partsCount));
1137  }
1138}
1139
1140APFloat::cmpResult
1141APFloat::compareAbsoluteValue(const APFloat &rhs) const
1142{
1143  int compare;
1144
1145  assert(semantics == rhs.semantics);
1146  assert(isFiniteNonZero());
1147  assert(rhs.isFiniteNonZero());
1148
1149  compare = exponent - rhs.exponent;
1150
1151  /* If exponents are equal, do an unsigned bignum comparison of the
1152     significands.  */
1153  if (compare == 0)
1154    compare = APInt::tcCompare(significandParts(), rhs.significandParts(),
1155                               partCount());
1156
1157  if (compare > 0)
1158    return cmpGreaterThan;
1159  else if (compare < 0)
1160    return cmpLessThan;
1161  else
1162    return cmpEqual;
1163}
1164
1165/* Handle overflow.  Sign is preserved.  We either become infinity or
1166   the largest finite number.  */
1167APFloat::opStatus
1168APFloat::handleOverflow(roundingMode rounding_mode)
1169{
1170  /* Infinity?  */
1171  if (rounding_mode == rmNearestTiesToEven ||
1172      rounding_mode == rmNearestTiesToAway ||
1173      (rounding_mode == rmTowardPositive && !sign) ||
1174      (rounding_mode == rmTowardNegative && sign)) {
1175    category = fcInfinity;
1176    return (opStatus) (opOverflow | opInexact);
1177  }
1178
1179  /* Otherwise we become the largest finite number.  */
1180  category = fcNormal;
1181  exponent = semantics->maxExponent;
1182  APInt::tcSetLeastSignificantBits(significandParts(), partCount(),
1183                                   semantics->precision);
1184
1185  return opInexact;
1186}
1187
1188/* Returns TRUE if, when truncating the current number, with BIT the
1189   new LSB, with the given lost fraction and rounding mode, the result
1190   would need to be rounded away from zero (i.e., by increasing the
1191   signficand).  This routine must work for fcZero of both signs, and
1192   fcNormal numbers.  */
1193bool
1194APFloat::roundAwayFromZero(roundingMode rounding_mode,
1195                           lostFraction lost_fraction,
1196                           unsigned int bit) const
1197{
1198  /* NaNs and infinities should not have lost fractions.  */
1199  assert(isFiniteNonZero() || category == fcZero);
1200
1201  /* Current callers never pass this so we don't handle it.  */
1202  assert(lost_fraction != lfExactlyZero);
1203
1204  switch (rounding_mode) {
1205  case rmNearestTiesToAway:
1206    return lost_fraction == lfExactlyHalf || lost_fraction == lfMoreThanHalf;
1207
1208  case rmNearestTiesToEven:
1209    if (lost_fraction == lfMoreThanHalf)
1210      return true;
1211
1212    /* Our zeroes don't have a significand to test.  */
1213    if (lost_fraction == lfExactlyHalf && category != fcZero)
1214      return APInt::tcExtractBit(significandParts(), bit);
1215
1216    return false;
1217
1218  case rmTowardZero:
1219    return false;
1220
1221  case rmTowardPositive:
1222    return sign == false;
1223
1224  case rmTowardNegative:
1225    return sign == true;
1226  }
1227  llvm_unreachable("Invalid rounding mode found");
1228}
1229
1230APFloat::opStatus
1231APFloat::normalize(roundingMode rounding_mode,
1232                   lostFraction lost_fraction)
1233{
1234  unsigned int omsb;                /* One, not zero, based MSB.  */
1235  int exponentChange;
1236
1237  if (!isFiniteNonZero())
1238    return opOK;
1239
1240  /* Before rounding normalize the exponent of fcNormal numbers.  */
1241  omsb = significandMSB() + 1;
1242
1243  if (omsb) {
1244    /* OMSB is numbered from 1.  We want to place it in the integer
1245       bit numbered PRECISION if possible, with a compensating change in
1246       the exponent.  */
1247    exponentChange = omsb - semantics->precision;
1248
1249    /* If the resulting exponent is too high, overflow according to
1250       the rounding mode.  */
1251    if (exponent + exponentChange > semantics->maxExponent)
1252      return handleOverflow(rounding_mode);
1253
1254    /* Subnormal numbers have exponent minExponent, and their MSB
1255       is forced based on that.  */
1256    if (exponent + exponentChange < semantics->minExponent)
1257      exponentChange = semantics->minExponent - exponent;
1258
1259    /* Shifting left is easy as we don't lose precision.  */
1260    if (exponentChange < 0) {
1261      assert(lost_fraction == lfExactlyZero);
1262
1263      shiftSignificandLeft(-exponentChange);
1264
1265      return opOK;
1266    }
1267
1268    if (exponentChange > 0) {
1269      lostFraction lf;
1270
1271      /* Shift right and capture any new lost fraction.  */
1272      lf = shiftSignificandRight(exponentChange);
1273
1274      lost_fraction = combineLostFractions(lf, lost_fraction);
1275
1276      /* Keep OMSB up-to-date.  */
1277      if (omsb > (unsigned) exponentChange)
1278        omsb -= exponentChange;
1279      else
1280        omsb = 0;
1281    }
1282  }
1283
1284  /* Now round the number according to rounding_mode given the lost
1285     fraction.  */
1286
1287  /* As specified in IEEE 754, since we do not trap we do not report
1288     underflow for exact results.  */
1289  if (lost_fraction == lfExactlyZero) {
1290    /* Canonicalize zeroes.  */
1291    if (omsb == 0)
1292      category = fcZero;
1293
1294    return opOK;
1295  }
1296
1297  /* Increment the significand if we're rounding away from zero.  */
1298  if (roundAwayFromZero(rounding_mode, lost_fraction, 0)) {
1299    if (omsb == 0)
1300      exponent = semantics->minExponent;
1301
1302    incrementSignificand();
1303    omsb = significandMSB() + 1;
1304
1305    /* Did the significand increment overflow?  */
1306    if (omsb == (unsigned) semantics->precision + 1) {
1307      /* Renormalize by incrementing the exponent and shifting our
1308         significand right one.  However if we already have the
1309         maximum exponent we overflow to infinity.  */
1310      if (exponent == semantics->maxExponent) {
1311        category = fcInfinity;
1312
1313        return (opStatus) (opOverflow | opInexact);
1314      }
1315
1316      shiftSignificandRight(1);
1317
1318      return opInexact;
1319    }
1320  }
1321
1322  /* The normal case - we were and are not denormal, and any
1323     significand increment above didn't overflow.  */
1324  if (omsb == semantics->precision)
1325    return opInexact;
1326
1327  /* We have a non-zero denormal.  */
1328  assert(omsb < semantics->precision);
1329
1330  /* Canonicalize zeroes.  */
1331  if (omsb == 0)
1332    category = fcZero;
1333
1334  /* The fcZero case is a denormal that underflowed to zero.  */
1335  return (opStatus) (opUnderflow | opInexact);
1336}
1337
1338APFloat::opStatus
1339APFloat::addOrSubtractSpecials(const APFloat &rhs, bool subtract)
1340{
1341  switch (PackCategoriesIntoKey(category, rhs.category)) {
1342  default:
1343    llvm_unreachable(0);
1344
1345  case PackCategoriesIntoKey(fcNaN, fcZero):
1346  case PackCategoriesIntoKey(fcNaN, fcNormal):
1347  case PackCategoriesIntoKey(fcNaN, fcInfinity):
1348  case PackCategoriesIntoKey(fcNaN, fcNaN):
1349  case PackCategoriesIntoKey(fcNormal, fcZero):
1350  case PackCategoriesIntoKey(fcInfinity, fcNormal):
1351  case PackCategoriesIntoKey(fcInfinity, fcZero):
1352    return opOK;
1353
1354  case PackCategoriesIntoKey(fcZero, fcNaN):
1355  case PackCategoriesIntoKey(fcNormal, fcNaN):
1356  case PackCategoriesIntoKey(fcInfinity, fcNaN):
1357    sign = false;
1358    category = fcNaN;
1359    copySignificand(rhs);
1360    return opOK;
1361
1362  case PackCategoriesIntoKey(fcNormal, fcInfinity):
1363  case PackCategoriesIntoKey(fcZero, fcInfinity):
1364    category = fcInfinity;
1365    sign = rhs.sign ^ subtract;
1366    return opOK;
1367
1368  case PackCategoriesIntoKey(fcZero, fcNormal):
1369    assign(rhs);
1370    sign = rhs.sign ^ subtract;
1371    return opOK;
1372
1373  case PackCategoriesIntoKey(fcZero, fcZero):
1374    /* Sign depends on rounding mode; handled by caller.  */
1375    return opOK;
1376
1377  case PackCategoriesIntoKey(fcInfinity, fcInfinity):
1378    /* Differently signed infinities can only be validly
1379       subtracted.  */
1380    if (((sign ^ rhs.sign)!=0) != subtract) {
1381      makeNaN();
1382      return opInvalidOp;
1383    }
1384
1385    return opOK;
1386
1387  case PackCategoriesIntoKey(fcNormal, fcNormal):
1388    return opDivByZero;
1389  }
1390}
1391
1392/* Add or subtract two normal numbers.  */
1393lostFraction
1394APFloat::addOrSubtractSignificand(const APFloat &rhs, bool subtract)
1395{
1396  integerPart carry;
1397  lostFraction lost_fraction;
1398  int bits;
1399
1400  /* Determine if the operation on the absolute values is effectively
1401     an addition or subtraction.  */
1402  subtract ^= (sign ^ rhs.sign) ? true : false;
1403
1404  /* Are we bigger exponent-wise than the RHS?  */
1405  bits = exponent - rhs.exponent;
1406
1407  /* Subtraction is more subtle than one might naively expect.  */
1408  if (subtract) {
1409    APFloat temp_rhs(rhs);
1410    bool reverse;
1411
1412    if (bits == 0) {
1413      reverse = compareAbsoluteValue(temp_rhs) == cmpLessThan;
1414      lost_fraction = lfExactlyZero;
1415    } else if (bits > 0) {
1416      lost_fraction = temp_rhs.shiftSignificandRight(bits - 1);
1417      shiftSignificandLeft(1);
1418      reverse = false;
1419    } else {
1420      lost_fraction = shiftSignificandRight(-bits - 1);
1421      temp_rhs.shiftSignificandLeft(1);
1422      reverse = true;
1423    }
1424
1425    if (reverse) {
1426      carry = temp_rhs.subtractSignificand
1427        (*this, lost_fraction != lfExactlyZero);
1428      copySignificand(temp_rhs);
1429      sign = !sign;
1430    } else {
1431      carry = subtractSignificand
1432        (temp_rhs, lost_fraction != lfExactlyZero);
1433    }
1434
1435    /* Invert the lost fraction - it was on the RHS and
1436       subtracted.  */
1437    if (lost_fraction == lfLessThanHalf)
1438      lost_fraction = lfMoreThanHalf;
1439    else if (lost_fraction == lfMoreThanHalf)
1440      lost_fraction = lfLessThanHalf;
1441
1442    /* The code above is intended to ensure that no borrow is
1443       necessary.  */
1444    assert(!carry);
1445    (void)carry;
1446  } else {
1447    if (bits > 0) {
1448      APFloat temp_rhs(rhs);
1449
1450      lost_fraction = temp_rhs.shiftSignificandRight(bits);
1451      carry = addSignificand(temp_rhs);
1452    } else {
1453      lost_fraction = shiftSignificandRight(-bits);
1454      carry = addSignificand(rhs);
1455    }
1456
1457    /* We have a guard bit; generating a carry cannot happen.  */
1458    assert(!carry);
1459    (void)carry;
1460  }
1461
1462  return lost_fraction;
1463}
1464
1465APFloat::opStatus
1466APFloat::multiplySpecials(const APFloat &rhs)
1467{
1468  switch (PackCategoriesIntoKey(category, rhs.category)) {
1469  default:
1470    llvm_unreachable(0);
1471
1472  case PackCategoriesIntoKey(fcNaN, fcZero):
1473  case PackCategoriesIntoKey(fcNaN, fcNormal):
1474  case PackCategoriesIntoKey(fcNaN, fcInfinity):
1475  case PackCategoriesIntoKey(fcNaN, fcNaN):
1476    sign = false;
1477    return opOK;
1478
1479  case PackCategoriesIntoKey(fcZero, fcNaN):
1480  case PackCategoriesIntoKey(fcNormal, fcNaN):
1481  case PackCategoriesIntoKey(fcInfinity, fcNaN):
1482    sign = false;
1483    category = fcNaN;
1484    copySignificand(rhs);
1485    return opOK;
1486
1487  case PackCategoriesIntoKey(fcNormal, fcInfinity):
1488  case PackCategoriesIntoKey(fcInfinity, fcNormal):
1489  case PackCategoriesIntoKey(fcInfinity, fcInfinity):
1490    category = fcInfinity;
1491    return opOK;
1492
1493  case PackCategoriesIntoKey(fcZero, fcNormal):
1494  case PackCategoriesIntoKey(fcNormal, fcZero):
1495  case PackCategoriesIntoKey(fcZero, fcZero):
1496    category = fcZero;
1497    return opOK;
1498
1499  case PackCategoriesIntoKey(fcZero, fcInfinity):
1500  case PackCategoriesIntoKey(fcInfinity, fcZero):
1501    makeNaN();
1502    return opInvalidOp;
1503
1504  case PackCategoriesIntoKey(fcNormal, fcNormal):
1505    return opOK;
1506  }
1507}
1508
1509APFloat::opStatus
1510APFloat::divideSpecials(const APFloat &rhs)
1511{
1512  switch (PackCategoriesIntoKey(category, rhs.category)) {
1513  default:
1514    llvm_unreachable(0);
1515
1516  case PackCategoriesIntoKey(fcZero, fcNaN):
1517  case PackCategoriesIntoKey(fcNormal, fcNaN):
1518  case PackCategoriesIntoKey(fcInfinity, fcNaN):
1519    category = fcNaN;
1520    copySignificand(rhs);
1521  case PackCategoriesIntoKey(fcNaN, fcZero):
1522  case PackCategoriesIntoKey(fcNaN, fcNormal):
1523  case PackCategoriesIntoKey(fcNaN, fcInfinity):
1524  case PackCategoriesIntoKey(fcNaN, fcNaN):
1525    sign = false;
1526  case PackCategoriesIntoKey(fcInfinity, fcZero):
1527  case PackCategoriesIntoKey(fcInfinity, fcNormal):
1528  case PackCategoriesIntoKey(fcZero, fcInfinity):
1529  case PackCategoriesIntoKey(fcZero, fcNormal):
1530    return opOK;
1531
1532  case PackCategoriesIntoKey(fcNormal, fcInfinity):
1533    category = fcZero;
1534    return opOK;
1535
1536  case PackCategoriesIntoKey(fcNormal, fcZero):
1537    category = fcInfinity;
1538    return opDivByZero;
1539
1540  case PackCategoriesIntoKey(fcInfinity, fcInfinity):
1541  case PackCategoriesIntoKey(fcZero, fcZero):
1542    makeNaN();
1543    return opInvalidOp;
1544
1545  case PackCategoriesIntoKey(fcNormal, fcNormal):
1546    return opOK;
1547  }
1548}
1549
1550APFloat::opStatus
1551APFloat::modSpecials(const APFloat &rhs)
1552{
1553  switch (PackCategoriesIntoKey(category, rhs.category)) {
1554  default:
1555    llvm_unreachable(0);
1556
1557  case PackCategoriesIntoKey(fcNaN, fcZero):
1558  case PackCategoriesIntoKey(fcNaN, fcNormal):
1559  case PackCategoriesIntoKey(fcNaN, fcInfinity):
1560  case PackCategoriesIntoKey(fcNaN, fcNaN):
1561  case PackCategoriesIntoKey(fcZero, fcInfinity):
1562  case PackCategoriesIntoKey(fcZero, fcNormal):
1563  case PackCategoriesIntoKey(fcNormal, fcInfinity):
1564    return opOK;
1565
1566  case PackCategoriesIntoKey(fcZero, fcNaN):
1567  case PackCategoriesIntoKey(fcNormal, fcNaN):
1568  case PackCategoriesIntoKey(fcInfinity, fcNaN):
1569    sign = false;
1570    category = fcNaN;
1571    copySignificand(rhs);
1572    return opOK;
1573
1574  case PackCategoriesIntoKey(fcNormal, fcZero):
1575  case PackCategoriesIntoKey(fcInfinity, fcZero):
1576  case PackCategoriesIntoKey(fcInfinity, fcNormal):
1577  case PackCategoriesIntoKey(fcInfinity, fcInfinity):
1578  case PackCategoriesIntoKey(fcZero, fcZero):
1579    makeNaN();
1580    return opInvalidOp;
1581
1582  case PackCategoriesIntoKey(fcNormal, fcNormal):
1583    return opOK;
1584  }
1585}
1586
1587/* Change sign.  */
1588void
1589APFloat::changeSign()
1590{
1591  /* Look mummy, this one's easy.  */
1592  sign = !sign;
1593}
1594
1595void
1596APFloat::clearSign()
1597{
1598  /* So is this one. */
1599  sign = 0;
1600}
1601
1602void
1603APFloat::copySign(const APFloat &rhs)
1604{
1605  /* And this one. */
1606  sign = rhs.sign;
1607}
1608
1609/* Normalized addition or subtraction.  */
1610APFloat::opStatus
1611APFloat::addOrSubtract(const APFloat &rhs, roundingMode rounding_mode,
1612                       bool subtract)
1613{
1614  opStatus fs;
1615
1616  fs = addOrSubtractSpecials(rhs, subtract);
1617
1618  /* This return code means it was not a simple case.  */
1619  if (fs == opDivByZero) {
1620    lostFraction lost_fraction;
1621
1622    lost_fraction = addOrSubtractSignificand(rhs, subtract);
1623    fs = normalize(rounding_mode, lost_fraction);
1624
1625    /* Can only be zero if we lost no fraction.  */
1626    assert(category != fcZero || lost_fraction == lfExactlyZero);
1627  }
1628
1629  /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1630     positive zero unless rounding to minus infinity, except that
1631     adding two like-signed zeroes gives that zero.  */
1632  if (category == fcZero) {
1633    if (rhs.category != fcZero || (sign == rhs.sign) == subtract)
1634      sign = (rounding_mode == rmTowardNegative);
1635  }
1636
1637  return fs;
1638}
1639
1640/* Normalized addition.  */
1641APFloat::opStatus
1642APFloat::add(const APFloat &rhs, roundingMode rounding_mode)
1643{
1644  return addOrSubtract(rhs, rounding_mode, false);
1645}
1646
1647/* Normalized subtraction.  */
1648APFloat::opStatus
1649APFloat::subtract(const APFloat &rhs, roundingMode rounding_mode)
1650{
1651  return addOrSubtract(rhs, rounding_mode, true);
1652}
1653
1654/* Normalized multiply.  */
1655APFloat::opStatus
1656APFloat::multiply(const APFloat &rhs, roundingMode rounding_mode)
1657{
1658  opStatus fs;
1659
1660  sign ^= rhs.sign;
1661  fs = multiplySpecials(rhs);
1662
1663  if (isFiniteNonZero()) {
1664    lostFraction lost_fraction = multiplySignificand(rhs, 0);
1665    fs = normalize(rounding_mode, lost_fraction);
1666    if (lost_fraction != lfExactlyZero)
1667      fs = (opStatus) (fs | opInexact);
1668  }
1669
1670  return fs;
1671}
1672
1673/* Normalized divide.  */
1674APFloat::opStatus
1675APFloat::divide(const APFloat &rhs, roundingMode rounding_mode)
1676{
1677  opStatus fs;
1678
1679  sign ^= rhs.sign;
1680  fs = divideSpecials(rhs);
1681
1682  if (isFiniteNonZero()) {
1683    lostFraction lost_fraction = divideSignificand(rhs);
1684    fs = normalize(rounding_mode, lost_fraction);
1685    if (lost_fraction != lfExactlyZero)
1686      fs = (opStatus) (fs | opInexact);
1687  }
1688
1689  return fs;
1690}
1691
1692/* Normalized remainder.  This is not currently correct in all cases.  */
1693APFloat::opStatus
1694APFloat::remainder(const APFloat &rhs)
1695{
1696  opStatus fs;
1697  APFloat V = *this;
1698  unsigned int origSign = sign;
1699
1700  fs = V.divide(rhs, rmNearestTiesToEven);
1701  if (fs == opDivByZero)
1702    return fs;
1703
1704  int parts = partCount();
1705  integerPart *x = new integerPart[parts];
1706  bool ignored;
1707  fs = V.convertToInteger(x, parts * integerPartWidth, true,
1708                          rmNearestTiesToEven, &ignored);
1709  if (fs==opInvalidOp)
1710    return fs;
1711
1712  fs = V.convertFromZeroExtendedInteger(x, parts * integerPartWidth, true,
1713                                        rmNearestTiesToEven);
1714  assert(fs==opOK);   // should always work
1715
1716  fs = V.multiply(rhs, rmNearestTiesToEven);
1717  assert(fs==opOK || fs==opInexact);   // should not overflow or underflow
1718
1719  fs = subtract(V, rmNearestTiesToEven);
1720  assert(fs==opOK || fs==opInexact);   // likewise
1721
1722  if (isZero())
1723    sign = origSign;    // IEEE754 requires this
1724  delete[] x;
1725  return fs;
1726}
1727
1728/* Normalized llvm frem (C fmod).
1729   This is not currently correct in all cases.  */
1730APFloat::opStatus
1731APFloat::mod(const APFloat &rhs, roundingMode rounding_mode)
1732{
1733  opStatus fs;
1734  fs = modSpecials(rhs);
1735
1736  if (isFiniteNonZero() && rhs.isFiniteNonZero()) {
1737    APFloat V = *this;
1738    unsigned int origSign = sign;
1739
1740    fs = V.divide(rhs, rmNearestTiesToEven);
1741    if (fs == opDivByZero)
1742      return fs;
1743
1744    int parts = partCount();
1745    integerPart *x = new integerPart[parts];
1746    bool ignored;
1747    fs = V.convertToInteger(x, parts * integerPartWidth, true,
1748                            rmTowardZero, &ignored);
1749    if (fs==opInvalidOp)
1750      return fs;
1751
1752    fs = V.convertFromZeroExtendedInteger(x, parts * integerPartWidth, true,
1753                                          rmNearestTiesToEven);
1754    assert(fs==opOK);   // should always work
1755
1756    fs = V.multiply(rhs, rounding_mode);
1757    assert(fs==opOK || fs==opInexact);   // should not overflow or underflow
1758
1759    fs = subtract(V, rounding_mode);
1760    assert(fs==opOK || fs==opInexact);   // likewise
1761
1762    if (isZero())
1763      sign = origSign;    // IEEE754 requires this
1764    delete[] x;
1765  }
1766  return fs;
1767}
1768
1769/* Normalized fused-multiply-add.  */
1770APFloat::opStatus
1771APFloat::fusedMultiplyAdd(const APFloat &multiplicand,
1772                          const APFloat &addend,
1773                          roundingMode rounding_mode)
1774{
1775  opStatus fs;
1776
1777  /* Post-multiplication sign, before addition.  */
1778  sign ^= multiplicand.sign;
1779
1780  /* If and only if all arguments are normal do we need to do an
1781     extended-precision calculation.  */
1782  if (isFiniteNonZero() &&
1783      multiplicand.isFiniteNonZero() &&
1784      addend.isFiniteNonZero()) {
1785    lostFraction lost_fraction;
1786
1787    lost_fraction = multiplySignificand(multiplicand, &addend);
1788    fs = normalize(rounding_mode, lost_fraction);
1789    if (lost_fraction != lfExactlyZero)
1790      fs = (opStatus) (fs | opInexact);
1791
1792    /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1793       positive zero unless rounding to minus infinity, except that
1794       adding two like-signed zeroes gives that zero.  */
1795    if (category == fcZero && sign != addend.sign)
1796      sign = (rounding_mode == rmTowardNegative);
1797  } else {
1798    fs = multiplySpecials(multiplicand);
1799
1800    /* FS can only be opOK or opInvalidOp.  There is no more work
1801       to do in the latter case.  The IEEE-754R standard says it is
1802       implementation-defined in this case whether, if ADDEND is a
1803       quiet NaN, we raise invalid op; this implementation does so.
1804
1805       If we need to do the addition we can do so with normal
1806       precision.  */
1807    if (fs == opOK)
1808      fs = addOrSubtract(addend, rounding_mode, false);
1809  }
1810
1811  return fs;
1812}
1813
1814/* Rounding-mode corrrect round to integral value.  */
1815APFloat::opStatus APFloat::roundToIntegral(roundingMode rounding_mode) {
1816  opStatus fs;
1817
1818  // If the exponent is large enough, we know that this value is already
1819  // integral, and the arithmetic below would potentially cause it to saturate
1820  // to +/-Inf.  Bail out early instead.
1821  if (isFiniteNonZero() && exponent+1 >= (int)semanticsPrecision(*semantics))
1822    return opOK;
1823
1824  // The algorithm here is quite simple: we add 2^(p-1), where p is the
1825  // precision of our format, and then subtract it back off again.  The choice
1826  // of rounding modes for the addition/subtraction determines the rounding mode
1827  // for our integral rounding as well.
1828  // NOTE: When the input value is negative, we do subtraction followed by
1829  // addition instead.
1830  APInt IntegerConstant(NextPowerOf2(semanticsPrecision(*semantics)), 1);
1831  IntegerConstant <<= semanticsPrecision(*semantics)-1;
1832  APFloat MagicConstant(*semantics);
1833  fs = MagicConstant.convertFromAPInt(IntegerConstant, false,
1834                                      rmNearestTiesToEven);
1835  MagicConstant.copySign(*this);
1836
1837  if (fs != opOK)
1838    return fs;
1839
1840  // Preserve the input sign so that we can handle 0.0/-0.0 cases correctly.
1841  bool inputSign = isNegative();
1842
1843  fs = add(MagicConstant, rounding_mode);
1844  if (fs != opOK && fs != opInexact)
1845    return fs;
1846
1847  fs = subtract(MagicConstant, rounding_mode);
1848
1849  // Restore the input sign.
1850  if (inputSign != isNegative())
1851    changeSign();
1852
1853  return fs;
1854}
1855
1856
1857/* Comparison requires normalized numbers.  */
1858APFloat::cmpResult
1859APFloat::compare(const APFloat &rhs) const
1860{
1861  cmpResult result;
1862
1863  assert(semantics == rhs.semantics);
1864
1865  switch (PackCategoriesIntoKey(category, rhs.category)) {
1866  default:
1867    llvm_unreachable(0);
1868
1869  case PackCategoriesIntoKey(fcNaN, fcZero):
1870  case PackCategoriesIntoKey(fcNaN, fcNormal):
1871  case PackCategoriesIntoKey(fcNaN, fcInfinity):
1872  case PackCategoriesIntoKey(fcNaN, fcNaN):
1873  case PackCategoriesIntoKey(fcZero, fcNaN):
1874  case PackCategoriesIntoKey(fcNormal, fcNaN):
1875  case PackCategoriesIntoKey(fcInfinity, fcNaN):
1876    return cmpUnordered;
1877
1878  case PackCategoriesIntoKey(fcInfinity, fcNormal):
1879  case PackCategoriesIntoKey(fcInfinity, fcZero):
1880  case PackCategoriesIntoKey(fcNormal, fcZero):
1881    if (sign)
1882      return cmpLessThan;
1883    else
1884      return cmpGreaterThan;
1885
1886  case PackCategoriesIntoKey(fcNormal, fcInfinity):
1887  case PackCategoriesIntoKey(fcZero, fcInfinity):
1888  case PackCategoriesIntoKey(fcZero, fcNormal):
1889    if (rhs.sign)
1890      return cmpGreaterThan;
1891    else
1892      return cmpLessThan;
1893
1894  case PackCategoriesIntoKey(fcInfinity, fcInfinity):
1895    if (sign == rhs.sign)
1896      return cmpEqual;
1897    else if (sign)
1898      return cmpLessThan;
1899    else
1900      return cmpGreaterThan;
1901
1902  case PackCategoriesIntoKey(fcZero, fcZero):
1903    return cmpEqual;
1904
1905  case PackCategoriesIntoKey(fcNormal, fcNormal):
1906    break;
1907  }
1908
1909  /* Two normal numbers.  Do they have the same sign?  */
1910  if (sign != rhs.sign) {
1911    if (sign)
1912      result = cmpLessThan;
1913    else
1914      result = cmpGreaterThan;
1915  } else {
1916    /* Compare absolute values; invert result if negative.  */
1917    result = compareAbsoluteValue(rhs);
1918
1919    if (sign) {
1920      if (result == cmpLessThan)
1921        result = cmpGreaterThan;
1922      else if (result == cmpGreaterThan)
1923        result = cmpLessThan;
1924    }
1925  }
1926
1927  return result;
1928}
1929
1930/// APFloat::convert - convert a value of one floating point type to another.
1931/// The return value corresponds to the IEEE754 exceptions.  *losesInfo
1932/// records whether the transformation lost information, i.e. whether
1933/// converting the result back to the original type will produce the
1934/// original value (this is almost the same as return value==fsOK, but there
1935/// are edge cases where this is not so).
1936
1937APFloat::opStatus
1938APFloat::convert(const fltSemantics &toSemantics,
1939                 roundingMode rounding_mode, bool *losesInfo)
1940{
1941  lostFraction lostFraction;
1942  unsigned int newPartCount, oldPartCount;
1943  opStatus fs;
1944  int shift;
1945  const fltSemantics &fromSemantics = *semantics;
1946
1947  lostFraction = lfExactlyZero;
1948  newPartCount = partCountForBits(toSemantics.precision + 1);
1949  oldPartCount = partCount();
1950  shift = toSemantics.precision - fromSemantics.precision;
1951
1952  bool X86SpecialNan = false;
1953  if (&fromSemantics == &APFloat::x87DoubleExtended &&
1954      &toSemantics != &APFloat::x87DoubleExtended && category == fcNaN &&
1955      (!(*significandParts() & 0x8000000000000000ULL) ||
1956       !(*significandParts() & 0x4000000000000000ULL))) {
1957    // x86 has some unusual NaNs which cannot be represented in any other
1958    // format; note them here.
1959    X86SpecialNan = true;
1960  }
1961
1962  // If this is a truncation of a denormal number, and the target semantics
1963  // has larger exponent range than the source semantics (this can happen
1964  // when truncating from PowerPC double-double to double format), the
1965  // right shift could lose result mantissa bits.  Adjust exponent instead
1966  // of performing excessive shift.
1967  if (shift < 0 && isFiniteNonZero()) {
1968    int exponentChange = significandMSB() + 1 - fromSemantics.precision;
1969    if (exponent + exponentChange < toSemantics.minExponent)
1970      exponentChange = toSemantics.minExponent - exponent;
1971    if (exponentChange < shift)
1972      exponentChange = shift;
1973    if (exponentChange < 0) {
1974      shift -= exponentChange;
1975      exponent += exponentChange;
1976    }
1977  }
1978
1979  // If this is a truncation, perform the shift before we narrow the storage.
1980  if (shift < 0 && (isFiniteNonZero() || category==fcNaN))
1981    lostFraction = shiftRight(significandParts(), oldPartCount, -shift);
1982
1983  // Fix the storage so it can hold to new value.
1984  if (newPartCount > oldPartCount) {
1985    // The new type requires more storage; make it available.
1986    integerPart *newParts;
1987    newParts = new integerPart[newPartCount];
1988    APInt::tcSet(newParts, 0, newPartCount);
1989    if (isFiniteNonZero() || category==fcNaN)
1990      APInt::tcAssign(newParts, significandParts(), oldPartCount);
1991    freeSignificand();
1992    significand.parts = newParts;
1993  } else if (newPartCount == 1 && oldPartCount != 1) {
1994    // Switch to built-in storage for a single part.
1995    integerPart newPart = 0;
1996    if (isFiniteNonZero() || category==fcNaN)
1997      newPart = significandParts()[0];
1998    freeSignificand();
1999    significand.part = newPart;
2000  }
2001
2002  // Now that we have the right storage, switch the semantics.
2003  semantics = &toSemantics;
2004
2005  // If this is an extension, perform the shift now that the storage is
2006  // available.
2007  if (shift > 0 && (isFiniteNonZero() || category==fcNaN))
2008    APInt::tcShiftLeft(significandParts(), newPartCount, shift);
2009
2010  if (isFiniteNonZero()) {
2011    fs = normalize(rounding_mode, lostFraction);
2012    *losesInfo = (fs != opOK);
2013  } else if (category == fcNaN) {
2014    *losesInfo = lostFraction != lfExactlyZero || X86SpecialNan;
2015
2016    // For x87 extended precision, we want to make a NaN, not a special NaN if
2017    // the input wasn't special either.
2018    if (!X86SpecialNan && semantics == &APFloat::x87DoubleExtended)
2019      APInt::tcSetBit(significandParts(), semantics->precision - 1);
2020
2021    // gcc forces the Quiet bit on, which means (float)(double)(float_sNan)
2022    // does not give you back the same bits.  This is dubious, and we
2023    // don't currently do it.  You're really supposed to get
2024    // an invalid operation signal at runtime, but nobody does that.
2025    fs = opOK;
2026  } else {
2027    *losesInfo = false;
2028    fs = opOK;
2029  }
2030
2031  return fs;
2032}
2033
2034/* Convert a floating point number to an integer according to the
2035   rounding mode.  If the rounded integer value is out of range this
2036   returns an invalid operation exception and the contents of the
2037   destination parts are unspecified.  If the rounded value is in
2038   range but the floating point number is not the exact integer, the C
2039   standard doesn't require an inexact exception to be raised.  IEEE
2040   854 does require it so we do that.
2041
2042   Note that for conversions to integer type the C standard requires
2043   round-to-zero to always be used.  */
2044APFloat::opStatus
2045APFloat::convertToSignExtendedInteger(integerPart *parts, unsigned int width,
2046                                      bool isSigned,
2047                                      roundingMode rounding_mode,
2048                                      bool *isExact) const
2049{
2050  lostFraction lost_fraction;
2051  const integerPart *src;
2052  unsigned int dstPartsCount, truncatedBits;
2053
2054  *isExact = false;
2055
2056  /* Handle the three special cases first.  */
2057  if (category == fcInfinity || category == fcNaN)
2058    return opInvalidOp;
2059
2060  dstPartsCount = partCountForBits(width);
2061
2062  if (category == fcZero) {
2063    APInt::tcSet(parts, 0, dstPartsCount);
2064    // Negative zero can't be represented as an int.
2065    *isExact = !sign;
2066    return opOK;
2067  }
2068
2069  src = significandParts();
2070
2071  /* Step 1: place our absolute value, with any fraction truncated, in
2072     the destination.  */
2073  if (exponent < 0) {
2074    /* Our absolute value is less than one; truncate everything.  */
2075    APInt::tcSet(parts, 0, dstPartsCount);
2076    /* For exponent -1 the integer bit represents .5, look at that.
2077       For smaller exponents leftmost truncated bit is 0. */
2078    truncatedBits = semantics->precision -1U - exponent;
2079  } else {
2080    /* We want the most significant (exponent + 1) bits; the rest are
2081       truncated.  */
2082    unsigned int bits = exponent + 1U;
2083
2084    /* Hopelessly large in magnitude?  */
2085    if (bits > width)
2086      return opInvalidOp;
2087
2088    if (bits < semantics->precision) {
2089      /* We truncate (semantics->precision - bits) bits.  */
2090      truncatedBits = semantics->precision - bits;
2091      APInt::tcExtract(parts, dstPartsCount, src, bits, truncatedBits);
2092    } else {
2093      /* We want at least as many bits as are available.  */
2094      APInt::tcExtract(parts, dstPartsCount, src, semantics->precision, 0);
2095      APInt::tcShiftLeft(parts, dstPartsCount, bits - semantics->precision);
2096      truncatedBits = 0;
2097    }
2098  }
2099
2100  /* Step 2: work out any lost fraction, and increment the absolute
2101     value if we would round away from zero.  */
2102  if (truncatedBits) {
2103    lost_fraction = lostFractionThroughTruncation(src, partCount(),
2104                                                  truncatedBits);
2105    if (lost_fraction != lfExactlyZero &&
2106        roundAwayFromZero(rounding_mode, lost_fraction, truncatedBits)) {
2107      if (APInt::tcIncrement(parts, dstPartsCount))
2108        return opInvalidOp;     /* Overflow.  */
2109    }
2110  } else {
2111    lost_fraction = lfExactlyZero;
2112  }
2113
2114  /* Step 3: check if we fit in the destination.  */
2115  unsigned int omsb = APInt::tcMSB(parts, dstPartsCount) + 1;
2116
2117  if (sign) {
2118    if (!isSigned) {
2119      /* Negative numbers cannot be represented as unsigned.  */
2120      if (omsb != 0)
2121        return opInvalidOp;
2122    } else {
2123      /* It takes omsb bits to represent the unsigned integer value.
2124         We lose a bit for the sign, but care is needed as the
2125         maximally negative integer is a special case.  */
2126      if (omsb == width && APInt::tcLSB(parts, dstPartsCount) + 1 != omsb)
2127        return opInvalidOp;
2128
2129      /* This case can happen because of rounding.  */
2130      if (omsb > width)
2131        return opInvalidOp;
2132    }
2133
2134    APInt::tcNegate (parts, dstPartsCount);
2135  } else {
2136    if (omsb >= width + !isSigned)
2137      return opInvalidOp;
2138  }
2139
2140  if (lost_fraction == lfExactlyZero) {
2141    *isExact = true;
2142    return opOK;
2143  } else
2144    return opInexact;
2145}
2146
2147/* Same as convertToSignExtendedInteger, except we provide
2148   deterministic values in case of an invalid operation exception,
2149   namely zero for NaNs and the minimal or maximal value respectively
2150   for underflow or overflow.
2151   The *isExact output tells whether the result is exact, in the sense
2152   that converting it back to the original floating point type produces
2153   the original value.  This is almost equivalent to result==opOK,
2154   except for negative zeroes.
2155*/
2156APFloat::opStatus
2157APFloat::convertToInteger(integerPart *parts, unsigned int width,
2158                          bool isSigned,
2159                          roundingMode rounding_mode, bool *isExact) const
2160{
2161  opStatus fs;
2162
2163  fs = convertToSignExtendedInteger(parts, width, isSigned, rounding_mode,
2164                                    isExact);
2165
2166  if (fs == opInvalidOp) {
2167    unsigned int bits, dstPartsCount;
2168
2169    dstPartsCount = partCountForBits(width);
2170
2171    if (category == fcNaN)
2172      bits = 0;
2173    else if (sign)
2174      bits = isSigned;
2175    else
2176      bits = width - isSigned;
2177
2178    APInt::tcSetLeastSignificantBits(parts, dstPartsCount, bits);
2179    if (sign && isSigned)
2180      APInt::tcShiftLeft(parts, dstPartsCount, width - 1);
2181  }
2182
2183  return fs;
2184}
2185
2186/* Same as convertToInteger(integerPart*, ...), except the result is returned in
2187   an APSInt, whose initial bit-width and signed-ness are used to determine the
2188   precision of the conversion.
2189 */
2190APFloat::opStatus
2191APFloat::convertToInteger(APSInt &result,
2192                          roundingMode rounding_mode, bool *isExact) const
2193{
2194  unsigned bitWidth = result.getBitWidth();
2195  SmallVector<uint64_t, 4> parts(result.getNumWords());
2196  opStatus status = convertToInteger(
2197    parts.data(), bitWidth, result.isSigned(), rounding_mode, isExact);
2198  // Keeps the original signed-ness.
2199  result = APInt(bitWidth, parts);
2200  return status;
2201}
2202
2203/* Convert an unsigned integer SRC to a floating point number,
2204   rounding according to ROUNDING_MODE.  The sign of the floating
2205   point number is not modified.  */
2206APFloat::opStatus
2207APFloat::convertFromUnsignedParts(const integerPart *src,
2208                                  unsigned int srcCount,
2209                                  roundingMode rounding_mode)
2210{
2211  unsigned int omsb, precision, dstCount;
2212  integerPart *dst;
2213  lostFraction lost_fraction;
2214
2215  category = fcNormal;
2216  omsb = APInt::tcMSB(src, srcCount) + 1;
2217  dst = significandParts();
2218  dstCount = partCount();
2219  precision = semantics->precision;
2220
2221  /* We want the most significant PRECISION bits of SRC.  There may not
2222     be that many; extract what we can.  */
2223  if (precision <= omsb) {
2224    exponent = omsb - 1;
2225    lost_fraction = lostFractionThroughTruncation(src, srcCount,
2226                                                  omsb - precision);
2227    APInt::tcExtract(dst, dstCount, src, precision, omsb - precision);
2228  } else {
2229    exponent = precision - 1;
2230    lost_fraction = lfExactlyZero;
2231    APInt::tcExtract(dst, dstCount, src, omsb, 0);
2232  }
2233
2234  return normalize(rounding_mode, lost_fraction);
2235}
2236
2237APFloat::opStatus
2238APFloat::convertFromAPInt(const APInt &Val,
2239                          bool isSigned,
2240                          roundingMode rounding_mode)
2241{
2242  unsigned int partCount = Val.getNumWords();
2243  APInt api = Val;
2244
2245  sign = false;
2246  if (isSigned && api.isNegative()) {
2247    sign = true;
2248    api = -api;
2249  }
2250
2251  return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode);
2252}
2253
2254/* Convert a two's complement integer SRC to a floating point number,
2255   rounding according to ROUNDING_MODE.  ISSIGNED is true if the
2256   integer is signed, in which case it must be sign-extended.  */
2257APFloat::opStatus
2258APFloat::convertFromSignExtendedInteger(const integerPart *src,
2259                                        unsigned int srcCount,
2260                                        bool isSigned,
2261                                        roundingMode rounding_mode)
2262{
2263  opStatus status;
2264
2265  if (isSigned &&
2266      APInt::tcExtractBit(src, srcCount * integerPartWidth - 1)) {
2267    integerPart *copy;
2268
2269    /* If we're signed and negative negate a copy.  */
2270    sign = true;
2271    copy = new integerPart[srcCount];
2272    APInt::tcAssign(copy, src, srcCount);
2273    APInt::tcNegate(copy, srcCount);
2274    status = convertFromUnsignedParts(copy, srcCount, rounding_mode);
2275    delete [] copy;
2276  } else {
2277    sign = false;
2278    status = convertFromUnsignedParts(src, srcCount, rounding_mode);
2279  }
2280
2281  return status;
2282}
2283
2284/* FIXME: should this just take a const APInt reference?  */
2285APFloat::opStatus
2286APFloat::convertFromZeroExtendedInteger(const integerPart *parts,
2287                                        unsigned int width, bool isSigned,
2288                                        roundingMode rounding_mode)
2289{
2290  unsigned int partCount = partCountForBits(width);
2291  APInt api = APInt(width, makeArrayRef(parts, partCount));
2292
2293  sign = false;
2294  if (isSigned && APInt::tcExtractBit(parts, width - 1)) {
2295    sign = true;
2296    api = -api;
2297  }
2298
2299  return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode);
2300}
2301
2302APFloat::opStatus
2303APFloat::convertFromHexadecimalString(StringRef s, roundingMode rounding_mode)
2304{
2305  lostFraction lost_fraction = lfExactlyZero;
2306
2307  category = fcNormal;
2308  zeroSignificand();
2309  exponent = 0;
2310
2311  integerPart *significand = significandParts();
2312  unsigned partsCount = partCount();
2313  unsigned bitPos = partsCount * integerPartWidth;
2314  bool computedTrailingFraction = false;
2315
2316  // Skip leading zeroes and any (hexa)decimal point.
2317  StringRef::iterator begin = s.begin();
2318  StringRef::iterator end = s.end();
2319  StringRef::iterator dot;
2320  StringRef::iterator p = skipLeadingZeroesAndAnyDot(begin, end, &dot);
2321  StringRef::iterator firstSignificantDigit = p;
2322
2323  while (p != end) {
2324    integerPart hex_value;
2325
2326    if (*p == '.') {
2327      assert(dot == end && "String contains multiple dots");
2328      dot = p++;
2329      continue;
2330    }
2331
2332    hex_value = hexDigitValue(*p);
2333    if (hex_value == -1U)
2334      break;
2335
2336    p++;
2337
2338    // Store the number while we have space.
2339    if (bitPos) {
2340      bitPos -= 4;
2341      hex_value <<= bitPos % integerPartWidth;
2342      significand[bitPos / integerPartWidth] |= hex_value;
2343    } else if (!computedTrailingFraction) {
2344      lost_fraction = trailingHexadecimalFraction(p, end, hex_value);
2345      computedTrailingFraction = true;
2346    }
2347  }
2348
2349  /* Hex floats require an exponent but not a hexadecimal point.  */
2350  assert(p != end && "Hex strings require an exponent");
2351  assert((*p == 'p' || *p == 'P') && "Invalid character in significand");
2352  assert(p != begin && "Significand has no digits");
2353  assert((dot == end || p - begin != 1) && "Significand has no digits");
2354
2355  /* Ignore the exponent if we are zero.  */
2356  if (p != firstSignificantDigit) {
2357    int expAdjustment;
2358
2359    /* Implicit hexadecimal point?  */
2360    if (dot == end)
2361      dot = p;
2362
2363    /* Calculate the exponent adjustment implicit in the number of
2364       significant digits.  */
2365    expAdjustment = static_cast<int>(dot - firstSignificantDigit);
2366    if (expAdjustment < 0)
2367      expAdjustment++;
2368    expAdjustment = expAdjustment * 4 - 1;
2369
2370    /* Adjust for writing the significand starting at the most
2371       significant nibble.  */
2372    expAdjustment += semantics->precision;
2373    expAdjustment -= partsCount * integerPartWidth;
2374
2375    /* Adjust for the given exponent.  */
2376    exponent = totalExponent(p + 1, end, expAdjustment);
2377  }
2378
2379  return normalize(rounding_mode, lost_fraction);
2380}
2381
2382APFloat::opStatus
2383APFloat::roundSignificandWithExponent(const integerPart *decSigParts,
2384                                      unsigned sigPartCount, int exp,
2385                                      roundingMode rounding_mode)
2386{
2387  unsigned int parts, pow5PartCount;
2388  fltSemantics calcSemantics = { 32767, -32767, 0 };
2389  integerPart pow5Parts[maxPowerOfFiveParts];
2390  bool isNearest;
2391
2392  isNearest = (rounding_mode == rmNearestTiesToEven ||
2393               rounding_mode == rmNearestTiesToAway);
2394
2395  parts = partCountForBits(semantics->precision + 11);
2396
2397  /* Calculate pow(5, abs(exp)).  */
2398  pow5PartCount = powerOf5(pow5Parts, exp >= 0 ? exp: -exp);
2399
2400  for (;; parts *= 2) {
2401    opStatus sigStatus, powStatus;
2402    unsigned int excessPrecision, truncatedBits;
2403
2404    calcSemantics.precision = parts * integerPartWidth - 1;
2405    excessPrecision = calcSemantics.precision - semantics->precision;
2406    truncatedBits = excessPrecision;
2407
2408    APFloat decSig = APFloat::getZero(calcSemantics, sign);
2409    APFloat pow5(calcSemantics);
2410
2411    sigStatus = decSig.convertFromUnsignedParts(decSigParts, sigPartCount,
2412                                                rmNearestTiesToEven);
2413    powStatus = pow5.convertFromUnsignedParts(pow5Parts, pow5PartCount,
2414                                              rmNearestTiesToEven);
2415    /* Add exp, as 10^n = 5^n * 2^n.  */
2416    decSig.exponent += exp;
2417
2418    lostFraction calcLostFraction;
2419    integerPart HUerr, HUdistance;
2420    unsigned int powHUerr;
2421
2422    if (exp >= 0) {
2423      /* multiplySignificand leaves the precision-th bit set to 1.  */
2424      calcLostFraction = decSig.multiplySignificand(pow5, NULL);
2425      powHUerr = powStatus != opOK;
2426    } else {
2427      calcLostFraction = decSig.divideSignificand(pow5);
2428      /* Denormal numbers have less precision.  */
2429      if (decSig.exponent < semantics->minExponent) {
2430        excessPrecision += (semantics->minExponent - decSig.exponent);
2431        truncatedBits = excessPrecision;
2432        if (excessPrecision > calcSemantics.precision)
2433          excessPrecision = calcSemantics.precision;
2434      }
2435      /* Extra half-ulp lost in reciprocal of exponent.  */
2436      powHUerr = (powStatus == opOK && calcLostFraction == lfExactlyZero) ? 0:2;
2437    }
2438
2439    /* Both multiplySignificand and divideSignificand return the
2440       result with the integer bit set.  */
2441    assert(APInt::tcExtractBit
2442           (decSig.significandParts(), calcSemantics.precision - 1) == 1);
2443
2444    HUerr = HUerrBound(calcLostFraction != lfExactlyZero, sigStatus != opOK,
2445                       powHUerr);
2446    HUdistance = 2 * ulpsFromBoundary(decSig.significandParts(),
2447                                      excessPrecision, isNearest);
2448
2449    /* Are we guaranteed to round correctly if we truncate?  */
2450    if (HUdistance >= HUerr) {
2451      APInt::tcExtract(significandParts(), partCount(), decSig.significandParts(),
2452                       calcSemantics.precision - excessPrecision,
2453                       excessPrecision);
2454      /* Take the exponent of decSig.  If we tcExtract-ed less bits
2455         above we must adjust our exponent to compensate for the
2456         implicit right shift.  */
2457      exponent = (decSig.exponent + semantics->precision
2458                  - (calcSemantics.precision - excessPrecision));
2459      calcLostFraction = lostFractionThroughTruncation(decSig.significandParts(),
2460                                                       decSig.partCount(),
2461                                                       truncatedBits);
2462      return normalize(rounding_mode, calcLostFraction);
2463    }
2464  }
2465}
2466
2467APFloat::opStatus
2468APFloat::convertFromDecimalString(StringRef str, roundingMode rounding_mode)
2469{
2470  decimalInfo D;
2471  opStatus fs;
2472
2473  /* Scan the text.  */
2474  StringRef::iterator p = str.begin();
2475  interpretDecimal(p, str.end(), &D);
2476
2477  /* Handle the quick cases.  First the case of no significant digits,
2478     i.e. zero, and then exponents that are obviously too large or too
2479     small.  Writing L for log 10 / log 2, a number d.ddddd*10^exp
2480     definitely overflows if
2481
2482           (exp - 1) * L >= maxExponent
2483
2484     and definitely underflows to zero where
2485
2486           (exp + 1) * L <= minExponent - precision
2487
2488     With integer arithmetic the tightest bounds for L are
2489
2490           93/28 < L < 196/59            [ numerator <= 256 ]
2491           42039/12655 < L < 28738/8651  [ numerator <= 65536 ]
2492  */
2493
2494  // Test if we have a zero number allowing for strings with no null terminators
2495  // and zero decimals with non-zero exponents.
2496  //
2497  // We computed firstSigDigit by ignoring all zeros and dots. Thus if
2498  // D->firstSigDigit equals str.end(), every digit must be a zero and there can
2499  // be at most one dot. On the other hand, if we have a zero with a non-zero
2500  // exponent, then we know that D.firstSigDigit will be non-numeric.
2501  if (D.firstSigDigit == str.end() || decDigitValue(*D.firstSigDigit) >= 10U) {
2502    category = fcZero;
2503    fs = opOK;
2504
2505  /* Check whether the normalized exponent is high enough to overflow
2506     max during the log-rebasing in the max-exponent check below. */
2507  } else if (D.normalizedExponent - 1 > INT_MAX / 42039) {
2508    fs = handleOverflow(rounding_mode);
2509
2510  /* If it wasn't, then it also wasn't high enough to overflow max
2511     during the log-rebasing in the min-exponent check.  Check that it
2512     won't overflow min in either check, then perform the min-exponent
2513     check. */
2514  } else if (D.normalizedExponent - 1 < INT_MIN / 42039 ||
2515             (D.normalizedExponent + 1) * 28738 <=
2516               8651 * (semantics->minExponent - (int) semantics->precision)) {
2517    /* Underflow to zero and round.  */
2518    category = fcNormal;
2519    zeroSignificand();
2520    fs = normalize(rounding_mode, lfLessThanHalf);
2521
2522  /* We can finally safely perform the max-exponent check. */
2523  } else if ((D.normalizedExponent - 1) * 42039
2524             >= 12655 * semantics->maxExponent) {
2525    /* Overflow and round.  */
2526    fs = handleOverflow(rounding_mode);
2527  } else {
2528    integerPart *decSignificand;
2529    unsigned int partCount;
2530
2531    /* A tight upper bound on number of bits required to hold an
2532       N-digit decimal integer is N * 196 / 59.  Allocate enough space
2533       to hold the full significand, and an extra part required by
2534       tcMultiplyPart.  */
2535    partCount = static_cast<unsigned int>(D.lastSigDigit - D.firstSigDigit) + 1;
2536    partCount = partCountForBits(1 + 196 * partCount / 59);
2537    decSignificand = new integerPart[partCount + 1];
2538    partCount = 0;
2539
2540    /* Convert to binary efficiently - we do almost all multiplication
2541       in an integerPart.  When this would overflow do we do a single
2542       bignum multiplication, and then revert again to multiplication
2543       in an integerPart.  */
2544    do {
2545      integerPart decValue, val, multiplier;
2546
2547      val = 0;
2548      multiplier = 1;
2549
2550      do {
2551        if (*p == '.') {
2552          p++;
2553          if (p == str.end()) {
2554            break;
2555          }
2556        }
2557        decValue = decDigitValue(*p++);
2558        assert(decValue < 10U && "Invalid character in significand");
2559        multiplier *= 10;
2560        val = val * 10 + decValue;
2561        /* The maximum number that can be multiplied by ten with any
2562           digit added without overflowing an integerPart.  */
2563      } while (p <= D.lastSigDigit && multiplier <= (~ (integerPart) 0 - 9) / 10);
2564
2565      /* Multiply out the current part.  */
2566      APInt::tcMultiplyPart(decSignificand, decSignificand, multiplier, val,
2567                            partCount, partCount + 1, false);
2568
2569      /* If we used another part (likely but not guaranteed), increase
2570         the count.  */
2571      if (decSignificand[partCount])
2572        partCount++;
2573    } while (p <= D.lastSigDigit);
2574
2575    category = fcNormal;
2576    fs = roundSignificandWithExponent(decSignificand, partCount,
2577                                      D.exponent, rounding_mode);
2578
2579    delete [] decSignificand;
2580  }
2581
2582  return fs;
2583}
2584
2585bool
2586APFloat::convertFromStringSpecials(StringRef str) {
2587  if (str.equals("inf") || str.equals("INFINITY")) {
2588    makeInf(false);
2589    return true;
2590  }
2591
2592  if (str.equals("-inf") || str.equals("-INFINITY")) {
2593    makeInf(true);
2594    return true;
2595  }
2596
2597  if (str.equals("nan") || str.equals("NaN")) {
2598    makeNaN(false, false);
2599    return true;
2600  }
2601
2602  if (str.equals("-nan") || str.equals("-NaN")) {
2603    makeNaN(false, true);
2604    return true;
2605  }
2606
2607  return false;
2608}
2609
2610APFloat::opStatus
2611APFloat::convertFromString(StringRef str, roundingMode rounding_mode)
2612{
2613  assert(!str.empty() && "Invalid string length");
2614
2615  // Handle special cases.
2616  if (convertFromStringSpecials(str))
2617    return opOK;
2618
2619  /* Handle a leading minus sign.  */
2620  StringRef::iterator p = str.begin();
2621  size_t slen = str.size();
2622  sign = *p == '-' ? 1 : 0;
2623  if (*p == '-' || *p == '+') {
2624    p++;
2625    slen--;
2626    assert(slen && "String has no digits");
2627  }
2628
2629  if (slen >= 2 && p[0] == '0' && (p[1] == 'x' || p[1] == 'X')) {
2630    assert(slen - 2 && "Invalid string");
2631    return convertFromHexadecimalString(StringRef(p + 2, slen - 2),
2632                                        rounding_mode);
2633  }
2634
2635  return convertFromDecimalString(StringRef(p, slen), rounding_mode);
2636}
2637
2638/* Write out a hexadecimal representation of the floating point value
2639   to DST, which must be of sufficient size, in the C99 form
2640   [-]0xh.hhhhp[+-]d.  Return the number of characters written,
2641   excluding the terminating NUL.
2642
2643   If UPPERCASE, the output is in upper case, otherwise in lower case.
2644
2645   HEXDIGITS digits appear altogether, rounding the value if
2646   necessary.  If HEXDIGITS is 0, the minimal precision to display the
2647   number precisely is used instead.  If nothing would appear after
2648   the decimal point it is suppressed.
2649
2650   The decimal exponent is always printed and has at least one digit.
2651   Zero values display an exponent of zero.  Infinities and NaNs
2652   appear as "infinity" or "nan" respectively.
2653
2654   The above rules are as specified by C99.  There is ambiguity about
2655   what the leading hexadecimal digit should be.  This implementation
2656   uses whatever is necessary so that the exponent is displayed as
2657   stored.  This implies the exponent will fall within the IEEE format
2658   range, and the leading hexadecimal digit will be 0 (for denormals),
2659   1 (normal numbers) or 2 (normal numbers rounded-away-from-zero with
2660   any other digits zero).
2661*/
2662unsigned int
2663APFloat::convertToHexString(char *dst, unsigned int hexDigits,
2664                            bool upperCase, roundingMode rounding_mode) const
2665{
2666  char *p;
2667
2668  p = dst;
2669  if (sign)
2670    *dst++ = '-';
2671
2672  switch (category) {
2673  case fcInfinity:
2674    memcpy (dst, upperCase ? infinityU: infinityL, sizeof infinityU - 1);
2675    dst += sizeof infinityL - 1;
2676    break;
2677
2678  case fcNaN:
2679    memcpy (dst, upperCase ? NaNU: NaNL, sizeof NaNU - 1);
2680    dst += sizeof NaNU - 1;
2681    break;
2682
2683  case fcZero:
2684    *dst++ = '0';
2685    *dst++ = upperCase ? 'X': 'x';
2686    *dst++ = '0';
2687    if (hexDigits > 1) {
2688      *dst++ = '.';
2689      memset (dst, '0', hexDigits - 1);
2690      dst += hexDigits - 1;
2691    }
2692    *dst++ = upperCase ? 'P': 'p';
2693    *dst++ = '0';
2694    break;
2695
2696  case fcNormal:
2697    dst = convertNormalToHexString (dst, hexDigits, upperCase, rounding_mode);
2698    break;
2699  }
2700
2701  *dst = 0;
2702
2703  return static_cast<unsigned int>(dst - p);
2704}
2705
2706/* Does the hard work of outputting the correctly rounded hexadecimal
2707   form of a normal floating point number with the specified number of
2708   hexadecimal digits.  If HEXDIGITS is zero the minimum number of
2709   digits necessary to print the value precisely is output.  */
2710char *
2711APFloat::convertNormalToHexString(char *dst, unsigned int hexDigits,
2712                                  bool upperCase,
2713                                  roundingMode rounding_mode) const
2714{
2715  unsigned int count, valueBits, shift, partsCount, outputDigits;
2716  const char *hexDigitChars;
2717  const integerPart *significand;
2718  char *p;
2719  bool roundUp;
2720
2721  *dst++ = '0';
2722  *dst++ = upperCase ? 'X': 'x';
2723
2724  roundUp = false;
2725  hexDigitChars = upperCase ? hexDigitsUpper: hexDigitsLower;
2726
2727  significand = significandParts();
2728  partsCount = partCount();
2729
2730  /* +3 because the first digit only uses the single integer bit, so
2731     we have 3 virtual zero most-significant-bits.  */
2732  valueBits = semantics->precision + 3;
2733  shift = integerPartWidth - valueBits % integerPartWidth;
2734
2735  /* The natural number of digits required ignoring trailing
2736     insignificant zeroes.  */
2737  outputDigits = (valueBits - significandLSB () + 3) / 4;
2738
2739  /* hexDigits of zero means use the required number for the
2740     precision.  Otherwise, see if we are truncating.  If we are,
2741     find out if we need to round away from zero.  */
2742  if (hexDigits) {
2743    if (hexDigits < outputDigits) {
2744      /* We are dropping non-zero bits, so need to check how to round.
2745         "bits" is the number of dropped bits.  */
2746      unsigned int bits;
2747      lostFraction fraction;
2748
2749      bits = valueBits - hexDigits * 4;
2750      fraction = lostFractionThroughTruncation (significand, partsCount, bits);
2751      roundUp = roundAwayFromZero(rounding_mode, fraction, bits);
2752    }
2753    outputDigits = hexDigits;
2754  }
2755
2756  /* Write the digits consecutively, and start writing in the location
2757     of the hexadecimal point.  We move the most significant digit
2758     left and add the hexadecimal point later.  */
2759  p = ++dst;
2760
2761  count = (valueBits + integerPartWidth - 1) / integerPartWidth;
2762
2763  while (outputDigits && count) {
2764    integerPart part;
2765
2766    /* Put the most significant integerPartWidth bits in "part".  */
2767    if (--count == partsCount)
2768      part = 0;  /* An imaginary higher zero part.  */
2769    else
2770      part = significand[count] << shift;
2771
2772    if (count && shift)
2773      part |= significand[count - 1] >> (integerPartWidth - shift);
2774
2775    /* Convert as much of "part" to hexdigits as we can.  */
2776    unsigned int curDigits = integerPartWidth / 4;
2777
2778    if (curDigits > outputDigits)
2779      curDigits = outputDigits;
2780    dst += partAsHex (dst, part, curDigits, hexDigitChars);
2781    outputDigits -= curDigits;
2782  }
2783
2784  if (roundUp) {
2785    char *q = dst;
2786
2787    /* Note that hexDigitChars has a trailing '0'.  */
2788    do {
2789      q--;
2790      *q = hexDigitChars[hexDigitValue (*q) + 1];
2791    } while (*q == '0');
2792    assert(q >= p);
2793  } else {
2794    /* Add trailing zeroes.  */
2795    memset (dst, '0', outputDigits);
2796    dst += outputDigits;
2797  }
2798
2799  /* Move the most significant digit to before the point, and if there
2800     is something after the decimal point add it.  This must come
2801     after rounding above.  */
2802  p[-1] = p[0];
2803  if (dst -1 == p)
2804    dst--;
2805  else
2806    p[0] = '.';
2807
2808  /* Finally output the exponent.  */
2809  *dst++ = upperCase ? 'P': 'p';
2810
2811  return writeSignedDecimal (dst, exponent);
2812}
2813
2814hash_code llvm::hash_value(const APFloat &Arg) {
2815  if (!Arg.isFiniteNonZero())
2816    return hash_combine((uint8_t)Arg.category,
2817                        // NaN has no sign, fix it at zero.
2818                        Arg.isNaN() ? (uint8_t)0 : (uint8_t)Arg.sign,
2819                        Arg.semantics->precision);
2820
2821  // Normal floats need their exponent and significand hashed.
2822  return hash_combine((uint8_t)Arg.category, (uint8_t)Arg.sign,
2823                      Arg.semantics->precision, Arg.exponent,
2824                      hash_combine_range(
2825                        Arg.significandParts(),
2826                        Arg.significandParts() + Arg.partCount()));
2827}
2828
2829// Conversion from APFloat to/from host float/double.  It may eventually be
2830// possible to eliminate these and have everybody deal with APFloats, but that
2831// will take a while.  This approach will not easily extend to long double.
2832// Current implementation requires integerPartWidth==64, which is correct at
2833// the moment but could be made more general.
2834
2835// Denormals have exponent minExponent in APFloat, but minExponent-1 in
2836// the actual IEEE respresentations.  We compensate for that here.
2837
2838APInt
2839APFloat::convertF80LongDoubleAPFloatToAPInt() const
2840{
2841  assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended);
2842  assert(partCount()==2);
2843
2844  uint64_t myexponent, mysignificand;
2845
2846  if (isFiniteNonZero()) {
2847    myexponent = exponent+16383; //bias
2848    mysignificand = significandParts()[0];
2849    if (myexponent==1 && !(mysignificand & 0x8000000000000000ULL))
2850      myexponent = 0;   // denormal
2851  } else if (category==fcZero) {
2852    myexponent = 0;
2853    mysignificand = 0;
2854  } else if (category==fcInfinity) {
2855    myexponent = 0x7fff;
2856    mysignificand = 0x8000000000000000ULL;
2857  } else {
2858    assert(category == fcNaN && "Unknown category");
2859    myexponent = 0x7fff;
2860    mysignificand = significandParts()[0];
2861  }
2862
2863  uint64_t words[2];
2864  words[0] = mysignificand;
2865  words[1] =  ((uint64_t)(sign & 1) << 15) |
2866              (myexponent & 0x7fffLL);
2867  return APInt(80, words);
2868}
2869
2870APInt
2871APFloat::convertPPCDoubleDoubleAPFloatToAPInt() const
2872{
2873  assert(semantics == (const llvm::fltSemantics*)&PPCDoubleDouble);
2874  assert(partCount()==2);
2875
2876  uint64_t words[2];
2877  opStatus fs;
2878  bool losesInfo;
2879
2880  // Convert number to double.  To avoid spurious underflows, we re-
2881  // normalize against the "double" minExponent first, and only *then*
2882  // truncate the mantissa.  The result of that second conversion
2883  // may be inexact, but should never underflow.
2884  // Declare fltSemantics before APFloat that uses it (and
2885  // saves pointer to it) to ensure correct destruction order.
2886  fltSemantics extendedSemantics = *semantics;
2887  extendedSemantics.minExponent = IEEEdouble.minExponent;
2888  APFloat extended(*this);
2889  fs = extended.convert(extendedSemantics, rmNearestTiesToEven, &losesInfo);
2890  assert(fs == opOK && !losesInfo);
2891  (void)fs;
2892
2893  APFloat u(extended);
2894  fs = u.convert(IEEEdouble, rmNearestTiesToEven, &losesInfo);
2895  assert(fs == opOK || fs == opInexact);
2896  (void)fs;
2897  words[0] = *u.convertDoubleAPFloatToAPInt().getRawData();
2898
2899  // If conversion was exact or resulted in a special case, we're done;
2900  // just set the second double to zero.  Otherwise, re-convert back to
2901  // the extended format and compute the difference.  This now should
2902  // convert exactly to double.
2903  if (u.isFiniteNonZero() && losesInfo) {
2904    fs = u.convert(extendedSemantics, rmNearestTiesToEven, &losesInfo);
2905    assert(fs == opOK && !losesInfo);
2906    (void)fs;
2907
2908    APFloat v(extended);
2909    v.subtract(u, rmNearestTiesToEven);
2910    fs = v.convert(IEEEdouble, rmNearestTiesToEven, &losesInfo);
2911    assert(fs == opOK && !losesInfo);
2912    (void)fs;
2913    words[1] = *v.convertDoubleAPFloatToAPInt().getRawData();
2914  } else {
2915    words[1] = 0;
2916  }
2917
2918  return APInt(128, words);
2919}
2920
2921APInt
2922APFloat::convertQuadrupleAPFloatToAPInt() const
2923{
2924  assert(semantics == (const llvm::fltSemantics*)&IEEEquad);
2925  assert(partCount()==2);
2926
2927  uint64_t myexponent, mysignificand, mysignificand2;
2928
2929  if (isFiniteNonZero()) {
2930    myexponent = exponent+16383; //bias
2931    mysignificand = significandParts()[0];
2932    mysignificand2 = significandParts()[1];
2933    if (myexponent==1 && !(mysignificand2 & 0x1000000000000LL))
2934      myexponent = 0;   // denormal
2935  } else if (category==fcZero) {
2936    myexponent = 0;
2937    mysignificand = mysignificand2 = 0;
2938  } else if (category==fcInfinity) {
2939    myexponent = 0x7fff;
2940    mysignificand = mysignificand2 = 0;
2941  } else {
2942    assert(category == fcNaN && "Unknown category!");
2943    myexponent = 0x7fff;
2944    mysignificand = significandParts()[0];
2945    mysignificand2 = significandParts()[1];
2946  }
2947
2948  uint64_t words[2];
2949  words[0] = mysignificand;
2950  words[1] = ((uint64_t)(sign & 1) << 63) |
2951             ((myexponent & 0x7fff) << 48) |
2952             (mysignificand2 & 0xffffffffffffLL);
2953
2954  return APInt(128, words);
2955}
2956
2957APInt
2958APFloat::convertDoubleAPFloatToAPInt() const
2959{
2960  assert(semantics == (const llvm::fltSemantics*)&IEEEdouble);
2961  assert(partCount()==1);
2962
2963  uint64_t myexponent, mysignificand;
2964
2965  if (isFiniteNonZero()) {
2966    myexponent = exponent+1023; //bias
2967    mysignificand = *significandParts();
2968    if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
2969      myexponent = 0;   // denormal
2970  } else if (category==fcZero) {
2971    myexponent = 0;
2972    mysignificand = 0;
2973  } else if (category==fcInfinity) {
2974    myexponent = 0x7ff;
2975    mysignificand = 0;
2976  } else {
2977    assert(category == fcNaN && "Unknown category!");
2978    myexponent = 0x7ff;
2979    mysignificand = *significandParts();
2980  }
2981
2982  return APInt(64, ((((uint64_t)(sign & 1) << 63) |
2983                     ((myexponent & 0x7ff) <<  52) |
2984                     (mysignificand & 0xfffffffffffffLL))));
2985}
2986
2987APInt
2988APFloat::convertFloatAPFloatToAPInt() const
2989{
2990  assert(semantics == (const llvm::fltSemantics*)&IEEEsingle);
2991  assert(partCount()==1);
2992
2993  uint32_t myexponent, mysignificand;
2994
2995  if (isFiniteNonZero()) {
2996    myexponent = exponent+127; //bias
2997    mysignificand = (uint32_t)*significandParts();
2998    if (myexponent == 1 && !(mysignificand & 0x800000))
2999      myexponent = 0;   // denormal
3000  } else if (category==fcZero) {
3001    myexponent = 0;
3002    mysignificand = 0;
3003  } else if (category==fcInfinity) {
3004    myexponent = 0xff;
3005    mysignificand = 0;
3006  } else {
3007    assert(category == fcNaN && "Unknown category!");
3008    myexponent = 0xff;
3009    mysignificand = (uint32_t)*significandParts();
3010  }
3011
3012  return APInt(32, (((sign&1) << 31) | ((myexponent&0xff) << 23) |
3013                    (mysignificand & 0x7fffff)));
3014}
3015
3016APInt
3017APFloat::convertHalfAPFloatToAPInt() const
3018{
3019  assert(semantics == (const llvm::fltSemantics*)&IEEEhalf);
3020  assert(partCount()==1);
3021
3022  uint32_t myexponent, mysignificand;
3023
3024  if (isFiniteNonZero()) {
3025    myexponent = exponent+15; //bias
3026    mysignificand = (uint32_t)*significandParts();
3027    if (myexponent == 1 && !(mysignificand & 0x400))
3028      myexponent = 0;   // denormal
3029  } else if (category==fcZero) {
3030    myexponent = 0;
3031    mysignificand = 0;
3032  } else if (category==fcInfinity) {
3033    myexponent = 0x1f;
3034    mysignificand = 0;
3035  } else {
3036    assert(category == fcNaN && "Unknown category!");
3037    myexponent = 0x1f;
3038    mysignificand = (uint32_t)*significandParts();
3039  }
3040
3041  return APInt(16, (((sign&1) << 15) | ((myexponent&0x1f) << 10) |
3042                    (mysignificand & 0x3ff)));
3043}
3044
3045// This function creates an APInt that is just a bit map of the floating
3046// point constant as it would appear in memory.  It is not a conversion,
3047// and treating the result as a normal integer is unlikely to be useful.
3048
3049APInt
3050APFloat::bitcastToAPInt() const
3051{
3052  if (semantics == (const llvm::fltSemantics*)&IEEEhalf)
3053    return convertHalfAPFloatToAPInt();
3054
3055  if (semantics == (const llvm::fltSemantics*)&IEEEsingle)
3056    return convertFloatAPFloatToAPInt();
3057
3058  if (semantics == (const llvm::fltSemantics*)&IEEEdouble)
3059    return convertDoubleAPFloatToAPInt();
3060
3061  if (semantics == (const llvm::fltSemantics*)&IEEEquad)
3062    return convertQuadrupleAPFloatToAPInt();
3063
3064  if (semantics == (const llvm::fltSemantics*)&PPCDoubleDouble)
3065    return convertPPCDoubleDoubleAPFloatToAPInt();
3066
3067  assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended &&
3068         "unknown format!");
3069  return convertF80LongDoubleAPFloatToAPInt();
3070}
3071
3072float
3073APFloat::convertToFloat() const
3074{
3075  assert(semantics == (const llvm::fltSemantics*)&IEEEsingle &&
3076         "Float semantics are not IEEEsingle");
3077  APInt api = bitcastToAPInt();
3078  return api.bitsToFloat();
3079}
3080
3081double
3082APFloat::convertToDouble() const
3083{
3084  assert(semantics == (const llvm::fltSemantics*)&IEEEdouble &&
3085         "Float semantics are not IEEEdouble");
3086  APInt api = bitcastToAPInt();
3087  return api.bitsToDouble();
3088}
3089
3090/// Integer bit is explicit in this format.  Intel hardware (387 and later)
3091/// does not support these bit patterns:
3092///  exponent = all 1's, integer bit 0, significand 0 ("pseudoinfinity")
3093///  exponent = all 1's, integer bit 0, significand nonzero ("pseudoNaN")
3094///  exponent = 0, integer bit 1 ("pseudodenormal")
3095///  exponent!=0 nor all 1's, integer bit 0 ("unnormal")
3096/// At the moment, the first two are treated as NaNs, the second two as Normal.
3097void
3098APFloat::initFromF80LongDoubleAPInt(const APInt &api)
3099{
3100  assert(api.getBitWidth()==80);
3101  uint64_t i1 = api.getRawData()[0];
3102  uint64_t i2 = api.getRawData()[1];
3103  uint64_t myexponent = (i2 & 0x7fff);
3104  uint64_t mysignificand = i1;
3105
3106  initialize(&APFloat::x87DoubleExtended);
3107  assert(partCount()==2);
3108
3109  sign = static_cast<unsigned int>(i2>>15);
3110  if (myexponent==0 && mysignificand==0) {
3111    // exponent, significand meaningless
3112    category = fcZero;
3113  } else if (myexponent==0x7fff && mysignificand==0x8000000000000000ULL) {
3114    // exponent, significand meaningless
3115    category = fcInfinity;
3116  } else if (myexponent==0x7fff && mysignificand!=0x8000000000000000ULL) {
3117    // exponent meaningless
3118    category = fcNaN;
3119    significandParts()[0] = mysignificand;
3120    significandParts()[1] = 0;
3121  } else {
3122    category = fcNormal;
3123    exponent = myexponent - 16383;
3124    significandParts()[0] = mysignificand;
3125    significandParts()[1] = 0;
3126    if (myexponent==0)          // denormal
3127      exponent = -16382;
3128  }
3129}
3130
3131void
3132APFloat::initFromPPCDoubleDoubleAPInt(const APInt &api)
3133{
3134  assert(api.getBitWidth()==128);
3135  uint64_t i1 = api.getRawData()[0];
3136  uint64_t i2 = api.getRawData()[1];
3137  opStatus fs;
3138  bool losesInfo;
3139
3140  // Get the first double and convert to our format.
3141  initFromDoubleAPInt(APInt(64, i1));
3142  fs = convert(PPCDoubleDouble, rmNearestTiesToEven, &losesInfo);
3143  assert(fs == opOK && !losesInfo);
3144  (void)fs;
3145
3146  // Unless we have a special case, add in second double.
3147  if (isFiniteNonZero()) {
3148    APFloat v(IEEEdouble, APInt(64, i2));
3149    fs = v.convert(PPCDoubleDouble, rmNearestTiesToEven, &losesInfo);
3150    assert(fs == opOK && !losesInfo);
3151    (void)fs;
3152
3153    add(v, rmNearestTiesToEven);
3154  }
3155}
3156
3157void
3158APFloat::initFromQuadrupleAPInt(const APInt &api)
3159{
3160  assert(api.getBitWidth()==128);
3161  uint64_t i1 = api.getRawData()[0];
3162  uint64_t i2 = api.getRawData()[1];
3163  uint64_t myexponent = (i2 >> 48) & 0x7fff;
3164  uint64_t mysignificand  = i1;
3165  uint64_t mysignificand2 = i2 & 0xffffffffffffLL;
3166
3167  initialize(&APFloat::IEEEquad);
3168  assert(partCount()==2);
3169
3170  sign = static_cast<unsigned int>(i2>>63);
3171  if (myexponent==0 &&
3172      (mysignificand==0 && mysignificand2==0)) {
3173    // exponent, significand meaningless
3174    category = fcZero;
3175  } else if (myexponent==0x7fff &&
3176             (mysignificand==0 && mysignificand2==0)) {
3177    // exponent, significand meaningless
3178    category = fcInfinity;
3179  } else if (myexponent==0x7fff &&
3180             (mysignificand!=0 || mysignificand2 !=0)) {
3181    // exponent meaningless
3182    category = fcNaN;
3183    significandParts()[0] = mysignificand;
3184    significandParts()[1] = mysignificand2;
3185  } else {
3186    category = fcNormal;
3187    exponent = myexponent - 16383;
3188    significandParts()[0] = mysignificand;
3189    significandParts()[1] = mysignificand2;
3190    if (myexponent==0)          // denormal
3191      exponent = -16382;
3192    else
3193      significandParts()[1] |= 0x1000000000000LL;  // integer bit
3194  }
3195}
3196
3197void
3198APFloat::initFromDoubleAPInt(const APInt &api)
3199{
3200  assert(api.getBitWidth()==64);
3201  uint64_t i = *api.getRawData();
3202  uint64_t myexponent = (i >> 52) & 0x7ff;
3203  uint64_t mysignificand = i & 0xfffffffffffffLL;
3204
3205  initialize(&APFloat::IEEEdouble);
3206  assert(partCount()==1);
3207
3208  sign = static_cast<unsigned int>(i>>63);
3209  if (myexponent==0 && mysignificand==0) {
3210    // exponent, significand meaningless
3211    category = fcZero;
3212  } else if (myexponent==0x7ff && mysignificand==0) {
3213    // exponent, significand meaningless
3214    category = fcInfinity;
3215  } else if (myexponent==0x7ff && mysignificand!=0) {
3216    // exponent meaningless
3217    category = fcNaN;
3218    *significandParts() = mysignificand;
3219  } else {
3220    category = fcNormal;
3221    exponent = myexponent - 1023;
3222    *significandParts() = mysignificand;
3223    if (myexponent==0)          // denormal
3224      exponent = -1022;
3225    else
3226      *significandParts() |= 0x10000000000000LL;  // integer bit
3227  }
3228}
3229
3230void
3231APFloat::initFromFloatAPInt(const APInt & api)
3232{
3233  assert(api.getBitWidth()==32);
3234  uint32_t i = (uint32_t)*api.getRawData();
3235  uint32_t myexponent = (i >> 23) & 0xff;
3236  uint32_t mysignificand = i & 0x7fffff;
3237
3238  initialize(&APFloat::IEEEsingle);
3239  assert(partCount()==1);
3240
3241  sign = i >> 31;
3242  if (myexponent==0 && mysignificand==0) {
3243    // exponent, significand meaningless
3244    category = fcZero;
3245  } else if (myexponent==0xff && mysignificand==0) {
3246    // exponent, significand meaningless
3247    category = fcInfinity;
3248  } else if (myexponent==0xff && mysignificand!=0) {
3249    // sign, exponent, significand meaningless
3250    category = fcNaN;
3251    *significandParts() = mysignificand;
3252  } else {
3253    category = fcNormal;
3254    exponent = myexponent - 127;  //bias
3255    *significandParts() = mysignificand;
3256    if (myexponent==0)    // denormal
3257      exponent = -126;
3258    else
3259      *significandParts() |= 0x800000; // integer bit
3260  }
3261}
3262
3263void
3264APFloat::initFromHalfAPInt(const APInt & api)
3265{
3266  assert(api.getBitWidth()==16);
3267  uint32_t i = (uint32_t)*api.getRawData();
3268  uint32_t myexponent = (i >> 10) & 0x1f;
3269  uint32_t mysignificand = i & 0x3ff;
3270
3271  initialize(&APFloat::IEEEhalf);
3272  assert(partCount()==1);
3273
3274  sign = i >> 15;
3275  if (myexponent==0 && mysignificand==0) {
3276    // exponent, significand meaningless
3277    category = fcZero;
3278  } else if (myexponent==0x1f && mysignificand==0) {
3279    // exponent, significand meaningless
3280    category = fcInfinity;
3281  } else if (myexponent==0x1f && mysignificand!=0) {
3282    // sign, exponent, significand meaningless
3283    category = fcNaN;
3284    *significandParts() = mysignificand;
3285  } else {
3286    category = fcNormal;
3287    exponent = myexponent - 15;  //bias
3288    *significandParts() = mysignificand;
3289    if (myexponent==0)    // denormal
3290      exponent = -14;
3291    else
3292      *significandParts() |= 0x400; // integer bit
3293  }
3294}
3295
3296/// Treat api as containing the bits of a floating point number.  Currently
3297/// we infer the floating point type from the size of the APInt.  The
3298/// isIEEE argument distinguishes between PPC128 and IEEE128 (not meaningful
3299/// when the size is anything else).
3300void
3301APFloat::initFromAPInt(const fltSemantics* Sem, const APInt& api)
3302{
3303  if (Sem == &IEEEhalf)
3304    return initFromHalfAPInt(api);
3305  if (Sem == &IEEEsingle)
3306    return initFromFloatAPInt(api);
3307  if (Sem == &IEEEdouble)
3308    return initFromDoubleAPInt(api);
3309  if (Sem == &x87DoubleExtended)
3310    return initFromF80LongDoubleAPInt(api);
3311  if (Sem == &IEEEquad)
3312    return initFromQuadrupleAPInt(api);
3313  if (Sem == &PPCDoubleDouble)
3314    return initFromPPCDoubleDoubleAPInt(api);
3315
3316  llvm_unreachable(0);
3317}
3318
3319APFloat
3320APFloat::getAllOnesValue(unsigned BitWidth, bool isIEEE)
3321{
3322  switch (BitWidth) {
3323  case 16:
3324    return APFloat(IEEEhalf, APInt::getAllOnesValue(BitWidth));
3325  case 32:
3326    return APFloat(IEEEsingle, APInt::getAllOnesValue(BitWidth));
3327  case 64:
3328    return APFloat(IEEEdouble, APInt::getAllOnesValue(BitWidth));
3329  case 80:
3330    return APFloat(x87DoubleExtended, APInt::getAllOnesValue(BitWidth));
3331  case 128:
3332    if (isIEEE)
3333      return APFloat(IEEEquad, APInt::getAllOnesValue(BitWidth));
3334    return APFloat(PPCDoubleDouble, APInt::getAllOnesValue(BitWidth));
3335  default:
3336    llvm_unreachable("Unknown floating bit width");
3337  }
3338}
3339
3340/// Make this number the largest magnitude normal number in the given
3341/// semantics.
3342void APFloat::makeLargest(bool Negative) {
3343  // We want (in interchange format):
3344  //   sign = {Negative}
3345  //   exponent = 1..10
3346  //   significand = 1..1
3347  category = fcNormal;
3348  sign = Negative;
3349  exponent = semantics->maxExponent;
3350
3351  // Use memset to set all but the highest integerPart to all ones.
3352  integerPart *significand = significandParts();
3353  unsigned PartCount = partCount();
3354  memset(significand, 0xFF, sizeof(integerPart)*(PartCount - 1));
3355
3356  // Set the high integerPart especially setting all unused top bits for
3357  // internal consistency.
3358  const unsigned NumUnusedHighBits =
3359    PartCount*integerPartWidth - semantics->precision;
3360  significand[PartCount - 1] = ~integerPart(0) >> NumUnusedHighBits;
3361}
3362
3363/// Make this number the smallest magnitude denormal number in the given
3364/// semantics.
3365void APFloat::makeSmallest(bool Negative) {
3366  // We want (in interchange format):
3367  //   sign = {Negative}
3368  //   exponent = 0..0
3369  //   significand = 0..01
3370  category = fcNormal;
3371  sign = Negative;
3372  exponent = semantics->minExponent;
3373  APInt::tcSet(significandParts(), 1, partCount());
3374}
3375
3376
3377APFloat APFloat::getLargest(const fltSemantics &Sem, bool Negative) {
3378  // We want (in interchange format):
3379  //   sign = {Negative}
3380  //   exponent = 1..10
3381  //   significand = 1..1
3382  APFloat Val(Sem, uninitialized);
3383  Val.makeLargest(Negative);
3384  return Val;
3385}
3386
3387APFloat APFloat::getSmallest(const fltSemantics &Sem, bool Negative) {
3388  // We want (in interchange format):
3389  //   sign = {Negative}
3390  //   exponent = 0..0
3391  //   significand = 0..01
3392  APFloat Val(Sem, uninitialized);
3393  Val.makeSmallest(Negative);
3394  return Val;
3395}
3396
3397APFloat APFloat::getSmallestNormalized(const fltSemantics &Sem, bool Negative) {
3398  APFloat Val(Sem, uninitialized);
3399
3400  // We want (in interchange format):
3401  //   sign = {Negative}
3402  //   exponent = 0..0
3403  //   significand = 10..0
3404
3405  Val.category = fcNormal;
3406  Val.zeroSignificand();
3407  Val.sign = Negative;
3408  Val.exponent = Sem.minExponent;
3409  Val.significandParts()[partCountForBits(Sem.precision)-1] |=
3410    (((integerPart) 1) << ((Sem.precision - 1) % integerPartWidth));
3411
3412  return Val;
3413}
3414
3415APFloat::APFloat(const fltSemantics &Sem, const APInt &API) {
3416  initFromAPInt(&Sem, API);
3417}
3418
3419APFloat::APFloat(float f) {
3420  initFromAPInt(&IEEEsingle, APInt::floatToBits(f));
3421}
3422
3423APFloat::APFloat(double d) {
3424  initFromAPInt(&IEEEdouble, APInt::doubleToBits(d));
3425}
3426
3427namespace {
3428  void append(SmallVectorImpl<char> &Buffer, StringRef Str) {
3429    Buffer.append(Str.begin(), Str.end());
3430  }
3431
3432  /// Removes data from the given significand until it is no more
3433  /// precise than is required for the desired precision.
3434  void AdjustToPrecision(APInt &significand,
3435                         int &exp, unsigned FormatPrecision) {
3436    unsigned bits = significand.getActiveBits();
3437
3438    // 196/59 is a very slight overestimate of lg_2(10).
3439    unsigned bitsRequired = (FormatPrecision * 196 + 58) / 59;
3440
3441    if (bits <= bitsRequired) return;
3442
3443    unsigned tensRemovable = (bits - bitsRequired) * 59 / 196;
3444    if (!tensRemovable) return;
3445
3446    exp += tensRemovable;
3447
3448    APInt divisor(significand.getBitWidth(), 1);
3449    APInt powten(significand.getBitWidth(), 10);
3450    while (true) {
3451      if (tensRemovable & 1)
3452        divisor *= powten;
3453      tensRemovable >>= 1;
3454      if (!tensRemovable) break;
3455      powten *= powten;
3456    }
3457
3458    significand = significand.udiv(divisor);
3459
3460    // Truncate the significand down to its active bit count.
3461    significand = significand.trunc(significand.getActiveBits());
3462  }
3463
3464
3465  void AdjustToPrecision(SmallVectorImpl<char> &buffer,
3466                         int &exp, unsigned FormatPrecision) {
3467    unsigned N = buffer.size();
3468    if (N <= FormatPrecision) return;
3469
3470    // The most significant figures are the last ones in the buffer.
3471    unsigned FirstSignificant = N - FormatPrecision;
3472
3473    // Round.
3474    // FIXME: this probably shouldn't use 'round half up'.
3475
3476    // Rounding down is just a truncation, except we also want to drop
3477    // trailing zeros from the new result.
3478    if (buffer[FirstSignificant - 1] < '5') {
3479      while (FirstSignificant < N && buffer[FirstSignificant] == '0')
3480        FirstSignificant++;
3481
3482      exp += FirstSignificant;
3483      buffer.erase(&buffer[0], &buffer[FirstSignificant]);
3484      return;
3485    }
3486
3487    // Rounding up requires a decimal add-with-carry.  If we continue
3488    // the carry, the newly-introduced zeros will just be truncated.
3489    for (unsigned I = FirstSignificant; I != N; ++I) {
3490      if (buffer[I] == '9') {
3491        FirstSignificant++;
3492      } else {
3493        buffer[I]++;
3494        break;
3495      }
3496    }
3497
3498    // If we carried through, we have exactly one digit of precision.
3499    if (FirstSignificant == N) {
3500      exp += FirstSignificant;
3501      buffer.clear();
3502      buffer.push_back('1');
3503      return;
3504    }
3505
3506    exp += FirstSignificant;
3507    buffer.erase(&buffer[0], &buffer[FirstSignificant]);
3508  }
3509}
3510
3511void APFloat::toString(SmallVectorImpl<char> &Str,
3512                       unsigned FormatPrecision,
3513                       unsigned FormatMaxPadding) const {
3514  switch (category) {
3515  case fcInfinity:
3516    if (isNegative())
3517      return append(Str, "-Inf");
3518    else
3519      return append(Str, "+Inf");
3520
3521  case fcNaN: return append(Str, "NaN");
3522
3523  case fcZero:
3524    if (isNegative())
3525      Str.push_back('-');
3526
3527    if (!FormatMaxPadding)
3528      append(Str, "0.0E+0");
3529    else
3530      Str.push_back('0');
3531    return;
3532
3533  case fcNormal:
3534    break;
3535  }
3536
3537  if (isNegative())
3538    Str.push_back('-');
3539
3540  // Decompose the number into an APInt and an exponent.
3541  int exp = exponent - ((int) semantics->precision - 1);
3542  APInt significand(semantics->precision,
3543                    makeArrayRef(significandParts(),
3544                                 partCountForBits(semantics->precision)));
3545
3546  // Set FormatPrecision if zero.  We want to do this before we
3547  // truncate trailing zeros, as those are part of the precision.
3548  if (!FormatPrecision) {
3549    // We use enough digits so the number can be round-tripped back to an
3550    // APFloat. The formula comes from "How to Print Floating-Point Numbers
3551    // Accurately" by Steele and White.
3552    // FIXME: Using a formula based purely on the precision is conservative;
3553    // we can print fewer digits depending on the actual value being printed.
3554
3555    // FormatPrecision = 2 + floor(significandBits / lg_2(10))
3556    FormatPrecision = 2 + semantics->precision * 59 / 196;
3557  }
3558
3559  // Ignore trailing binary zeros.
3560  int trailingZeros = significand.countTrailingZeros();
3561  exp += trailingZeros;
3562  significand = significand.lshr(trailingZeros);
3563
3564  // Change the exponent from 2^e to 10^e.
3565  if (exp == 0) {
3566    // Nothing to do.
3567  } else if (exp > 0) {
3568    // Just shift left.
3569    significand = significand.zext(semantics->precision + exp);
3570    significand <<= exp;
3571    exp = 0;
3572  } else { /* exp < 0 */
3573    int texp = -exp;
3574
3575    // We transform this using the identity:
3576    //   (N)(2^-e) == (N)(5^e)(10^-e)
3577    // This means we have to multiply N (the significand) by 5^e.
3578    // To avoid overflow, we have to operate on numbers large
3579    // enough to store N * 5^e:
3580    //   log2(N * 5^e) == log2(N) + e * log2(5)
3581    //                 <= semantics->precision + e * 137 / 59
3582    //   (log_2(5) ~ 2.321928 < 2.322034 ~ 137/59)
3583
3584    unsigned precision = semantics->precision + (137 * texp + 136) / 59;
3585
3586    // Multiply significand by 5^e.
3587    //   N * 5^0101 == N * 5^(1*1) * 5^(0*2) * 5^(1*4) * 5^(0*8)
3588    significand = significand.zext(precision);
3589    APInt five_to_the_i(precision, 5);
3590    while (true) {
3591      if (texp & 1) significand *= five_to_the_i;
3592
3593      texp >>= 1;
3594      if (!texp) break;
3595      five_to_the_i *= five_to_the_i;
3596    }
3597  }
3598
3599  AdjustToPrecision(significand, exp, FormatPrecision);
3600
3601  SmallVector<char, 256> buffer;
3602
3603  // Fill the buffer.
3604  unsigned precision = significand.getBitWidth();
3605  APInt ten(precision, 10);
3606  APInt digit(precision, 0);
3607
3608  bool inTrail = true;
3609  while (significand != 0) {
3610    // digit <- significand % 10
3611    // significand <- significand / 10
3612    APInt::udivrem(significand, ten, significand, digit);
3613
3614    unsigned d = digit.getZExtValue();
3615
3616    // Drop trailing zeros.
3617    if (inTrail && !d) exp++;
3618    else {
3619      buffer.push_back((char) ('0' + d));
3620      inTrail = false;
3621    }
3622  }
3623
3624  assert(!buffer.empty() && "no characters in buffer!");
3625
3626  // Drop down to FormatPrecision.
3627  // TODO: don't do more precise calculations above than are required.
3628  AdjustToPrecision(buffer, exp, FormatPrecision);
3629
3630  unsigned NDigits = buffer.size();
3631
3632  // Check whether we should use scientific notation.
3633  bool FormatScientific;
3634  if (!FormatMaxPadding)
3635    FormatScientific = true;
3636  else {
3637    if (exp >= 0) {
3638      // 765e3 --> 765000
3639      //              ^^^
3640      // But we shouldn't make the number look more precise than it is.
3641      FormatScientific = ((unsigned) exp > FormatMaxPadding ||
3642                          NDigits + (unsigned) exp > FormatPrecision);
3643    } else {
3644      // Power of the most significant digit.
3645      int MSD = exp + (int) (NDigits - 1);
3646      if (MSD >= 0) {
3647        // 765e-2 == 7.65
3648        FormatScientific = false;
3649      } else {
3650        // 765e-5 == 0.00765
3651        //           ^ ^^
3652        FormatScientific = ((unsigned) -MSD) > FormatMaxPadding;
3653      }
3654    }
3655  }
3656
3657  // Scientific formatting is pretty straightforward.
3658  if (FormatScientific) {
3659    exp += (NDigits - 1);
3660
3661    Str.push_back(buffer[NDigits-1]);
3662    Str.push_back('.');
3663    if (NDigits == 1)
3664      Str.push_back('0');
3665    else
3666      for (unsigned I = 1; I != NDigits; ++I)
3667        Str.push_back(buffer[NDigits-1-I]);
3668    Str.push_back('E');
3669
3670    Str.push_back(exp >= 0 ? '+' : '-');
3671    if (exp < 0) exp = -exp;
3672    SmallVector<char, 6> expbuf;
3673    do {
3674      expbuf.push_back((char) ('0' + (exp % 10)));
3675      exp /= 10;
3676    } while (exp);
3677    for (unsigned I = 0, E = expbuf.size(); I != E; ++I)
3678      Str.push_back(expbuf[E-1-I]);
3679    return;
3680  }
3681
3682  // Non-scientific, positive exponents.
3683  if (exp >= 0) {
3684    for (unsigned I = 0; I != NDigits; ++I)
3685      Str.push_back(buffer[NDigits-1-I]);
3686    for (unsigned I = 0; I != (unsigned) exp; ++I)
3687      Str.push_back('0');
3688    return;
3689  }
3690
3691  // Non-scientific, negative exponents.
3692
3693  // The number of digits to the left of the decimal point.
3694  int NWholeDigits = exp + (int) NDigits;
3695
3696  unsigned I = 0;
3697  if (NWholeDigits > 0) {
3698    for (; I != (unsigned) NWholeDigits; ++I)
3699      Str.push_back(buffer[NDigits-I-1]);
3700    Str.push_back('.');
3701  } else {
3702    unsigned NZeros = 1 + (unsigned) -NWholeDigits;
3703
3704    Str.push_back('0');
3705    Str.push_back('.');
3706    for (unsigned Z = 1; Z != NZeros; ++Z)
3707      Str.push_back('0');
3708  }
3709
3710  for (; I != NDigits; ++I)
3711    Str.push_back(buffer[NDigits-I-1]);
3712}
3713
3714bool APFloat::getExactInverse(APFloat *inv) const {
3715  // Special floats and denormals have no exact inverse.
3716  if (!isFiniteNonZero())
3717    return false;
3718
3719  // Check that the number is a power of two by making sure that only the
3720  // integer bit is set in the significand.
3721  if (significandLSB() != semantics->precision - 1)
3722    return false;
3723
3724  // Get the inverse.
3725  APFloat reciprocal(*semantics, 1ULL);
3726  if (reciprocal.divide(*this, rmNearestTiesToEven) != opOK)
3727    return false;
3728
3729  // Avoid multiplication with a denormal, it is not safe on all platforms and
3730  // may be slower than a normal division.
3731  if (reciprocal.isDenormal())
3732    return false;
3733
3734  assert(reciprocal.isFiniteNonZero() &&
3735         reciprocal.significandLSB() == reciprocal.semantics->precision - 1);
3736
3737  if (inv)
3738    *inv = reciprocal;
3739
3740  return true;
3741}
3742
3743bool APFloat::isSignaling() const {
3744  if (!isNaN())
3745    return false;
3746
3747  // IEEE-754R 2008 6.2.1: A signaling NaN bit string should be encoded with the
3748  // first bit of the trailing significand being 0.
3749  return !APInt::tcExtractBit(significandParts(), semantics->precision - 2);
3750}
3751
3752/// IEEE-754R 2008 5.3.1: nextUp/nextDown.
3753///
3754/// *NOTE* since nextDown(x) = -nextUp(-x), we only implement nextUp with
3755/// appropriate sign switching before/after the computation.
3756APFloat::opStatus APFloat::next(bool nextDown) {
3757  // If we are performing nextDown, swap sign so we have -x.
3758  if (nextDown)
3759    changeSign();
3760
3761  // Compute nextUp(x)
3762  opStatus result = opOK;
3763
3764  // Handle each float category separately.
3765  switch (category) {
3766  case fcInfinity:
3767    // nextUp(+inf) = +inf
3768    if (!isNegative())
3769      break;
3770    // nextUp(-inf) = -getLargest()
3771    makeLargest(true);
3772    break;
3773  case fcNaN:
3774    // IEEE-754R 2008 6.2 Par 2: nextUp(sNaN) = qNaN. Set Invalid flag.
3775    // IEEE-754R 2008 6.2: nextUp(qNaN) = qNaN. Must be identity so we do not
3776    //                     change the payload.
3777    if (isSignaling()) {
3778      result = opInvalidOp;
3779      // For consistency, propogate the sign of the sNaN to the qNaN.
3780      makeNaN(false, isNegative(), 0);
3781    }
3782    break;
3783  case fcZero:
3784    // nextUp(pm 0) = +getSmallest()
3785    makeSmallest(false);
3786    break;
3787  case fcNormal:
3788    // nextUp(-getSmallest()) = -0
3789    if (isSmallest() && isNegative()) {
3790      APInt::tcSet(significandParts(), 0, partCount());
3791      category = fcZero;
3792      exponent = 0;
3793      break;
3794    }
3795
3796    // nextUp(getLargest()) == INFINITY
3797    if (isLargest() && !isNegative()) {
3798      APInt::tcSet(significandParts(), 0, partCount());
3799      category = fcInfinity;
3800      exponent = semantics->maxExponent + 1;
3801      break;
3802    }
3803
3804    // nextUp(normal) == normal + inc.
3805    if (isNegative()) {
3806      // If we are negative, we need to decrement the significand.
3807
3808      // We only cross a binade boundary that requires adjusting the exponent
3809      // if:
3810      //   1. exponent != semantics->minExponent. This implies we are not in the
3811      //   smallest binade or are dealing with denormals.
3812      //   2. Our significand excluding the integral bit is all zeros.
3813      bool WillCrossBinadeBoundary =
3814        exponent != semantics->minExponent && isSignificandAllZeros();
3815
3816      // Decrement the significand.
3817      //
3818      // We always do this since:
3819      //   1. If we are dealing with a non binade decrement, by definition we
3820      //   just decrement the significand.
3821      //   2. If we are dealing with a normal -> normal binade decrement, since
3822      //   we have an explicit integral bit the fact that all bits but the
3823      //   integral bit are zero implies that subtracting one will yield a
3824      //   significand with 0 integral bit and 1 in all other spots. Thus we
3825      //   must just adjust the exponent and set the integral bit to 1.
3826      //   3. If we are dealing with a normal -> denormal binade decrement,
3827      //   since we set the integral bit to 0 when we represent denormals, we
3828      //   just decrement the significand.
3829      integerPart *Parts = significandParts();
3830      APInt::tcDecrement(Parts, partCount());
3831
3832      if (WillCrossBinadeBoundary) {
3833        // Our result is a normal number. Do the following:
3834        // 1. Set the integral bit to 1.
3835        // 2. Decrement the exponent.
3836        APInt::tcSetBit(Parts, semantics->precision - 1);
3837        exponent--;
3838      }
3839    } else {
3840      // If we are positive, we need to increment the significand.
3841
3842      // We only cross a binade boundary that requires adjusting the exponent if
3843      // the input is not a denormal and all of said input's significand bits
3844      // are set. If all of said conditions are true: clear the significand, set
3845      // the integral bit to 1, and increment the exponent. If we have a
3846      // denormal always increment since moving denormals and the numbers in the
3847      // smallest normal binade have the same exponent in our representation.
3848      bool WillCrossBinadeBoundary = !isDenormal() && isSignificandAllOnes();
3849
3850      if (WillCrossBinadeBoundary) {
3851        integerPart *Parts = significandParts();
3852        APInt::tcSet(Parts, 0, partCount());
3853        APInt::tcSetBit(Parts, semantics->precision - 1);
3854        assert(exponent != semantics->maxExponent &&
3855               "We can not increment an exponent beyond the maxExponent allowed"
3856               " by the given floating point semantics.");
3857        exponent++;
3858      } else {
3859        incrementSignificand();
3860      }
3861    }
3862    break;
3863  }
3864
3865  // If we are performing nextDown, swap sign so we have -nextUp(-x)
3866  if (nextDown)
3867    changeSign();
3868
3869  return result;
3870}
3871
3872void
3873APFloat::makeInf(bool Negative) {
3874  category = fcInfinity;
3875  sign = Negative;
3876  exponent = semantics->maxExponent + 1;
3877  APInt::tcSet(significandParts(), 0, partCount());
3878}
3879
3880void
3881APFloat::makeZero(bool Negative) {
3882  category = fcZero;
3883  sign = Negative;
3884  exponent = semantics->minExponent-1;
3885  APInt::tcSet(significandParts(), 0, partCount());
3886}
3887