1//===-- APFloat.cpp - Implement APFloat class -----------------------------===//
2//
3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4// See https://llvm.org/LICENSE.txt for license information.
5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6//
7//===----------------------------------------------------------------------===//
8//
9// This file implements a class to represent arbitrary precision floating
10// point values and provide a variety of arithmetic operations on them.
11//
12//===----------------------------------------------------------------------===//
13
14#include "llvm/ADT/APFloat.h"
15#include "llvm/ADT/APSInt.h"
16#include "llvm/ADT/ArrayRef.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/Config/llvm-config.h"
22#include "llvm/Support/Debug.h"
23#include "llvm/Support/Error.h"
24#include "llvm/Support/MathExtras.h"
25#include "llvm/Support/raw_ostream.h"
26#include <cstring>
27#include <limits.h>
28
29#define APFLOAT_DISPATCH_ON_SEMANTICS(METHOD_CALL)                             \
30  do {                                                                         \
31    if (usesLayout<IEEEFloat>(getSemantics()))                                 \
32      return U.IEEE.METHOD_CALL;                                               \
33    if (usesLayout<DoubleAPFloat>(getSemantics()))                             \
34      return U.Double.METHOD_CALL;                                             \
35    llvm_unreachable("Unexpected semantics");                                  \
36  } while (false)
37
38using namespace llvm;
39
40/// A macro used to combine two fcCategory enums into one key which can be used
41/// in a switch statement to classify how the interaction of two APFloat's
42/// categories affects an operation.
43///
44/// TODO: If clang source code is ever allowed to use constexpr in its own
45/// codebase, change this into a static inline function.
46#define PackCategoriesIntoKey(_lhs, _rhs) ((_lhs) * 4 + (_rhs))
47
48/* Assumed in hexadecimal significand parsing, and conversion to
49   hexadecimal strings.  */
50static_assert(APFloatBase::integerPartWidth % 4 == 0, "Part width must be divisible by 4!");
51
52namespace llvm {
53  /* Represents floating point arithmetic semantics.  */
54  struct fltSemantics {
55    /* The largest E such that 2^E is representable; this matches the
56       definition of IEEE 754.  */
57    APFloatBase::ExponentType maxExponent;
58
59    /* The smallest E such that 2^E is a normalized number; this
60       matches the definition of IEEE 754.  */
61    APFloatBase::ExponentType minExponent;
62
63    /* Number of bits in the significand.  This includes the integer
64       bit.  */
65    unsigned int precision;
66
67    /* Number of bits actually used in the semantics. */
68    unsigned int sizeInBits;
69  };
70
71  static const fltSemantics semIEEEhalf = {15, -14, 11, 16};
72  static const fltSemantics semBFloat = {127, -126, 8, 16};
73  static const fltSemantics semIEEEsingle = {127, -126, 24, 32};
74  static const fltSemantics semIEEEdouble = {1023, -1022, 53, 64};
75  static const fltSemantics semIEEEquad = {16383, -16382, 113, 128};
76  static const fltSemantics semX87DoubleExtended = {16383, -16382, 64, 80};
77  static const fltSemantics semBogus = {0, 0, 0, 0};
78
79  /* The IBM double-double semantics. Such a number consists of a pair of IEEE
80     64-bit doubles (Hi, Lo), where |Hi| > |Lo|, and if normal,
81     (double)(Hi + Lo) == Hi. The numeric value it's modeling is Hi + Lo.
82     Therefore it has two 53-bit mantissa parts that aren't necessarily adjacent
83     to each other, and two 11-bit exponents.
84
85     Note: we need to make the value different from semBogus as otherwise
86     an unsafe optimization may collapse both values to a single address,
87     and we heavily rely on them having distinct addresses.             */
88  static const fltSemantics semPPCDoubleDouble = {-1, 0, 0, 0};
89
90  /* These are legacy semantics for the fallback, inaccrurate implementation of
91     IBM double-double, if the accurate semPPCDoubleDouble doesn't handle the
92     operation. It's equivalent to having an IEEE number with consecutive 106
93     bits of mantissa and 11 bits of exponent.
94
95     It's not equivalent to IBM double-double. For example, a legit IBM
96     double-double, 1 + epsilon:
97
98       1 + epsilon = 1 + (1 >> 1076)
99
100     is not representable by a consecutive 106 bits of mantissa.
101
102     Currently, these semantics are used in the following way:
103
104       semPPCDoubleDouble -> (IEEEdouble, IEEEdouble) ->
105       (64-bit APInt, 64-bit APInt) -> (128-bit APInt) ->
106       semPPCDoubleDoubleLegacy -> IEEE operations
107
108     We use bitcastToAPInt() to get the bit representation (in APInt) of the
109     underlying IEEEdouble, then use the APInt constructor to construct the
110     legacy IEEE float.
111
112     TODO: Implement all operations in semPPCDoubleDouble, and delete these
113     semantics.  */
114  static const fltSemantics semPPCDoubleDoubleLegacy = {1023, -1022 + 53,
115                                                        53 + 53, 128};
116
117  const llvm::fltSemantics &APFloatBase::EnumToSemantics(Semantics S) {
118    switch (S) {
119    case S_IEEEhalf:
120      return IEEEhalf();
121    case S_BFloat:
122      return BFloat();
123    case S_IEEEsingle:
124      return IEEEsingle();
125    case S_IEEEdouble:
126      return IEEEdouble();
127    case S_x87DoubleExtended:
128      return x87DoubleExtended();
129    case S_IEEEquad:
130      return IEEEquad();
131    case S_PPCDoubleDouble:
132      return PPCDoubleDouble();
133    }
134    llvm_unreachable("Unrecognised floating semantics");
135  }
136
137  APFloatBase::Semantics
138  APFloatBase::SemanticsToEnum(const llvm::fltSemantics &Sem) {
139    if (&Sem == &llvm::APFloat::IEEEhalf())
140      return S_IEEEhalf;
141    else if (&Sem == &llvm::APFloat::BFloat())
142      return S_BFloat;
143    else if (&Sem == &llvm::APFloat::IEEEsingle())
144      return S_IEEEsingle;
145    else if (&Sem == &llvm::APFloat::IEEEdouble())
146      return S_IEEEdouble;
147    else if (&Sem == &llvm::APFloat::x87DoubleExtended())
148      return S_x87DoubleExtended;
149    else if (&Sem == &llvm::APFloat::IEEEquad())
150      return S_IEEEquad;
151    else if (&Sem == &llvm::APFloat::PPCDoubleDouble())
152      return S_PPCDoubleDouble;
153    else
154      llvm_unreachable("Unknown floating semantics");
155  }
156
157  const fltSemantics &APFloatBase::IEEEhalf() {
158    return semIEEEhalf;
159  }
160  const fltSemantics &APFloatBase::BFloat() {
161    return semBFloat;
162  }
163  const fltSemantics &APFloatBase::IEEEsingle() {
164    return semIEEEsingle;
165  }
166  const fltSemantics &APFloatBase::IEEEdouble() {
167    return semIEEEdouble;
168  }
169  const fltSemantics &APFloatBase::IEEEquad() {
170    return semIEEEquad;
171  }
172  const fltSemantics &APFloatBase::x87DoubleExtended() {
173    return semX87DoubleExtended;
174  }
175  const fltSemantics &APFloatBase::Bogus() {
176    return semBogus;
177  }
178  const fltSemantics &APFloatBase::PPCDoubleDouble() {
179    return semPPCDoubleDouble;
180  }
181
182  constexpr RoundingMode APFloatBase::rmNearestTiesToEven;
183  constexpr RoundingMode APFloatBase::rmTowardPositive;
184  constexpr RoundingMode APFloatBase::rmTowardNegative;
185  constexpr RoundingMode APFloatBase::rmTowardZero;
186  constexpr RoundingMode APFloatBase::rmNearestTiesToAway;
187
188  /* A tight upper bound on number of parts required to hold the value
189     pow(5, power) is
190
191       power * 815 / (351 * integerPartWidth) + 1
192
193     However, whilst the result may require only this many parts,
194     because we are multiplying two values to get it, the
195     multiplication may require an extra part with the excess part
196     being zero (consider the trivial case of 1 * 1, tcFullMultiply
197     requires two parts to hold the single-part result).  So we add an
198     extra one to guarantee enough space whilst multiplying.  */
199  const unsigned int maxExponent = 16383;
200  const unsigned int maxPrecision = 113;
201  const unsigned int maxPowerOfFiveExponent = maxExponent + maxPrecision - 1;
202  const unsigned int maxPowerOfFiveParts = 2 + ((maxPowerOfFiveExponent * 815) / (351 * APFloatBase::integerPartWidth));
203
204  unsigned int APFloatBase::semanticsPrecision(const fltSemantics &semantics) {
205    return semantics.precision;
206  }
207  APFloatBase::ExponentType
208  APFloatBase::semanticsMaxExponent(const fltSemantics &semantics) {
209    return semantics.maxExponent;
210  }
211  APFloatBase::ExponentType
212  APFloatBase::semanticsMinExponent(const fltSemantics &semantics) {
213    return semantics.minExponent;
214  }
215  unsigned int APFloatBase::semanticsSizeInBits(const fltSemantics &semantics) {
216    return semantics.sizeInBits;
217  }
218
219  unsigned APFloatBase::getSizeInBits(const fltSemantics &Sem) {
220    return Sem.sizeInBits;
221}
222
223/* A bunch of private, handy routines.  */
224
225static inline Error createError(const Twine &Err) {
226  return make_error<StringError>(Err, inconvertibleErrorCode());
227}
228
229static inline unsigned int
230partCountForBits(unsigned int bits)
231{
232  return ((bits) + APFloatBase::integerPartWidth - 1) / APFloatBase::integerPartWidth;
233}
234
235/* Returns 0U-9U.  Return values >= 10U are not digits.  */
236static inline unsigned int
237decDigitValue(unsigned int c)
238{
239  return c - '0';
240}
241
242/* Return the value of a decimal exponent of the form
243   [+-]ddddddd.
244
245   If the exponent overflows, returns a large exponent with the
246   appropriate sign.  */
247static Expected<int> readExponent(StringRef::iterator begin,
248                                  StringRef::iterator end) {
249  bool isNegative;
250  unsigned int absExponent;
251  const unsigned int overlargeExponent = 24000;  /* FIXME.  */
252  StringRef::iterator p = begin;
253
254  // Treat no exponent as 0 to match binutils
255  if (p == end || ((*p == '-' || *p == '+') && (p + 1) == end)) {
256    return 0;
257  }
258
259  isNegative = (*p == '-');
260  if (*p == '-' || *p == '+') {
261    p++;
262    if (p == end)
263      return createError("Exponent has no digits");
264  }
265
266  absExponent = decDigitValue(*p++);
267  if (absExponent >= 10U)
268    return createError("Invalid character in exponent");
269
270  for (; p != end; ++p) {
271    unsigned int value;
272
273    value = decDigitValue(*p);
274    if (value >= 10U)
275      return createError("Invalid character in exponent");
276
277    absExponent = absExponent * 10U + value;
278    if (absExponent >= overlargeExponent) {
279      absExponent = overlargeExponent;
280      break;
281    }
282  }
283
284  if (isNegative)
285    return -(int) absExponent;
286  else
287    return (int) absExponent;
288}
289
290/* This is ugly and needs cleaning up, but I don't immediately see
291   how whilst remaining safe.  */
292static Expected<int> totalExponent(StringRef::iterator p,
293                                   StringRef::iterator end,
294                                   int exponentAdjustment) {
295  int unsignedExponent;
296  bool negative, overflow;
297  int exponent = 0;
298
299  if (p == end)
300    return createError("Exponent has no digits");
301
302  negative = *p == '-';
303  if (*p == '-' || *p == '+') {
304    p++;
305    if (p == end)
306      return createError("Exponent has no digits");
307  }
308
309  unsignedExponent = 0;
310  overflow = false;
311  for (; p != end; ++p) {
312    unsigned int value;
313
314    value = decDigitValue(*p);
315    if (value >= 10U)
316      return createError("Invalid character in exponent");
317
318    unsignedExponent = unsignedExponent * 10 + value;
319    if (unsignedExponent > 32767) {
320      overflow = true;
321      break;
322    }
323  }
324
325  if (exponentAdjustment > 32767 || exponentAdjustment < -32768)
326    overflow = true;
327
328  if (!overflow) {
329    exponent = unsignedExponent;
330    if (negative)
331      exponent = -exponent;
332    exponent += exponentAdjustment;
333    if (exponent > 32767 || exponent < -32768)
334      overflow = true;
335  }
336
337  if (overflow)
338    exponent = negative ? -32768: 32767;
339
340  return exponent;
341}
342
343static Expected<StringRef::iterator>
344skipLeadingZeroesAndAnyDot(StringRef::iterator begin, StringRef::iterator end,
345                           StringRef::iterator *dot) {
346  StringRef::iterator p = begin;
347  *dot = end;
348  while (p != end && *p == '0')
349    p++;
350
351  if (p != end && *p == '.') {
352    *dot = p++;
353
354    if (end - begin == 1)
355      return createError("Significand has no digits");
356
357    while (p != end && *p == '0')
358      p++;
359  }
360
361  return p;
362}
363
364/* Given a normal decimal floating point number of the form
365
366     dddd.dddd[eE][+-]ddd
367
368   where the decimal point and exponent are optional, fill out the
369   structure D.  Exponent is appropriate if the significand is
370   treated as an integer, and normalizedExponent if the significand
371   is taken to have the decimal point after a single leading
372   non-zero digit.
373
374   If the value is zero, V->firstSigDigit points to a non-digit, and
375   the return exponent is zero.
376*/
377struct decimalInfo {
378  const char *firstSigDigit;
379  const char *lastSigDigit;
380  int exponent;
381  int normalizedExponent;
382};
383
384static Error interpretDecimal(StringRef::iterator begin,
385                              StringRef::iterator end, decimalInfo *D) {
386  StringRef::iterator dot = end;
387
388  auto PtrOrErr = skipLeadingZeroesAndAnyDot(begin, end, &dot);
389  if (!PtrOrErr)
390    return PtrOrErr.takeError();
391  StringRef::iterator p = *PtrOrErr;
392
393  D->firstSigDigit = p;
394  D->exponent = 0;
395  D->normalizedExponent = 0;
396
397  for (; p != end; ++p) {
398    if (*p == '.') {
399      if (dot != end)
400        return createError("String contains multiple dots");
401      dot = p++;
402      if (p == end)
403        break;
404    }
405    if (decDigitValue(*p) >= 10U)
406      break;
407  }
408
409  if (p != end) {
410    if (*p != 'e' && *p != 'E')
411      return createError("Invalid character in significand");
412    if (p == begin)
413      return createError("Significand has no digits");
414    if (dot != end && p - begin == 1)
415      return createError("Significand has no digits");
416
417    /* p points to the first non-digit in the string */
418    auto ExpOrErr = readExponent(p + 1, end);
419    if (!ExpOrErr)
420      return ExpOrErr.takeError();
421    D->exponent = *ExpOrErr;
422
423    /* Implied decimal point?  */
424    if (dot == end)
425      dot = p;
426  }
427
428  /* If number is all zeroes accept any exponent.  */
429  if (p != D->firstSigDigit) {
430    /* Drop insignificant trailing zeroes.  */
431    if (p != begin) {
432      do
433        do
434          p--;
435        while (p != begin && *p == '0');
436      while (p != begin && *p == '.');
437    }
438
439    /* Adjust the exponents for any decimal point.  */
440    D->exponent += static_cast<APFloat::ExponentType>((dot - p) - (dot > p));
441    D->normalizedExponent = (D->exponent +
442              static_cast<APFloat::ExponentType>((p - D->firstSigDigit)
443                                      - (dot > D->firstSigDigit && dot < p)));
444  }
445
446  D->lastSigDigit = p;
447  return Error::success();
448}
449
450/* Return the trailing fraction of a hexadecimal number.
451   DIGITVALUE is the first hex digit of the fraction, P points to
452   the next digit.  */
453static Expected<lostFraction>
454trailingHexadecimalFraction(StringRef::iterator p, StringRef::iterator end,
455                            unsigned int digitValue) {
456  unsigned int hexDigit;
457
458  /* If the first trailing digit isn't 0 or 8 we can work out the
459     fraction immediately.  */
460  if (digitValue > 8)
461    return lfMoreThanHalf;
462  else if (digitValue < 8 && digitValue > 0)
463    return lfLessThanHalf;
464
465  // Otherwise we need to find the first non-zero digit.
466  while (p != end && (*p == '0' || *p == '.'))
467    p++;
468
469  if (p == end)
470    return createError("Invalid trailing hexadecimal fraction!");
471
472  hexDigit = hexDigitValue(*p);
473
474  /* If we ran off the end it is exactly zero or one-half, otherwise
475     a little more.  */
476  if (hexDigit == -1U)
477    return digitValue == 0 ? lfExactlyZero: lfExactlyHalf;
478  else
479    return digitValue == 0 ? lfLessThanHalf: lfMoreThanHalf;
480}
481
482/* Return the fraction lost were a bignum truncated losing the least
483   significant BITS bits.  */
484static lostFraction
485lostFractionThroughTruncation(const APFloatBase::integerPart *parts,
486                              unsigned int partCount,
487                              unsigned int bits)
488{
489  unsigned int lsb;
490
491  lsb = APInt::tcLSB(parts, partCount);
492
493  /* Note this is guaranteed true if bits == 0, or LSB == -1U.  */
494  if (bits <= lsb)
495    return lfExactlyZero;
496  if (bits == lsb + 1)
497    return lfExactlyHalf;
498  if (bits <= partCount * APFloatBase::integerPartWidth &&
499      APInt::tcExtractBit(parts, bits - 1))
500    return lfMoreThanHalf;
501
502  return lfLessThanHalf;
503}
504
505/* Shift DST right BITS bits noting lost fraction.  */
506static lostFraction
507shiftRight(APFloatBase::integerPart *dst, unsigned int parts, unsigned int bits)
508{
509  lostFraction lost_fraction;
510
511  lost_fraction = lostFractionThroughTruncation(dst, parts, bits);
512
513  APInt::tcShiftRight(dst, parts, bits);
514
515  return lost_fraction;
516}
517
518/* Combine the effect of two lost fractions.  */
519static lostFraction
520combineLostFractions(lostFraction moreSignificant,
521                     lostFraction lessSignificant)
522{
523  if (lessSignificant != lfExactlyZero) {
524    if (moreSignificant == lfExactlyZero)
525      moreSignificant = lfLessThanHalf;
526    else if (moreSignificant == lfExactlyHalf)
527      moreSignificant = lfMoreThanHalf;
528  }
529
530  return moreSignificant;
531}
532
533/* The error from the true value, in half-ulps, on multiplying two
534   floating point numbers, which differ from the value they
535   approximate by at most HUE1 and HUE2 half-ulps, is strictly less
536   than the returned value.
537
538   See "How to Read Floating Point Numbers Accurately" by William D
539   Clinger.  */
540static unsigned int
541HUerrBound(bool inexactMultiply, unsigned int HUerr1, unsigned int HUerr2)
542{
543  assert(HUerr1 < 2 || HUerr2 < 2 || (HUerr1 + HUerr2 < 8));
544
545  if (HUerr1 + HUerr2 == 0)
546    return inexactMultiply * 2;  /* <= inexactMultiply half-ulps.  */
547  else
548    return inexactMultiply + 2 * (HUerr1 + HUerr2);
549}
550
551/* The number of ulps from the boundary (zero, or half if ISNEAREST)
552   when the least significant BITS are truncated.  BITS cannot be
553   zero.  */
554static APFloatBase::integerPart
555ulpsFromBoundary(const APFloatBase::integerPart *parts, unsigned int bits,
556                 bool isNearest) {
557  unsigned int count, partBits;
558  APFloatBase::integerPart part, boundary;
559
560  assert(bits != 0);
561
562  bits--;
563  count = bits / APFloatBase::integerPartWidth;
564  partBits = bits % APFloatBase::integerPartWidth + 1;
565
566  part = parts[count] & (~(APFloatBase::integerPart) 0 >> (APFloatBase::integerPartWidth - partBits));
567
568  if (isNearest)
569    boundary = (APFloatBase::integerPart) 1 << (partBits - 1);
570  else
571    boundary = 0;
572
573  if (count == 0) {
574    if (part - boundary <= boundary - part)
575      return part - boundary;
576    else
577      return boundary - part;
578  }
579
580  if (part == boundary) {
581    while (--count)
582      if (parts[count])
583        return ~(APFloatBase::integerPart) 0; /* A lot.  */
584
585    return parts[0];
586  } else if (part == boundary - 1) {
587    while (--count)
588      if (~parts[count])
589        return ~(APFloatBase::integerPart) 0; /* A lot.  */
590
591    return -parts[0];
592  }
593
594  return ~(APFloatBase::integerPart) 0; /* A lot.  */
595}
596
597/* Place pow(5, power) in DST, and return the number of parts used.
598   DST must be at least one part larger than size of the answer.  */
599static unsigned int
600powerOf5(APFloatBase::integerPart *dst, unsigned int power) {
601  static const APFloatBase::integerPart firstEightPowers[] = { 1, 5, 25, 125, 625, 3125, 15625, 78125 };
602  APFloatBase::integerPart pow5s[maxPowerOfFiveParts * 2 + 5];
603  pow5s[0] = 78125 * 5;
604
605  unsigned int partsCount[16] = { 1 };
606  APFloatBase::integerPart scratch[maxPowerOfFiveParts], *p1, *p2, *pow5;
607  unsigned int result;
608  assert(power <= maxExponent);
609
610  p1 = dst;
611  p2 = scratch;
612
613  *p1 = firstEightPowers[power & 7];
614  power >>= 3;
615
616  result = 1;
617  pow5 = pow5s;
618
619  for (unsigned int n = 0; power; power >>= 1, n++) {
620    unsigned int pc;
621
622    pc = partsCount[n];
623
624    /* Calculate pow(5,pow(2,n+3)) if we haven't yet.  */
625    if (pc == 0) {
626      pc = partsCount[n - 1];
627      APInt::tcFullMultiply(pow5, pow5 - pc, pow5 - pc, pc, pc);
628      pc *= 2;
629      if (pow5[pc - 1] == 0)
630        pc--;
631      partsCount[n] = pc;
632    }
633
634    if (power & 1) {
635      APFloatBase::integerPart *tmp;
636
637      APInt::tcFullMultiply(p2, p1, pow5, result, pc);
638      result += pc;
639      if (p2[result - 1] == 0)
640        result--;
641
642      /* Now result is in p1 with partsCount parts and p2 is scratch
643         space.  */
644      tmp = p1;
645      p1 = p2;
646      p2 = tmp;
647    }
648
649    pow5 += pc;
650  }
651
652  if (p1 != dst)
653    APInt::tcAssign(dst, p1, result);
654
655  return result;
656}
657
658/* Zero at the end to avoid modular arithmetic when adding one; used
659   when rounding up during hexadecimal output.  */
660static const char hexDigitsLower[] = "0123456789abcdef0";
661static const char hexDigitsUpper[] = "0123456789ABCDEF0";
662static const char infinityL[] = "infinity";
663static const char infinityU[] = "INFINITY";
664static const char NaNL[] = "nan";
665static const char NaNU[] = "NAN";
666
667/* Write out an integerPart in hexadecimal, starting with the most
668   significant nibble.  Write out exactly COUNT hexdigits, return
669   COUNT.  */
670static unsigned int
671partAsHex (char *dst, APFloatBase::integerPart part, unsigned int count,
672           const char *hexDigitChars)
673{
674  unsigned int result = count;
675
676  assert(count != 0 && count <= APFloatBase::integerPartWidth / 4);
677
678  part >>= (APFloatBase::integerPartWidth - 4 * count);
679  while (count--) {
680    dst[count] = hexDigitChars[part & 0xf];
681    part >>= 4;
682  }
683
684  return result;
685}
686
687/* Write out an unsigned decimal integer.  */
688static char *
689writeUnsignedDecimal (char *dst, unsigned int n)
690{
691  char buff[40], *p;
692
693  p = buff;
694  do
695    *p++ = '0' + n % 10;
696  while (n /= 10);
697
698  do
699    *dst++ = *--p;
700  while (p != buff);
701
702  return dst;
703}
704
705/* Write out a signed decimal integer.  */
706static char *
707writeSignedDecimal (char *dst, int value)
708{
709  if (value < 0) {
710    *dst++ = '-';
711    dst = writeUnsignedDecimal(dst, -(unsigned) value);
712  } else
713    dst = writeUnsignedDecimal(dst, value);
714
715  return dst;
716}
717
718namespace detail {
719/* Constructors.  */
720void IEEEFloat::initialize(const fltSemantics *ourSemantics) {
721  unsigned int count;
722
723  semantics = ourSemantics;
724  count = partCount();
725  if (count > 1)
726    significand.parts = new integerPart[count];
727}
728
729void IEEEFloat::freeSignificand() {
730  if (needsCleanup())
731    delete [] significand.parts;
732}
733
734void IEEEFloat::assign(const IEEEFloat &rhs) {
735  assert(semantics == rhs.semantics);
736
737  sign = rhs.sign;
738  category = rhs.category;
739  exponent = rhs.exponent;
740  if (isFiniteNonZero() || category == fcNaN)
741    copySignificand(rhs);
742}
743
744void IEEEFloat::copySignificand(const IEEEFloat &rhs) {
745  assert(isFiniteNonZero() || category == fcNaN);
746  assert(rhs.partCount() >= partCount());
747
748  APInt::tcAssign(significandParts(), rhs.significandParts(),
749                  partCount());
750}
751
752/* Make this number a NaN, with an arbitrary but deterministic value
753   for the significand.  If double or longer, this is a signalling NaN,
754   which may not be ideal.  If float, this is QNaN(0).  */
755void IEEEFloat::makeNaN(bool SNaN, bool Negative, const APInt *fill) {
756  category = fcNaN;
757  sign = Negative;
758
759  integerPart *significand = significandParts();
760  unsigned numParts = partCount();
761
762  // Set the significand bits to the fill.
763  if (!fill || fill->getNumWords() < numParts)
764    APInt::tcSet(significand, 0, numParts);
765  if (fill) {
766    APInt::tcAssign(significand, fill->getRawData(),
767                    std::min(fill->getNumWords(), numParts));
768
769    // Zero out the excess bits of the significand.
770    unsigned bitsToPreserve = semantics->precision - 1;
771    unsigned part = bitsToPreserve / 64;
772    bitsToPreserve %= 64;
773    significand[part] &= ((1ULL << bitsToPreserve) - 1);
774    for (part++; part != numParts; ++part)
775      significand[part] = 0;
776  }
777
778  unsigned QNaNBit = semantics->precision - 2;
779
780  if (SNaN) {
781    // We always have to clear the QNaN bit to make it an SNaN.
782    APInt::tcClearBit(significand, QNaNBit);
783
784    // If there are no bits set in the payload, we have to set
785    // *something* to make it a NaN instead of an infinity;
786    // conventionally, this is the next bit down from the QNaN bit.
787    if (APInt::tcIsZero(significand, numParts))
788      APInt::tcSetBit(significand, QNaNBit - 1);
789  } else {
790    // We always have to set the QNaN bit to make it a QNaN.
791    APInt::tcSetBit(significand, QNaNBit);
792  }
793
794  // For x87 extended precision, we want to make a NaN, not a
795  // pseudo-NaN.  Maybe we should expose the ability to make
796  // pseudo-NaNs?
797  if (semantics == &semX87DoubleExtended)
798    APInt::tcSetBit(significand, QNaNBit + 1);
799}
800
801IEEEFloat &IEEEFloat::operator=(const IEEEFloat &rhs) {
802  if (this != &rhs) {
803    if (semantics != rhs.semantics) {
804      freeSignificand();
805      initialize(rhs.semantics);
806    }
807    assign(rhs);
808  }
809
810  return *this;
811}
812
813IEEEFloat &IEEEFloat::operator=(IEEEFloat &&rhs) {
814  freeSignificand();
815
816  semantics = rhs.semantics;
817  significand = rhs.significand;
818  exponent = rhs.exponent;
819  category = rhs.category;
820  sign = rhs.sign;
821
822  rhs.semantics = &semBogus;
823  return *this;
824}
825
826bool IEEEFloat::isDenormal() const {
827  return isFiniteNonZero() && (exponent == semantics->minExponent) &&
828         (APInt::tcExtractBit(significandParts(),
829                              semantics->precision - 1) == 0);
830}
831
832bool IEEEFloat::isSmallest() const {
833  // The smallest number by magnitude in our format will be the smallest
834  // denormal, i.e. the floating point number with exponent being minimum
835  // exponent and significand bitwise equal to 1 (i.e. with MSB equal to 0).
836  return isFiniteNonZero() && exponent == semantics->minExponent &&
837    significandMSB() == 0;
838}
839
840bool IEEEFloat::isSignificandAllOnes() const {
841  // Test if the significand excluding the integral bit is all ones. This allows
842  // us to test for binade boundaries.
843  const integerPart *Parts = significandParts();
844  const unsigned PartCount = partCount();
845  for (unsigned i = 0; i < PartCount - 1; i++)
846    if (~Parts[i])
847      return false;
848
849  // Set the unused high bits to all ones when we compare.
850  const unsigned NumHighBits =
851    PartCount*integerPartWidth - semantics->precision + 1;
852  assert(NumHighBits <= integerPartWidth && "Can not have more high bits to "
853         "fill than integerPartWidth");
854  const integerPart HighBitFill =
855    ~integerPart(0) << (integerPartWidth - NumHighBits);
856  if (~(Parts[PartCount - 1] | HighBitFill))
857    return false;
858
859  return true;
860}
861
862bool IEEEFloat::isSignificandAllZeros() const {
863  // Test if the significand excluding the integral bit is all zeros. This
864  // allows us to test for binade boundaries.
865  const integerPart *Parts = significandParts();
866  const unsigned PartCount = partCount();
867
868  for (unsigned i = 0; i < PartCount - 1; i++)
869    if (Parts[i])
870      return false;
871
872  const unsigned NumHighBits =
873    PartCount*integerPartWidth - semantics->precision + 1;
874  assert(NumHighBits <= integerPartWidth && "Can not have more high bits to "
875         "clear than integerPartWidth");
876  const integerPart HighBitMask = ~integerPart(0) >> NumHighBits;
877
878  if (Parts[PartCount - 1] & HighBitMask)
879    return false;
880
881  return true;
882}
883
884bool IEEEFloat::isLargest() const {
885  // The largest number by magnitude in our format will be the floating point
886  // number with maximum exponent and with significand that is all ones.
887  return isFiniteNonZero() && exponent == semantics->maxExponent
888    && isSignificandAllOnes();
889}
890
891bool IEEEFloat::isInteger() const {
892  // This could be made more efficient; I'm going for obviously correct.
893  if (!isFinite()) return false;
894  IEEEFloat truncated = *this;
895  truncated.roundToIntegral(rmTowardZero);
896  return compare(truncated) == cmpEqual;
897}
898
899bool IEEEFloat::bitwiseIsEqual(const IEEEFloat &rhs) const {
900  if (this == &rhs)
901    return true;
902  if (semantics != rhs.semantics ||
903      category != rhs.category ||
904      sign != rhs.sign)
905    return false;
906  if (category==fcZero || category==fcInfinity)
907    return true;
908
909  if (isFiniteNonZero() && exponent != rhs.exponent)
910    return false;
911
912  return std::equal(significandParts(), significandParts() + partCount(),
913                    rhs.significandParts());
914}
915
916IEEEFloat::IEEEFloat(const fltSemantics &ourSemantics, integerPart value) {
917  initialize(&ourSemantics);
918  sign = 0;
919  category = fcNormal;
920  zeroSignificand();
921  exponent = ourSemantics.precision - 1;
922  significandParts()[0] = value;
923  normalize(rmNearestTiesToEven, lfExactlyZero);
924}
925
926IEEEFloat::IEEEFloat(const fltSemantics &ourSemantics) {
927  initialize(&ourSemantics);
928  category = fcZero;
929  sign = false;
930}
931
932// Delegate to the previous constructor, because later copy constructor may
933// actually inspects category, which can't be garbage.
934IEEEFloat::IEEEFloat(const fltSemantics &ourSemantics, uninitializedTag tag)
935    : IEEEFloat(ourSemantics) {}
936
937IEEEFloat::IEEEFloat(const IEEEFloat &rhs) {
938  initialize(rhs.semantics);
939  assign(rhs);
940}
941
942IEEEFloat::IEEEFloat(IEEEFloat &&rhs) : semantics(&semBogus) {
943  *this = std::move(rhs);
944}
945
946IEEEFloat::~IEEEFloat() { freeSignificand(); }
947
948unsigned int IEEEFloat::partCount() const {
949  return partCountForBits(semantics->precision + 1);
950}
951
952const IEEEFloat::integerPart *IEEEFloat::significandParts() const {
953  return const_cast<IEEEFloat *>(this)->significandParts();
954}
955
956IEEEFloat::integerPart *IEEEFloat::significandParts() {
957  if (partCount() > 1)
958    return significand.parts;
959  else
960    return &significand.part;
961}
962
963void IEEEFloat::zeroSignificand() {
964  APInt::tcSet(significandParts(), 0, partCount());
965}
966
967/* Increment an fcNormal floating point number's significand.  */
968void IEEEFloat::incrementSignificand() {
969  integerPart carry;
970
971  carry = APInt::tcIncrement(significandParts(), partCount());
972
973  /* Our callers should never cause us to overflow.  */
974  assert(carry == 0);
975  (void)carry;
976}
977
978/* Add the significand of the RHS.  Returns the carry flag.  */
979IEEEFloat::integerPart IEEEFloat::addSignificand(const IEEEFloat &rhs) {
980  integerPart *parts;
981
982  parts = significandParts();
983
984  assert(semantics == rhs.semantics);
985  assert(exponent == rhs.exponent);
986
987  return APInt::tcAdd(parts, rhs.significandParts(), 0, partCount());
988}
989
990/* Subtract the significand of the RHS with a borrow flag.  Returns
991   the borrow flag.  */
992IEEEFloat::integerPart IEEEFloat::subtractSignificand(const IEEEFloat &rhs,
993                                                      integerPart borrow) {
994  integerPart *parts;
995
996  parts = significandParts();
997
998  assert(semantics == rhs.semantics);
999  assert(exponent == rhs.exponent);
1000
1001  return APInt::tcSubtract(parts, rhs.significandParts(), borrow,
1002                           partCount());
1003}
1004
1005/* Multiply the significand of the RHS.  If ADDEND is non-NULL, add it
1006   on to the full-precision result of the multiplication.  Returns the
1007   lost fraction.  */
1008lostFraction IEEEFloat::multiplySignificand(const IEEEFloat &rhs,
1009                                            IEEEFloat addend) {
1010  unsigned int omsb;        // One, not zero, based MSB.
1011  unsigned int partsCount, newPartsCount, precision;
1012  integerPart *lhsSignificand;
1013  integerPart scratch[4];
1014  integerPart *fullSignificand;
1015  lostFraction lost_fraction;
1016  bool ignored;
1017
1018  assert(semantics == rhs.semantics);
1019
1020  precision = semantics->precision;
1021
1022  // Allocate space for twice as many bits as the original significand, plus one
1023  // extra bit for the addition to overflow into.
1024  newPartsCount = partCountForBits(precision * 2 + 1);
1025
1026  if (newPartsCount > 4)
1027    fullSignificand = new integerPart[newPartsCount];
1028  else
1029    fullSignificand = scratch;
1030
1031  lhsSignificand = significandParts();
1032  partsCount = partCount();
1033
1034  APInt::tcFullMultiply(fullSignificand, lhsSignificand,
1035                        rhs.significandParts(), partsCount, partsCount);
1036
1037  lost_fraction = lfExactlyZero;
1038  omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
1039  exponent += rhs.exponent;
1040
1041  // Assume the operands involved in the multiplication are single-precision
1042  // FP, and the two multiplicants are:
1043  //   *this = a23 . a22 ... a0 * 2^e1
1044  //     rhs = b23 . b22 ... b0 * 2^e2
1045  // the result of multiplication is:
1046  //   *this = c48 c47 c46 . c45 ... c0 * 2^(e1+e2)
1047  // Note that there are three significant bits at the left-hand side of the
1048  // radix point: two for the multiplication, and an overflow bit for the
1049  // addition (that will always be zero at this point). Move the radix point
1050  // toward left by two bits, and adjust exponent accordingly.
1051  exponent += 2;
1052
1053  if (addend.isNonZero()) {
1054    // The intermediate result of the multiplication has "2 * precision"
1055    // signicant bit; adjust the addend to be consistent with mul result.
1056    //
1057    Significand savedSignificand = significand;
1058    const fltSemantics *savedSemantics = semantics;
1059    fltSemantics extendedSemantics;
1060    opStatus status;
1061    unsigned int extendedPrecision;
1062
1063    // Normalize our MSB to one below the top bit to allow for overflow.
1064    extendedPrecision = 2 * precision + 1;
1065    if (omsb != extendedPrecision - 1) {
1066      assert(extendedPrecision > omsb);
1067      APInt::tcShiftLeft(fullSignificand, newPartsCount,
1068                         (extendedPrecision - 1) - omsb);
1069      exponent -= (extendedPrecision - 1) - omsb;
1070    }
1071
1072    /* Create new semantics.  */
1073    extendedSemantics = *semantics;
1074    extendedSemantics.precision = extendedPrecision;
1075
1076    if (newPartsCount == 1)
1077      significand.part = fullSignificand[0];
1078    else
1079      significand.parts = fullSignificand;
1080    semantics = &extendedSemantics;
1081
1082    // Make a copy so we can convert it to the extended semantics.
1083    // Note that we cannot convert the addend directly, as the extendedSemantics
1084    // is a local variable (which we take a reference to).
1085    IEEEFloat extendedAddend(addend);
1086    status = extendedAddend.convert(extendedSemantics, rmTowardZero, &ignored);
1087    assert(status == opOK);
1088    (void)status;
1089
1090    // Shift the significand of the addend right by one bit. This guarantees
1091    // that the high bit of the significand is zero (same as fullSignificand),
1092    // so the addition will overflow (if it does overflow at all) into the top bit.
1093    lost_fraction = extendedAddend.shiftSignificandRight(1);
1094    assert(lost_fraction == lfExactlyZero &&
1095           "Lost precision while shifting addend for fused-multiply-add.");
1096
1097    lost_fraction = addOrSubtractSignificand(extendedAddend, false);
1098
1099    /* Restore our state.  */
1100    if (newPartsCount == 1)
1101      fullSignificand[0] = significand.part;
1102    significand = savedSignificand;
1103    semantics = savedSemantics;
1104
1105    omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
1106  }
1107
1108  // Convert the result having "2 * precision" significant-bits back to the one
1109  // having "precision" significant-bits. First, move the radix point from
1110  // poision "2*precision - 1" to "precision - 1". The exponent need to be
1111  // adjusted by "2*precision - 1" - "precision - 1" = "precision".
1112  exponent -= precision + 1;
1113
1114  // In case MSB resides at the left-hand side of radix point, shift the
1115  // mantissa right by some amount to make sure the MSB reside right before
1116  // the radix point (i.e. "MSB . rest-significant-bits").
1117  //
1118  // Note that the result is not normalized when "omsb < precision". So, the
1119  // caller needs to call IEEEFloat::normalize() if normalized value is
1120  // expected.
1121  if (omsb > precision) {
1122    unsigned int bits, significantParts;
1123    lostFraction lf;
1124
1125    bits = omsb - precision;
1126    significantParts = partCountForBits(omsb);
1127    lf = shiftRight(fullSignificand, significantParts, bits);
1128    lost_fraction = combineLostFractions(lf, lost_fraction);
1129    exponent += bits;
1130  }
1131
1132  APInt::tcAssign(lhsSignificand, fullSignificand, partsCount);
1133
1134  if (newPartsCount > 4)
1135    delete [] fullSignificand;
1136
1137  return lost_fraction;
1138}
1139
1140lostFraction IEEEFloat::multiplySignificand(const IEEEFloat &rhs) {
1141  return multiplySignificand(rhs, IEEEFloat(*semantics));
1142}
1143
1144/* Multiply the significands of LHS and RHS to DST.  */
1145lostFraction IEEEFloat::divideSignificand(const IEEEFloat &rhs) {
1146  unsigned int bit, i, partsCount;
1147  const integerPart *rhsSignificand;
1148  integerPart *lhsSignificand, *dividend, *divisor;
1149  integerPart scratch[4];
1150  lostFraction lost_fraction;
1151
1152  assert(semantics == rhs.semantics);
1153
1154  lhsSignificand = significandParts();
1155  rhsSignificand = rhs.significandParts();
1156  partsCount = partCount();
1157
1158  if (partsCount > 2)
1159    dividend = new integerPart[partsCount * 2];
1160  else
1161    dividend = scratch;
1162
1163  divisor = dividend + partsCount;
1164
1165  /* Copy the dividend and divisor as they will be modified in-place.  */
1166  for (i = 0; i < partsCount; i++) {
1167    dividend[i] = lhsSignificand[i];
1168    divisor[i] = rhsSignificand[i];
1169    lhsSignificand[i] = 0;
1170  }
1171
1172  exponent -= rhs.exponent;
1173
1174  unsigned int precision = semantics->precision;
1175
1176  /* Normalize the divisor.  */
1177  bit = precision - APInt::tcMSB(divisor, partsCount) - 1;
1178  if (bit) {
1179    exponent += bit;
1180    APInt::tcShiftLeft(divisor, partsCount, bit);
1181  }
1182
1183  /* Normalize the dividend.  */
1184  bit = precision - APInt::tcMSB(dividend, partsCount) - 1;
1185  if (bit) {
1186    exponent -= bit;
1187    APInt::tcShiftLeft(dividend, partsCount, bit);
1188  }
1189
1190  /* Ensure the dividend >= divisor initially for the loop below.
1191     Incidentally, this means that the division loop below is
1192     guaranteed to set the integer bit to one.  */
1193  if (APInt::tcCompare(dividend, divisor, partsCount) < 0) {
1194    exponent--;
1195    APInt::tcShiftLeft(dividend, partsCount, 1);
1196    assert(APInt::tcCompare(dividend, divisor, partsCount) >= 0);
1197  }
1198
1199  /* Long division.  */
1200  for (bit = precision; bit; bit -= 1) {
1201    if (APInt::tcCompare(dividend, divisor, partsCount) >= 0) {
1202      APInt::tcSubtract(dividend, divisor, 0, partsCount);
1203      APInt::tcSetBit(lhsSignificand, bit - 1);
1204    }
1205
1206    APInt::tcShiftLeft(dividend, partsCount, 1);
1207  }
1208
1209  /* Figure out the lost fraction.  */
1210  int cmp = APInt::tcCompare(dividend, divisor, partsCount);
1211
1212  if (cmp > 0)
1213    lost_fraction = lfMoreThanHalf;
1214  else if (cmp == 0)
1215    lost_fraction = lfExactlyHalf;
1216  else if (APInt::tcIsZero(dividend, partsCount))
1217    lost_fraction = lfExactlyZero;
1218  else
1219    lost_fraction = lfLessThanHalf;
1220
1221  if (partsCount > 2)
1222    delete [] dividend;
1223
1224  return lost_fraction;
1225}
1226
1227unsigned int IEEEFloat::significandMSB() const {
1228  return APInt::tcMSB(significandParts(), partCount());
1229}
1230
1231unsigned int IEEEFloat::significandLSB() const {
1232  return APInt::tcLSB(significandParts(), partCount());
1233}
1234
1235/* Note that a zero result is NOT normalized to fcZero.  */
1236lostFraction IEEEFloat::shiftSignificandRight(unsigned int bits) {
1237  /* Our exponent should not overflow.  */
1238  assert((ExponentType) (exponent + bits) >= exponent);
1239
1240  exponent += bits;
1241
1242  return shiftRight(significandParts(), partCount(), bits);
1243}
1244
1245/* Shift the significand left BITS bits, subtract BITS from its exponent.  */
1246void IEEEFloat::shiftSignificandLeft(unsigned int bits) {
1247  assert(bits < semantics->precision);
1248
1249  if (bits) {
1250    unsigned int partsCount = partCount();
1251
1252    APInt::tcShiftLeft(significandParts(), partsCount, bits);
1253    exponent -= bits;
1254
1255    assert(!APInt::tcIsZero(significandParts(), partsCount));
1256  }
1257}
1258
1259IEEEFloat::cmpResult
1260IEEEFloat::compareAbsoluteValue(const IEEEFloat &rhs) const {
1261  int compare;
1262
1263  assert(semantics == rhs.semantics);
1264  assert(isFiniteNonZero());
1265  assert(rhs.isFiniteNonZero());
1266
1267  compare = exponent - rhs.exponent;
1268
1269  /* If exponents are equal, do an unsigned bignum comparison of the
1270     significands.  */
1271  if (compare == 0)
1272    compare = APInt::tcCompare(significandParts(), rhs.significandParts(),
1273                               partCount());
1274
1275  if (compare > 0)
1276    return cmpGreaterThan;
1277  else if (compare < 0)
1278    return cmpLessThan;
1279  else
1280    return cmpEqual;
1281}
1282
1283/* Handle overflow.  Sign is preserved.  We either become infinity or
1284   the largest finite number.  */
1285IEEEFloat::opStatus IEEEFloat::handleOverflow(roundingMode rounding_mode) {
1286  /* Infinity?  */
1287  if (rounding_mode == rmNearestTiesToEven ||
1288      rounding_mode == rmNearestTiesToAway ||
1289      (rounding_mode == rmTowardPositive && !sign) ||
1290      (rounding_mode == rmTowardNegative && sign)) {
1291    category = fcInfinity;
1292    return (opStatus) (opOverflow | opInexact);
1293  }
1294
1295  /* Otherwise we become the largest finite number.  */
1296  category = fcNormal;
1297  exponent = semantics->maxExponent;
1298  APInt::tcSetLeastSignificantBits(significandParts(), partCount(),
1299                                   semantics->precision);
1300
1301  return opInexact;
1302}
1303
1304/* Returns TRUE if, when truncating the current number, with BIT the
1305   new LSB, with the given lost fraction and rounding mode, the result
1306   would need to be rounded away from zero (i.e., by increasing the
1307   signficand).  This routine must work for fcZero of both signs, and
1308   fcNormal numbers.  */
1309bool IEEEFloat::roundAwayFromZero(roundingMode rounding_mode,
1310                                  lostFraction lost_fraction,
1311                                  unsigned int bit) const {
1312  /* NaNs and infinities should not have lost fractions.  */
1313  assert(isFiniteNonZero() || category == fcZero);
1314
1315  /* Current callers never pass this so we don't handle it.  */
1316  assert(lost_fraction != lfExactlyZero);
1317
1318  switch (rounding_mode) {
1319  case rmNearestTiesToAway:
1320    return lost_fraction == lfExactlyHalf || lost_fraction == lfMoreThanHalf;
1321
1322  case rmNearestTiesToEven:
1323    if (lost_fraction == lfMoreThanHalf)
1324      return true;
1325
1326    /* Our zeroes don't have a significand to test.  */
1327    if (lost_fraction == lfExactlyHalf && category != fcZero)
1328      return APInt::tcExtractBit(significandParts(), bit);
1329
1330    return false;
1331
1332  case rmTowardZero:
1333    return false;
1334
1335  case rmTowardPositive:
1336    return !sign;
1337
1338  case rmTowardNegative:
1339    return sign;
1340
1341  default:
1342    break;
1343  }
1344  llvm_unreachable("Invalid rounding mode found");
1345}
1346
1347IEEEFloat::opStatus IEEEFloat::normalize(roundingMode rounding_mode,
1348                                         lostFraction lost_fraction) {
1349  unsigned int omsb;                /* One, not zero, based MSB.  */
1350  int exponentChange;
1351
1352  if (!isFiniteNonZero())
1353    return opOK;
1354
1355  /* Before rounding normalize the exponent of fcNormal numbers.  */
1356  omsb = significandMSB() + 1;
1357
1358  if (omsb) {
1359    /* OMSB is numbered from 1.  We want to place it in the integer
1360       bit numbered PRECISION if possible, with a compensating change in
1361       the exponent.  */
1362    exponentChange = omsb - semantics->precision;
1363
1364    /* If the resulting exponent is too high, overflow according to
1365       the rounding mode.  */
1366    if (exponent + exponentChange > semantics->maxExponent)
1367      return handleOverflow(rounding_mode);
1368
1369    /* Subnormal numbers have exponent minExponent, and their MSB
1370       is forced based on that.  */
1371    if (exponent + exponentChange < semantics->minExponent)
1372      exponentChange = semantics->minExponent - exponent;
1373
1374    /* Shifting left is easy as we don't lose precision.  */
1375    if (exponentChange < 0) {
1376      assert(lost_fraction == lfExactlyZero);
1377
1378      shiftSignificandLeft(-exponentChange);
1379
1380      return opOK;
1381    }
1382
1383    if (exponentChange > 0) {
1384      lostFraction lf;
1385
1386      /* Shift right and capture any new lost fraction.  */
1387      lf = shiftSignificandRight(exponentChange);
1388
1389      lost_fraction = combineLostFractions(lf, lost_fraction);
1390
1391      /* Keep OMSB up-to-date.  */
1392      if (omsb > (unsigned) exponentChange)
1393        omsb -= exponentChange;
1394      else
1395        omsb = 0;
1396    }
1397  }
1398
1399  /* Now round the number according to rounding_mode given the lost
1400     fraction.  */
1401
1402  /* As specified in IEEE 754, since we do not trap we do not report
1403     underflow for exact results.  */
1404  if (lost_fraction == lfExactlyZero) {
1405    /* Canonicalize zeroes.  */
1406    if (omsb == 0)
1407      category = fcZero;
1408
1409    return opOK;
1410  }
1411
1412  /* Increment the significand if we're rounding away from zero.  */
1413  if (roundAwayFromZero(rounding_mode, lost_fraction, 0)) {
1414    if (omsb == 0)
1415      exponent = semantics->minExponent;
1416
1417    incrementSignificand();
1418    omsb = significandMSB() + 1;
1419
1420    /* Did the significand increment overflow?  */
1421    if (omsb == (unsigned) semantics->precision + 1) {
1422      /* Renormalize by incrementing the exponent and shifting our
1423         significand right one.  However if we already have the
1424         maximum exponent we overflow to infinity.  */
1425      if (exponent == semantics->maxExponent) {
1426        category = fcInfinity;
1427
1428        return (opStatus) (opOverflow | opInexact);
1429      }
1430
1431      shiftSignificandRight(1);
1432
1433      return opInexact;
1434    }
1435  }
1436
1437  /* The normal case - we were and are not denormal, and any
1438     significand increment above didn't overflow.  */
1439  if (omsb == semantics->precision)
1440    return opInexact;
1441
1442  /* We have a non-zero denormal.  */
1443  assert(omsb < semantics->precision);
1444
1445  /* Canonicalize zeroes.  */
1446  if (omsb == 0)
1447    category = fcZero;
1448
1449  /* The fcZero case is a denormal that underflowed to zero.  */
1450  return (opStatus) (opUnderflow | opInexact);
1451}
1452
1453IEEEFloat::opStatus IEEEFloat::addOrSubtractSpecials(const IEEEFloat &rhs,
1454                                                     bool subtract) {
1455  switch (PackCategoriesIntoKey(category, rhs.category)) {
1456  default:
1457    llvm_unreachable(nullptr);
1458
1459  case PackCategoriesIntoKey(fcZero, fcNaN):
1460  case PackCategoriesIntoKey(fcNormal, fcNaN):
1461  case PackCategoriesIntoKey(fcInfinity, fcNaN):
1462    assign(rhs);
1463    LLVM_FALLTHROUGH;
1464  case PackCategoriesIntoKey(fcNaN, fcZero):
1465  case PackCategoriesIntoKey(fcNaN, fcNormal):
1466  case PackCategoriesIntoKey(fcNaN, fcInfinity):
1467  case PackCategoriesIntoKey(fcNaN, fcNaN):
1468    if (isSignaling()) {
1469      makeQuiet();
1470      return opInvalidOp;
1471    }
1472    return rhs.isSignaling() ? opInvalidOp : opOK;
1473
1474  case PackCategoriesIntoKey(fcNormal, fcZero):
1475  case PackCategoriesIntoKey(fcInfinity, fcNormal):
1476  case PackCategoriesIntoKey(fcInfinity, fcZero):
1477    return opOK;
1478
1479  case PackCategoriesIntoKey(fcNormal, fcInfinity):
1480  case PackCategoriesIntoKey(fcZero, fcInfinity):
1481    category = fcInfinity;
1482    sign = rhs.sign ^ subtract;
1483    return opOK;
1484
1485  case PackCategoriesIntoKey(fcZero, fcNormal):
1486    assign(rhs);
1487    sign = rhs.sign ^ subtract;
1488    return opOK;
1489
1490  case PackCategoriesIntoKey(fcZero, fcZero):
1491    /* Sign depends on rounding mode; handled by caller.  */
1492    return opOK;
1493
1494  case PackCategoriesIntoKey(fcInfinity, fcInfinity):
1495    /* Differently signed infinities can only be validly
1496       subtracted.  */
1497    if (((sign ^ rhs.sign)!=0) != subtract) {
1498      makeNaN();
1499      return opInvalidOp;
1500    }
1501
1502    return opOK;
1503
1504  case PackCategoriesIntoKey(fcNormal, fcNormal):
1505    return opDivByZero;
1506  }
1507}
1508
1509/* Add or subtract two normal numbers.  */
1510lostFraction IEEEFloat::addOrSubtractSignificand(const IEEEFloat &rhs,
1511                                                 bool subtract) {
1512  integerPart carry;
1513  lostFraction lost_fraction;
1514  int bits;
1515
1516  /* Determine if the operation on the absolute values is effectively
1517     an addition or subtraction.  */
1518  subtract ^= static_cast<bool>(sign ^ rhs.sign);
1519
1520  /* Are we bigger exponent-wise than the RHS?  */
1521  bits = exponent - rhs.exponent;
1522
1523  /* Subtraction is more subtle than one might naively expect.  */
1524  if (subtract) {
1525    IEEEFloat temp_rhs(rhs);
1526
1527    if (bits == 0)
1528      lost_fraction = lfExactlyZero;
1529    else if (bits > 0) {
1530      lost_fraction = temp_rhs.shiftSignificandRight(bits - 1);
1531      shiftSignificandLeft(1);
1532    } else {
1533      lost_fraction = shiftSignificandRight(-bits - 1);
1534      temp_rhs.shiftSignificandLeft(1);
1535    }
1536
1537    // Should we reverse the subtraction.
1538    if (compareAbsoluteValue(temp_rhs) == cmpLessThan) {
1539      carry = temp_rhs.subtractSignificand
1540        (*this, lost_fraction != lfExactlyZero);
1541      copySignificand(temp_rhs);
1542      sign = !sign;
1543    } else {
1544      carry = subtractSignificand
1545        (temp_rhs, lost_fraction != lfExactlyZero);
1546    }
1547
1548    /* Invert the lost fraction - it was on the RHS and
1549       subtracted.  */
1550    if (lost_fraction == lfLessThanHalf)
1551      lost_fraction = lfMoreThanHalf;
1552    else if (lost_fraction == lfMoreThanHalf)
1553      lost_fraction = lfLessThanHalf;
1554
1555    /* The code above is intended to ensure that no borrow is
1556       necessary.  */
1557    assert(!carry);
1558    (void)carry;
1559  } else {
1560    if (bits > 0) {
1561      IEEEFloat temp_rhs(rhs);
1562
1563      lost_fraction = temp_rhs.shiftSignificandRight(bits);
1564      carry = addSignificand(temp_rhs);
1565    } else {
1566      lost_fraction = shiftSignificandRight(-bits);
1567      carry = addSignificand(rhs);
1568    }
1569
1570    /* We have a guard bit; generating a carry cannot happen.  */
1571    assert(!carry);
1572    (void)carry;
1573  }
1574
1575  return lost_fraction;
1576}
1577
1578IEEEFloat::opStatus IEEEFloat::multiplySpecials(const IEEEFloat &rhs) {
1579  switch (PackCategoriesIntoKey(category, rhs.category)) {
1580  default:
1581    llvm_unreachable(nullptr);
1582
1583  case PackCategoriesIntoKey(fcZero, fcNaN):
1584  case PackCategoriesIntoKey(fcNormal, fcNaN):
1585  case PackCategoriesIntoKey(fcInfinity, fcNaN):
1586    assign(rhs);
1587    sign = false;
1588    LLVM_FALLTHROUGH;
1589  case PackCategoriesIntoKey(fcNaN, fcZero):
1590  case PackCategoriesIntoKey(fcNaN, fcNormal):
1591  case PackCategoriesIntoKey(fcNaN, fcInfinity):
1592  case PackCategoriesIntoKey(fcNaN, fcNaN):
1593    sign ^= rhs.sign; // restore the original sign
1594    if (isSignaling()) {
1595      makeQuiet();
1596      return opInvalidOp;
1597    }
1598    return rhs.isSignaling() ? opInvalidOp : opOK;
1599
1600  case PackCategoriesIntoKey(fcNormal, fcInfinity):
1601  case PackCategoriesIntoKey(fcInfinity, fcNormal):
1602  case PackCategoriesIntoKey(fcInfinity, fcInfinity):
1603    category = fcInfinity;
1604    return opOK;
1605
1606  case PackCategoriesIntoKey(fcZero, fcNormal):
1607  case PackCategoriesIntoKey(fcNormal, fcZero):
1608  case PackCategoriesIntoKey(fcZero, fcZero):
1609    category = fcZero;
1610    return opOK;
1611
1612  case PackCategoriesIntoKey(fcZero, fcInfinity):
1613  case PackCategoriesIntoKey(fcInfinity, fcZero):
1614    makeNaN();
1615    return opInvalidOp;
1616
1617  case PackCategoriesIntoKey(fcNormal, fcNormal):
1618    return opOK;
1619  }
1620}
1621
1622IEEEFloat::opStatus IEEEFloat::divideSpecials(const IEEEFloat &rhs) {
1623  switch (PackCategoriesIntoKey(category, rhs.category)) {
1624  default:
1625    llvm_unreachable(nullptr);
1626
1627  case PackCategoriesIntoKey(fcZero, fcNaN):
1628  case PackCategoriesIntoKey(fcNormal, fcNaN):
1629  case PackCategoriesIntoKey(fcInfinity, fcNaN):
1630    assign(rhs);
1631    sign = false;
1632    LLVM_FALLTHROUGH;
1633  case PackCategoriesIntoKey(fcNaN, fcZero):
1634  case PackCategoriesIntoKey(fcNaN, fcNormal):
1635  case PackCategoriesIntoKey(fcNaN, fcInfinity):
1636  case PackCategoriesIntoKey(fcNaN, fcNaN):
1637    sign ^= rhs.sign; // restore the original sign
1638    if (isSignaling()) {
1639      makeQuiet();
1640      return opInvalidOp;
1641    }
1642    return rhs.isSignaling() ? opInvalidOp : opOK;
1643
1644  case PackCategoriesIntoKey(fcInfinity, fcZero):
1645  case PackCategoriesIntoKey(fcInfinity, fcNormal):
1646  case PackCategoriesIntoKey(fcZero, fcInfinity):
1647  case PackCategoriesIntoKey(fcZero, fcNormal):
1648    return opOK;
1649
1650  case PackCategoriesIntoKey(fcNormal, fcInfinity):
1651    category = fcZero;
1652    return opOK;
1653
1654  case PackCategoriesIntoKey(fcNormal, fcZero):
1655    category = fcInfinity;
1656    return opDivByZero;
1657
1658  case PackCategoriesIntoKey(fcInfinity, fcInfinity):
1659  case PackCategoriesIntoKey(fcZero, fcZero):
1660    makeNaN();
1661    return opInvalidOp;
1662
1663  case PackCategoriesIntoKey(fcNormal, fcNormal):
1664    return opOK;
1665  }
1666}
1667
1668IEEEFloat::opStatus IEEEFloat::modSpecials(const IEEEFloat &rhs) {
1669  switch (PackCategoriesIntoKey(category, rhs.category)) {
1670  default:
1671    llvm_unreachable(nullptr);
1672
1673  case PackCategoriesIntoKey(fcZero, fcNaN):
1674  case PackCategoriesIntoKey(fcNormal, fcNaN):
1675  case PackCategoriesIntoKey(fcInfinity, fcNaN):
1676    assign(rhs);
1677    LLVM_FALLTHROUGH;
1678  case PackCategoriesIntoKey(fcNaN, fcZero):
1679  case PackCategoriesIntoKey(fcNaN, fcNormal):
1680  case PackCategoriesIntoKey(fcNaN, fcInfinity):
1681  case PackCategoriesIntoKey(fcNaN, fcNaN):
1682    if (isSignaling()) {
1683      makeQuiet();
1684      return opInvalidOp;
1685    }
1686    return rhs.isSignaling() ? opInvalidOp : opOK;
1687
1688  case PackCategoriesIntoKey(fcZero, fcInfinity):
1689  case PackCategoriesIntoKey(fcZero, fcNormal):
1690  case PackCategoriesIntoKey(fcNormal, fcInfinity):
1691    return opOK;
1692
1693  case PackCategoriesIntoKey(fcNormal, fcZero):
1694  case PackCategoriesIntoKey(fcInfinity, fcZero):
1695  case PackCategoriesIntoKey(fcInfinity, fcNormal):
1696  case PackCategoriesIntoKey(fcInfinity, fcInfinity):
1697  case PackCategoriesIntoKey(fcZero, fcZero):
1698    makeNaN();
1699    return opInvalidOp;
1700
1701  case PackCategoriesIntoKey(fcNormal, fcNormal):
1702    return opOK;
1703  }
1704}
1705
1706IEEEFloat::opStatus IEEEFloat::remainderSpecials(const IEEEFloat &rhs) {
1707  switch (PackCategoriesIntoKey(category, rhs.category)) {
1708  default:
1709    llvm_unreachable(nullptr);
1710
1711  case PackCategoriesIntoKey(fcZero, fcNaN):
1712  case PackCategoriesIntoKey(fcNormal, fcNaN):
1713  case PackCategoriesIntoKey(fcInfinity, fcNaN):
1714    assign(rhs);
1715    LLVM_FALLTHROUGH;
1716  case PackCategoriesIntoKey(fcNaN, fcZero):
1717  case PackCategoriesIntoKey(fcNaN, fcNormal):
1718  case PackCategoriesIntoKey(fcNaN, fcInfinity):
1719  case PackCategoriesIntoKey(fcNaN, fcNaN):
1720    if (isSignaling()) {
1721      makeQuiet();
1722      return opInvalidOp;
1723    }
1724    return rhs.isSignaling() ? opInvalidOp : opOK;
1725
1726  case PackCategoriesIntoKey(fcZero, fcInfinity):
1727  case PackCategoriesIntoKey(fcZero, fcNormal):
1728  case PackCategoriesIntoKey(fcNormal, fcInfinity):
1729    return opOK;
1730
1731  case PackCategoriesIntoKey(fcNormal, fcZero):
1732  case PackCategoriesIntoKey(fcInfinity, fcZero):
1733  case PackCategoriesIntoKey(fcInfinity, fcNormal):
1734  case PackCategoriesIntoKey(fcInfinity, fcInfinity):
1735  case PackCategoriesIntoKey(fcZero, fcZero):
1736    makeNaN();
1737    return opInvalidOp;
1738
1739  case PackCategoriesIntoKey(fcNormal, fcNormal):
1740    return opDivByZero; // fake status, indicating this is not a special case
1741  }
1742}
1743
1744/* Change sign.  */
1745void IEEEFloat::changeSign() {
1746  /* Look mummy, this one's easy.  */
1747  sign = !sign;
1748}
1749
1750/* Normalized addition or subtraction.  */
1751IEEEFloat::opStatus IEEEFloat::addOrSubtract(const IEEEFloat &rhs,
1752                                             roundingMode rounding_mode,
1753                                             bool subtract) {
1754  opStatus fs;
1755
1756  fs = addOrSubtractSpecials(rhs, subtract);
1757
1758  /* This return code means it was not a simple case.  */
1759  if (fs == opDivByZero) {
1760    lostFraction lost_fraction;
1761
1762    lost_fraction = addOrSubtractSignificand(rhs, subtract);
1763    fs = normalize(rounding_mode, lost_fraction);
1764
1765    /* Can only be zero if we lost no fraction.  */
1766    assert(category != fcZero || lost_fraction == lfExactlyZero);
1767  }
1768
1769  /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1770     positive zero unless rounding to minus infinity, except that
1771     adding two like-signed zeroes gives that zero.  */
1772  if (category == fcZero) {
1773    if (rhs.category != fcZero || (sign == rhs.sign) == subtract)
1774      sign = (rounding_mode == rmTowardNegative);
1775  }
1776
1777  return fs;
1778}
1779
1780/* Normalized addition.  */
1781IEEEFloat::opStatus IEEEFloat::add(const IEEEFloat &rhs,
1782                                   roundingMode rounding_mode) {
1783  return addOrSubtract(rhs, rounding_mode, false);
1784}
1785
1786/* Normalized subtraction.  */
1787IEEEFloat::opStatus IEEEFloat::subtract(const IEEEFloat &rhs,
1788                                        roundingMode rounding_mode) {
1789  return addOrSubtract(rhs, rounding_mode, true);
1790}
1791
1792/* Normalized multiply.  */
1793IEEEFloat::opStatus IEEEFloat::multiply(const IEEEFloat &rhs,
1794                                        roundingMode rounding_mode) {
1795  opStatus fs;
1796
1797  sign ^= rhs.sign;
1798  fs = multiplySpecials(rhs);
1799
1800  if (isFiniteNonZero()) {
1801    lostFraction lost_fraction = multiplySignificand(rhs);
1802    fs = normalize(rounding_mode, lost_fraction);
1803    if (lost_fraction != lfExactlyZero)
1804      fs = (opStatus) (fs | opInexact);
1805  }
1806
1807  return fs;
1808}
1809
1810/* Normalized divide.  */
1811IEEEFloat::opStatus IEEEFloat::divide(const IEEEFloat &rhs,
1812                                      roundingMode rounding_mode) {
1813  opStatus fs;
1814
1815  sign ^= rhs.sign;
1816  fs = divideSpecials(rhs);
1817
1818  if (isFiniteNonZero()) {
1819    lostFraction lost_fraction = divideSignificand(rhs);
1820    fs = normalize(rounding_mode, lost_fraction);
1821    if (lost_fraction != lfExactlyZero)
1822      fs = (opStatus) (fs | opInexact);
1823  }
1824
1825  return fs;
1826}
1827
1828/* Normalized remainder.  */
1829IEEEFloat::opStatus IEEEFloat::remainder(const IEEEFloat &rhs) {
1830  opStatus fs;
1831  unsigned int origSign = sign;
1832
1833  // First handle the special cases.
1834  fs = remainderSpecials(rhs);
1835  if (fs != opDivByZero)
1836    return fs;
1837
1838  fs = opOK;
1839
1840  // Make sure the current value is less than twice the denom. If the addition
1841  // did not succeed (an overflow has happened), which means that the finite
1842  // value we currently posses must be less than twice the denom (as we are
1843  // using the same semantics).
1844  IEEEFloat P2 = rhs;
1845  if (P2.add(rhs, rmNearestTiesToEven) == opOK) {
1846    fs = mod(P2);
1847    assert(fs == opOK);
1848  }
1849
1850  // Lets work with absolute numbers.
1851  IEEEFloat P = rhs;
1852  P.sign = false;
1853  sign = false;
1854
1855  //
1856  // To calculate the remainder we use the following scheme.
1857  //
1858  // The remainder is defained as follows:
1859  //
1860  // remainder = numer - rquot * denom = x - r * p
1861  //
1862  // Where r is the result of: x/p, rounded toward the nearest integral value
1863  // (with halfway cases rounded toward the even number).
1864  //
1865  // Currently, (after x mod 2p):
1866  // r is the number of 2p's present inside x, which is inherently, an even
1867  // number of p's.
1868  //
1869  // We may split the remaining calculation into 4 options:
1870  // - if x < 0.5p then we round to the nearest number with is 0, and are done.
1871  // - if x == 0.5p then we round to the nearest even number which is 0, and we
1872  //   are done as well.
1873  // - if 0.5p < x < p then we round to nearest number which is 1, and we have
1874  //   to subtract 1p at least once.
1875  // - if x >= p then we must subtract p at least once, as x must be a
1876  //   remainder.
1877  //
1878  // By now, we were done, or we added 1 to r, which in turn, now an odd number.
1879  //
1880  // We can now split the remaining calculation to the following 3 options:
1881  // - if x < 0.5p then we round to the nearest number with is 0, and are done.
1882  // - if x == 0.5p then we round to the nearest even number. As r is odd, we
1883  //   must round up to the next even number. so we must subtract p once more.
1884  // - if x > 0.5p (and inherently x < p) then we must round r up to the next
1885  //   integral, and subtract p once more.
1886  //
1887
1888  // Extend the semantics to prevent an overflow/underflow or inexact result.
1889  bool losesInfo;
1890  fltSemantics extendedSemantics = *semantics;
1891  extendedSemantics.maxExponent++;
1892  extendedSemantics.minExponent--;
1893  extendedSemantics.precision += 2;
1894
1895  IEEEFloat VEx = *this;
1896  fs = VEx.convert(extendedSemantics, rmNearestTiesToEven, &losesInfo);
1897  assert(fs == opOK && !losesInfo);
1898  IEEEFloat PEx = P;
1899  fs = PEx.convert(extendedSemantics, rmNearestTiesToEven, &losesInfo);
1900  assert(fs == opOK && !losesInfo);
1901
1902  // It is simpler to work with 2x instead of 0.5p, and we do not need to lose
1903  // any fraction.
1904  fs = VEx.add(VEx, rmNearestTiesToEven);
1905  assert(fs == opOK);
1906
1907  if (VEx.compare(PEx) == cmpGreaterThan) {
1908    fs = subtract(P, rmNearestTiesToEven);
1909    assert(fs == opOK);
1910
1911    // Make VEx = this.add(this), but because we have different semantics, we do
1912    // not want to `convert` again, so we just subtract PEx twice (which equals
1913    // to the desired value).
1914    fs = VEx.subtract(PEx, rmNearestTiesToEven);
1915    assert(fs == opOK);
1916    fs = VEx.subtract(PEx, rmNearestTiesToEven);
1917    assert(fs == opOK);
1918
1919    cmpResult result = VEx.compare(PEx);
1920    if (result == cmpGreaterThan || result == cmpEqual) {
1921      fs = subtract(P, rmNearestTiesToEven);
1922      assert(fs == opOK);
1923    }
1924  }
1925
1926  if (isZero())
1927    sign = origSign;    // IEEE754 requires this
1928  else
1929    sign ^= origSign;
1930  return fs;
1931}
1932
1933/* Normalized llvm frem (C fmod). */
1934IEEEFloat::opStatus IEEEFloat::mod(const IEEEFloat &rhs) {
1935  opStatus fs;
1936  fs = modSpecials(rhs);
1937  unsigned int origSign = sign;
1938
1939  while (isFiniteNonZero() && rhs.isFiniteNonZero() &&
1940         compareAbsoluteValue(rhs) != cmpLessThan) {
1941    IEEEFloat V = scalbn(rhs, ilogb(*this) - ilogb(rhs), rmNearestTiesToEven);
1942    if (compareAbsoluteValue(V) == cmpLessThan)
1943      V = scalbn(V, -1, rmNearestTiesToEven);
1944    V.sign = sign;
1945
1946    fs = subtract(V, rmNearestTiesToEven);
1947    assert(fs==opOK);
1948  }
1949  if (isZero())
1950    sign = origSign; // fmod requires this
1951  return fs;
1952}
1953
1954/* Normalized fused-multiply-add.  */
1955IEEEFloat::opStatus IEEEFloat::fusedMultiplyAdd(const IEEEFloat &multiplicand,
1956                                                const IEEEFloat &addend,
1957                                                roundingMode rounding_mode) {
1958  opStatus fs;
1959
1960  /* Post-multiplication sign, before addition.  */
1961  sign ^= multiplicand.sign;
1962
1963  /* If and only if all arguments are normal do we need to do an
1964     extended-precision calculation.  */
1965  if (isFiniteNonZero() &&
1966      multiplicand.isFiniteNonZero() &&
1967      addend.isFinite()) {
1968    lostFraction lost_fraction;
1969
1970    lost_fraction = multiplySignificand(multiplicand, addend);
1971    fs = normalize(rounding_mode, lost_fraction);
1972    if (lost_fraction != lfExactlyZero)
1973      fs = (opStatus) (fs | opInexact);
1974
1975    /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1976       positive zero unless rounding to minus infinity, except that
1977       adding two like-signed zeroes gives that zero.  */
1978    if (category == fcZero && !(fs & opUnderflow) && sign != addend.sign)
1979      sign = (rounding_mode == rmTowardNegative);
1980  } else {
1981    fs = multiplySpecials(multiplicand);
1982
1983    /* FS can only be opOK or opInvalidOp.  There is no more work
1984       to do in the latter case.  The IEEE-754R standard says it is
1985       implementation-defined in this case whether, if ADDEND is a
1986       quiet NaN, we raise invalid op; this implementation does so.
1987
1988       If we need to do the addition we can do so with normal
1989       precision.  */
1990    if (fs == opOK)
1991      fs = addOrSubtract(addend, rounding_mode, false);
1992  }
1993
1994  return fs;
1995}
1996
1997/* Rounding-mode correct round to integral value.  */
1998IEEEFloat::opStatus IEEEFloat::roundToIntegral(roundingMode rounding_mode) {
1999  opStatus fs;
2000
2001  if (isInfinity())
2002    // [IEEE Std 754-2008 6.1]:
2003    // The behavior of infinity in floating-point arithmetic is derived from the
2004    // limiting cases of real arithmetic with operands of arbitrarily
2005    // large magnitude, when such a limit exists.
2006    // ...
2007    // Operations on infinite operands are usually exact and therefore signal no
2008    // exceptions ...
2009    return opOK;
2010
2011  if (isNaN()) {
2012    if (isSignaling()) {
2013      // [IEEE Std 754-2008 6.2]:
2014      // Under default exception handling, any operation signaling an invalid
2015      // operation exception and for which a floating-point result is to be
2016      // delivered shall deliver a quiet NaN.
2017      makeQuiet();
2018      // [IEEE Std 754-2008 6.2]:
2019      // Signaling NaNs shall be reserved operands that, under default exception
2020      // handling, signal the invalid operation exception(see 7.2) for every
2021      // general-computational and signaling-computational operation except for
2022      // the conversions described in 5.12.
2023      return opInvalidOp;
2024    } else {
2025      // [IEEE Std 754-2008 6.2]:
2026      // For an operation with quiet NaN inputs, other than maximum and minimum
2027      // operations, if a floating-point result is to be delivered the result
2028      // shall be a quiet NaN which should be one of the input NaNs.
2029      // ...
2030      // Every general-computational and quiet-computational operation involving
2031      // one or more input NaNs, none of them signaling, shall signal no
2032      // exception, except fusedMultiplyAdd might signal the invalid operation
2033      // exception(see 7.2).
2034      return opOK;
2035    }
2036  }
2037
2038  if (isZero()) {
2039    // [IEEE Std 754-2008 6.3]:
2040    // ... the sign of the result of conversions, the quantize operation, the
2041    // roundToIntegral operations, and the roundToIntegralExact(see 5.3.1) is
2042    // the sign of the first or only operand.
2043    return opOK;
2044  }
2045
2046  // If the exponent is large enough, we know that this value is already
2047  // integral, and the arithmetic below would potentially cause it to saturate
2048  // to +/-Inf.  Bail out early instead.
2049  if (exponent+1 >= (int)semanticsPrecision(*semantics))
2050    return opOK;
2051
2052  // The algorithm here is quite simple: we add 2^(p-1), where p is the
2053  // precision of our format, and then subtract it back off again.  The choice
2054  // of rounding modes for the addition/subtraction determines the rounding mode
2055  // for our integral rounding as well.
2056  // NOTE: When the input value is negative, we do subtraction followed by
2057  // addition instead.
2058  APInt IntegerConstant(NextPowerOf2(semanticsPrecision(*semantics)), 1);
2059  IntegerConstant <<= semanticsPrecision(*semantics)-1;
2060  IEEEFloat MagicConstant(*semantics);
2061  fs = MagicConstant.convertFromAPInt(IntegerConstant, false,
2062                                      rmNearestTiesToEven);
2063  assert(fs == opOK);
2064  MagicConstant.sign = sign;
2065
2066  // Preserve the input sign so that we can handle the case of zero result
2067  // correctly.
2068  bool inputSign = isNegative();
2069
2070  fs = add(MagicConstant, rounding_mode);
2071
2072  // Current value and 'MagicConstant' are both integers, so the result of the
2073  // subtraction is always exact according to Sterbenz' lemma.
2074  subtract(MagicConstant, rounding_mode);
2075
2076  // Restore the input sign.
2077  if (inputSign != isNegative())
2078    changeSign();
2079
2080  return fs;
2081}
2082
2083
2084/* Comparison requires normalized numbers.  */
2085IEEEFloat::cmpResult IEEEFloat::compare(const IEEEFloat &rhs) const {
2086  cmpResult result;
2087
2088  assert(semantics == rhs.semantics);
2089
2090  switch (PackCategoriesIntoKey(category, rhs.category)) {
2091  default:
2092    llvm_unreachable(nullptr);
2093
2094  case PackCategoriesIntoKey(fcNaN, fcZero):
2095  case PackCategoriesIntoKey(fcNaN, fcNormal):
2096  case PackCategoriesIntoKey(fcNaN, fcInfinity):
2097  case PackCategoriesIntoKey(fcNaN, fcNaN):
2098  case PackCategoriesIntoKey(fcZero, fcNaN):
2099  case PackCategoriesIntoKey(fcNormal, fcNaN):
2100  case PackCategoriesIntoKey(fcInfinity, fcNaN):
2101    return cmpUnordered;
2102
2103  case PackCategoriesIntoKey(fcInfinity, fcNormal):
2104  case PackCategoriesIntoKey(fcInfinity, fcZero):
2105  case PackCategoriesIntoKey(fcNormal, fcZero):
2106    if (sign)
2107      return cmpLessThan;
2108    else
2109      return cmpGreaterThan;
2110
2111  case PackCategoriesIntoKey(fcNormal, fcInfinity):
2112  case PackCategoriesIntoKey(fcZero, fcInfinity):
2113  case PackCategoriesIntoKey(fcZero, fcNormal):
2114    if (rhs.sign)
2115      return cmpGreaterThan;
2116    else
2117      return cmpLessThan;
2118
2119  case PackCategoriesIntoKey(fcInfinity, fcInfinity):
2120    if (sign == rhs.sign)
2121      return cmpEqual;
2122    else if (sign)
2123      return cmpLessThan;
2124    else
2125      return cmpGreaterThan;
2126
2127  case PackCategoriesIntoKey(fcZero, fcZero):
2128    return cmpEqual;
2129
2130  case PackCategoriesIntoKey(fcNormal, fcNormal):
2131    break;
2132  }
2133
2134  /* Two normal numbers.  Do they have the same sign?  */
2135  if (sign != rhs.sign) {
2136    if (sign)
2137      result = cmpLessThan;
2138    else
2139      result = cmpGreaterThan;
2140  } else {
2141    /* Compare absolute values; invert result if negative.  */
2142    result = compareAbsoluteValue(rhs);
2143
2144    if (sign) {
2145      if (result == cmpLessThan)
2146        result = cmpGreaterThan;
2147      else if (result == cmpGreaterThan)
2148        result = cmpLessThan;
2149    }
2150  }
2151
2152  return result;
2153}
2154
2155/// IEEEFloat::convert - convert a value of one floating point type to another.
2156/// The return value corresponds to the IEEE754 exceptions.  *losesInfo
2157/// records whether the transformation lost information, i.e. whether
2158/// converting the result back to the original type will produce the
2159/// original value (this is almost the same as return value==fsOK, but there
2160/// are edge cases where this is not so).
2161
2162IEEEFloat::opStatus IEEEFloat::convert(const fltSemantics &toSemantics,
2163                                       roundingMode rounding_mode,
2164                                       bool *losesInfo) {
2165  lostFraction lostFraction;
2166  unsigned int newPartCount, oldPartCount;
2167  opStatus fs;
2168  int shift;
2169  const fltSemantics &fromSemantics = *semantics;
2170
2171  lostFraction = lfExactlyZero;
2172  newPartCount = partCountForBits(toSemantics.precision + 1);
2173  oldPartCount = partCount();
2174  shift = toSemantics.precision - fromSemantics.precision;
2175
2176  bool X86SpecialNan = false;
2177  if (&fromSemantics == &semX87DoubleExtended &&
2178      &toSemantics != &semX87DoubleExtended && category == fcNaN &&
2179      (!(*significandParts() & 0x8000000000000000ULL) ||
2180       !(*significandParts() & 0x4000000000000000ULL))) {
2181    // x86 has some unusual NaNs which cannot be represented in any other
2182    // format; note them here.
2183    X86SpecialNan = true;
2184  }
2185
2186  // If this is a truncation of a denormal number, and the target semantics
2187  // has larger exponent range than the source semantics (this can happen
2188  // when truncating from PowerPC double-double to double format), the
2189  // right shift could lose result mantissa bits.  Adjust exponent instead
2190  // of performing excessive shift.
2191  if (shift < 0 && isFiniteNonZero()) {
2192    int exponentChange = significandMSB() + 1 - fromSemantics.precision;
2193    if (exponent + exponentChange < toSemantics.minExponent)
2194      exponentChange = toSemantics.minExponent - exponent;
2195    if (exponentChange < shift)
2196      exponentChange = shift;
2197    if (exponentChange < 0) {
2198      shift -= exponentChange;
2199      exponent += exponentChange;
2200    }
2201  }
2202
2203  // If this is a truncation, perform the shift before we narrow the storage.
2204  if (shift < 0 && (isFiniteNonZero() || category==fcNaN))
2205    lostFraction = shiftRight(significandParts(), oldPartCount, -shift);
2206
2207  // Fix the storage so it can hold to new value.
2208  if (newPartCount > oldPartCount) {
2209    // The new type requires more storage; make it available.
2210    integerPart *newParts;
2211    newParts = new integerPart[newPartCount];
2212    APInt::tcSet(newParts, 0, newPartCount);
2213    if (isFiniteNonZero() || category==fcNaN)
2214      APInt::tcAssign(newParts, significandParts(), oldPartCount);
2215    freeSignificand();
2216    significand.parts = newParts;
2217  } else if (newPartCount == 1 && oldPartCount != 1) {
2218    // Switch to built-in storage for a single part.
2219    integerPart newPart = 0;
2220    if (isFiniteNonZero() || category==fcNaN)
2221      newPart = significandParts()[0];
2222    freeSignificand();
2223    significand.part = newPart;
2224  }
2225
2226  // Now that we have the right storage, switch the semantics.
2227  semantics = &toSemantics;
2228
2229  // If this is an extension, perform the shift now that the storage is
2230  // available.
2231  if (shift > 0 && (isFiniteNonZero() || category==fcNaN))
2232    APInt::tcShiftLeft(significandParts(), newPartCount, shift);
2233
2234  if (isFiniteNonZero()) {
2235    fs = normalize(rounding_mode, lostFraction);
2236    *losesInfo = (fs != opOK);
2237  } else if (category == fcNaN) {
2238    *losesInfo = lostFraction != lfExactlyZero || X86SpecialNan;
2239
2240    // For x87 extended precision, we want to make a NaN, not a special NaN if
2241    // the input wasn't special either.
2242    if (!X86SpecialNan && semantics == &semX87DoubleExtended)
2243      APInt::tcSetBit(significandParts(), semantics->precision - 1);
2244
2245    // If we are truncating NaN, it is possible that we shifted out all of the
2246    // set bits in a signalling NaN payload. But NaN must remain NaN, so some
2247    // bit in the significand must be set (otherwise it is Inf).
2248    // This can only happen with sNaN. Set the 1st bit after the quiet bit,
2249    // so that we still have an sNaN.
2250    // FIXME: Set quiet and return opInvalidOp (on convert of any sNaN).
2251    //        But this requires fixing LLVM to parse 32-bit hex FP or ignoring
2252    //        conversions while parsing IR.
2253    if (APInt::tcIsZero(significandParts(), newPartCount)) {
2254      assert(shift < 0 && "Should not lose NaN payload on extend");
2255      assert(semantics->precision >= 3 && "Unexpectedly narrow significand");
2256      assert(*losesInfo && "Missing payload should have set lost info");
2257      APInt::tcSetBit(significandParts(), semantics->precision - 3);
2258    }
2259
2260    // gcc forces the Quiet bit on, which means (float)(double)(float_sNan)
2261    // does not give you back the same bits.  This is dubious, and we
2262    // don't currently do it.  You're really supposed to get
2263    // an invalid operation signal at runtime, but nobody does that.
2264    fs = opOK;
2265  } else {
2266    *losesInfo = false;
2267    fs = opOK;
2268  }
2269
2270  return fs;
2271}
2272
2273/* Convert a floating point number to an integer according to the
2274   rounding mode.  If the rounded integer value is out of range this
2275   returns an invalid operation exception and the contents of the
2276   destination parts are unspecified.  If the rounded value is in
2277   range but the floating point number is not the exact integer, the C
2278   standard doesn't require an inexact exception to be raised.  IEEE
2279   854 does require it so we do that.
2280
2281   Note that for conversions to integer type the C standard requires
2282   round-to-zero to always be used.  */
2283IEEEFloat::opStatus IEEEFloat::convertToSignExtendedInteger(
2284    MutableArrayRef<integerPart> parts, unsigned int width, bool isSigned,
2285    roundingMode rounding_mode, bool *isExact) const {
2286  lostFraction lost_fraction;
2287  const integerPart *src;
2288  unsigned int dstPartsCount, truncatedBits;
2289
2290  *isExact = false;
2291
2292  /* Handle the three special cases first.  */
2293  if (category == fcInfinity || category == fcNaN)
2294    return opInvalidOp;
2295
2296  dstPartsCount = partCountForBits(width);
2297  assert(dstPartsCount <= parts.size() && "Integer too big");
2298
2299  if (category == fcZero) {
2300    APInt::tcSet(parts.data(), 0, dstPartsCount);
2301    // Negative zero can't be represented as an int.
2302    *isExact = !sign;
2303    return opOK;
2304  }
2305
2306  src = significandParts();
2307
2308  /* Step 1: place our absolute value, with any fraction truncated, in
2309     the destination.  */
2310  if (exponent < 0) {
2311    /* Our absolute value is less than one; truncate everything.  */
2312    APInt::tcSet(parts.data(), 0, dstPartsCount);
2313    /* For exponent -1 the integer bit represents .5, look at that.
2314       For smaller exponents leftmost truncated bit is 0. */
2315    truncatedBits = semantics->precision -1U - exponent;
2316  } else {
2317    /* We want the most significant (exponent + 1) bits; the rest are
2318       truncated.  */
2319    unsigned int bits = exponent + 1U;
2320
2321    /* Hopelessly large in magnitude?  */
2322    if (bits > width)
2323      return opInvalidOp;
2324
2325    if (bits < semantics->precision) {
2326      /* We truncate (semantics->precision - bits) bits.  */
2327      truncatedBits = semantics->precision - bits;
2328      APInt::tcExtract(parts.data(), dstPartsCount, src, bits, truncatedBits);
2329    } else {
2330      /* We want at least as many bits as are available.  */
2331      APInt::tcExtract(parts.data(), dstPartsCount, src, semantics->precision,
2332                       0);
2333      APInt::tcShiftLeft(parts.data(), dstPartsCount,
2334                         bits - semantics->precision);
2335      truncatedBits = 0;
2336    }
2337  }
2338
2339  /* Step 2: work out any lost fraction, and increment the absolute
2340     value if we would round away from zero.  */
2341  if (truncatedBits) {
2342    lost_fraction = lostFractionThroughTruncation(src, partCount(),
2343                                                  truncatedBits);
2344    if (lost_fraction != lfExactlyZero &&
2345        roundAwayFromZero(rounding_mode, lost_fraction, truncatedBits)) {
2346      if (APInt::tcIncrement(parts.data(), dstPartsCount))
2347        return opInvalidOp;     /* Overflow.  */
2348    }
2349  } else {
2350    lost_fraction = lfExactlyZero;
2351  }
2352
2353  /* Step 3: check if we fit in the destination.  */
2354  unsigned int omsb = APInt::tcMSB(parts.data(), dstPartsCount) + 1;
2355
2356  if (sign) {
2357    if (!isSigned) {
2358      /* Negative numbers cannot be represented as unsigned.  */
2359      if (omsb != 0)
2360        return opInvalidOp;
2361    } else {
2362      /* It takes omsb bits to represent the unsigned integer value.
2363         We lose a bit for the sign, but care is needed as the
2364         maximally negative integer is a special case.  */
2365      if (omsb == width &&
2366          APInt::tcLSB(parts.data(), dstPartsCount) + 1 != omsb)
2367        return opInvalidOp;
2368
2369      /* This case can happen because of rounding.  */
2370      if (omsb > width)
2371        return opInvalidOp;
2372    }
2373
2374    APInt::tcNegate (parts.data(), dstPartsCount);
2375  } else {
2376    if (omsb >= width + !isSigned)
2377      return opInvalidOp;
2378  }
2379
2380  if (lost_fraction == lfExactlyZero) {
2381    *isExact = true;
2382    return opOK;
2383  } else
2384    return opInexact;
2385}
2386
2387/* Same as convertToSignExtendedInteger, except we provide
2388   deterministic values in case of an invalid operation exception,
2389   namely zero for NaNs and the minimal or maximal value respectively
2390   for underflow or overflow.
2391   The *isExact output tells whether the result is exact, in the sense
2392   that converting it back to the original floating point type produces
2393   the original value.  This is almost equivalent to result==opOK,
2394   except for negative zeroes.
2395*/
2396IEEEFloat::opStatus
2397IEEEFloat::convertToInteger(MutableArrayRef<integerPart> parts,
2398                            unsigned int width, bool isSigned,
2399                            roundingMode rounding_mode, bool *isExact) const {
2400  opStatus fs;
2401
2402  fs = convertToSignExtendedInteger(parts, width, isSigned, rounding_mode,
2403                                    isExact);
2404
2405  if (fs == opInvalidOp) {
2406    unsigned int bits, dstPartsCount;
2407
2408    dstPartsCount = partCountForBits(width);
2409    assert(dstPartsCount <= parts.size() && "Integer too big");
2410
2411    if (category == fcNaN)
2412      bits = 0;
2413    else if (sign)
2414      bits = isSigned;
2415    else
2416      bits = width - isSigned;
2417
2418    APInt::tcSetLeastSignificantBits(parts.data(), dstPartsCount, bits);
2419    if (sign && isSigned)
2420      APInt::tcShiftLeft(parts.data(), dstPartsCount, width - 1);
2421  }
2422
2423  return fs;
2424}
2425
2426/* Convert an unsigned integer SRC to a floating point number,
2427   rounding according to ROUNDING_MODE.  The sign of the floating
2428   point number is not modified.  */
2429IEEEFloat::opStatus IEEEFloat::convertFromUnsignedParts(
2430    const integerPart *src, unsigned int srcCount, roundingMode rounding_mode) {
2431  unsigned int omsb, precision, dstCount;
2432  integerPart *dst;
2433  lostFraction lost_fraction;
2434
2435  category = fcNormal;
2436  omsb = APInt::tcMSB(src, srcCount) + 1;
2437  dst = significandParts();
2438  dstCount = partCount();
2439  precision = semantics->precision;
2440
2441  /* We want the most significant PRECISION bits of SRC.  There may not
2442     be that many; extract what we can.  */
2443  if (precision <= omsb) {
2444    exponent = omsb - 1;
2445    lost_fraction = lostFractionThroughTruncation(src, srcCount,
2446                                                  omsb - precision);
2447    APInt::tcExtract(dst, dstCount, src, precision, omsb - precision);
2448  } else {
2449    exponent = precision - 1;
2450    lost_fraction = lfExactlyZero;
2451    APInt::tcExtract(dst, dstCount, src, omsb, 0);
2452  }
2453
2454  return normalize(rounding_mode, lost_fraction);
2455}
2456
2457IEEEFloat::opStatus IEEEFloat::convertFromAPInt(const APInt &Val, bool isSigned,
2458                                                roundingMode rounding_mode) {
2459  unsigned int partCount = Val.getNumWords();
2460  APInt api = Val;
2461
2462  sign = false;
2463  if (isSigned && api.isNegative()) {
2464    sign = true;
2465    api = -api;
2466  }
2467
2468  return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode);
2469}
2470
2471/* Convert a two's complement integer SRC to a floating point number,
2472   rounding according to ROUNDING_MODE.  ISSIGNED is true if the
2473   integer is signed, in which case it must be sign-extended.  */
2474IEEEFloat::opStatus
2475IEEEFloat::convertFromSignExtendedInteger(const integerPart *src,
2476                                          unsigned int srcCount, bool isSigned,
2477                                          roundingMode rounding_mode) {
2478  opStatus status;
2479
2480  if (isSigned &&
2481      APInt::tcExtractBit(src, srcCount * integerPartWidth - 1)) {
2482    integerPart *copy;
2483
2484    /* If we're signed and negative negate a copy.  */
2485    sign = true;
2486    copy = new integerPart[srcCount];
2487    APInt::tcAssign(copy, src, srcCount);
2488    APInt::tcNegate(copy, srcCount);
2489    status = convertFromUnsignedParts(copy, srcCount, rounding_mode);
2490    delete [] copy;
2491  } else {
2492    sign = false;
2493    status = convertFromUnsignedParts(src, srcCount, rounding_mode);
2494  }
2495
2496  return status;
2497}
2498
2499/* FIXME: should this just take a const APInt reference?  */
2500IEEEFloat::opStatus
2501IEEEFloat::convertFromZeroExtendedInteger(const integerPart *parts,
2502                                          unsigned int width, bool isSigned,
2503                                          roundingMode rounding_mode) {
2504  unsigned int partCount = partCountForBits(width);
2505  APInt api = APInt(width, makeArrayRef(parts, partCount));
2506
2507  sign = false;
2508  if (isSigned && APInt::tcExtractBit(parts, width - 1)) {
2509    sign = true;
2510    api = -api;
2511  }
2512
2513  return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode);
2514}
2515
2516Expected<IEEEFloat::opStatus>
2517IEEEFloat::convertFromHexadecimalString(StringRef s,
2518                                        roundingMode rounding_mode) {
2519  lostFraction lost_fraction = lfExactlyZero;
2520
2521  category = fcNormal;
2522  zeroSignificand();
2523  exponent = 0;
2524
2525  integerPart *significand = significandParts();
2526  unsigned partsCount = partCount();
2527  unsigned bitPos = partsCount * integerPartWidth;
2528  bool computedTrailingFraction = false;
2529
2530  // Skip leading zeroes and any (hexa)decimal point.
2531  StringRef::iterator begin = s.begin();
2532  StringRef::iterator end = s.end();
2533  StringRef::iterator dot;
2534  auto PtrOrErr = skipLeadingZeroesAndAnyDot(begin, end, &dot);
2535  if (!PtrOrErr)
2536    return PtrOrErr.takeError();
2537  StringRef::iterator p = *PtrOrErr;
2538  StringRef::iterator firstSignificantDigit = p;
2539
2540  while (p != end) {
2541    integerPart hex_value;
2542
2543    if (*p == '.') {
2544      if (dot != end)
2545        return createError("String contains multiple dots");
2546      dot = p++;
2547      continue;
2548    }
2549
2550    hex_value = hexDigitValue(*p);
2551    if (hex_value == -1U)
2552      break;
2553
2554    p++;
2555
2556    // Store the number while we have space.
2557    if (bitPos) {
2558      bitPos -= 4;
2559      hex_value <<= bitPos % integerPartWidth;
2560      significand[bitPos / integerPartWidth] |= hex_value;
2561    } else if (!computedTrailingFraction) {
2562      auto FractOrErr = trailingHexadecimalFraction(p, end, hex_value);
2563      if (!FractOrErr)
2564        return FractOrErr.takeError();
2565      lost_fraction = *FractOrErr;
2566      computedTrailingFraction = true;
2567    }
2568  }
2569
2570  /* Hex floats require an exponent but not a hexadecimal point.  */
2571  if (p == end)
2572    return createError("Hex strings require an exponent");
2573  if (*p != 'p' && *p != 'P')
2574    return createError("Invalid character in significand");
2575  if (p == begin)
2576    return createError("Significand has no digits");
2577  if (dot != end && p - begin == 1)
2578    return createError("Significand has no digits");
2579
2580  /* Ignore the exponent if we are zero.  */
2581  if (p != firstSignificantDigit) {
2582    int expAdjustment;
2583
2584    /* Implicit hexadecimal point?  */
2585    if (dot == end)
2586      dot = p;
2587
2588    /* Calculate the exponent adjustment implicit in the number of
2589       significant digits.  */
2590    expAdjustment = static_cast<int>(dot - firstSignificantDigit);
2591    if (expAdjustment < 0)
2592      expAdjustment++;
2593    expAdjustment = expAdjustment * 4 - 1;
2594
2595    /* Adjust for writing the significand starting at the most
2596       significant nibble.  */
2597    expAdjustment += semantics->precision;
2598    expAdjustment -= partsCount * integerPartWidth;
2599
2600    /* Adjust for the given exponent.  */
2601    auto ExpOrErr = totalExponent(p + 1, end, expAdjustment);
2602    if (!ExpOrErr)
2603      return ExpOrErr.takeError();
2604    exponent = *ExpOrErr;
2605  }
2606
2607  return normalize(rounding_mode, lost_fraction);
2608}
2609
2610IEEEFloat::opStatus
2611IEEEFloat::roundSignificandWithExponent(const integerPart *decSigParts,
2612                                        unsigned sigPartCount, int exp,
2613                                        roundingMode rounding_mode) {
2614  unsigned int parts, pow5PartCount;
2615  fltSemantics calcSemantics = { 32767, -32767, 0, 0 };
2616  integerPart pow5Parts[maxPowerOfFiveParts];
2617  bool isNearest;
2618
2619  isNearest = (rounding_mode == rmNearestTiesToEven ||
2620               rounding_mode == rmNearestTiesToAway);
2621
2622  parts = partCountForBits(semantics->precision + 11);
2623
2624  /* Calculate pow(5, abs(exp)).  */
2625  pow5PartCount = powerOf5(pow5Parts, exp >= 0 ? exp: -exp);
2626
2627  for (;; parts *= 2) {
2628    opStatus sigStatus, powStatus;
2629    unsigned int excessPrecision, truncatedBits;
2630
2631    calcSemantics.precision = parts * integerPartWidth - 1;
2632    excessPrecision = calcSemantics.precision - semantics->precision;
2633    truncatedBits = excessPrecision;
2634
2635    IEEEFloat decSig(calcSemantics, uninitialized);
2636    decSig.makeZero(sign);
2637    IEEEFloat pow5(calcSemantics);
2638
2639    sigStatus = decSig.convertFromUnsignedParts(decSigParts, sigPartCount,
2640                                                rmNearestTiesToEven);
2641    powStatus = pow5.convertFromUnsignedParts(pow5Parts, pow5PartCount,
2642                                              rmNearestTiesToEven);
2643    /* Add exp, as 10^n = 5^n * 2^n.  */
2644    decSig.exponent += exp;
2645
2646    lostFraction calcLostFraction;
2647    integerPart HUerr, HUdistance;
2648    unsigned int powHUerr;
2649
2650    if (exp >= 0) {
2651      /* multiplySignificand leaves the precision-th bit set to 1.  */
2652      calcLostFraction = decSig.multiplySignificand(pow5);
2653      powHUerr = powStatus != opOK;
2654    } else {
2655      calcLostFraction = decSig.divideSignificand(pow5);
2656      /* Denormal numbers have less precision.  */
2657      if (decSig.exponent < semantics->minExponent) {
2658        excessPrecision += (semantics->minExponent - decSig.exponent);
2659        truncatedBits = excessPrecision;
2660        if (excessPrecision > calcSemantics.precision)
2661          excessPrecision = calcSemantics.precision;
2662      }
2663      /* Extra half-ulp lost in reciprocal of exponent.  */
2664      powHUerr = (powStatus == opOK && calcLostFraction == lfExactlyZero) ? 0:2;
2665    }
2666
2667    /* Both multiplySignificand and divideSignificand return the
2668       result with the integer bit set.  */
2669    assert(APInt::tcExtractBit
2670           (decSig.significandParts(), calcSemantics.precision - 1) == 1);
2671
2672    HUerr = HUerrBound(calcLostFraction != lfExactlyZero, sigStatus != opOK,
2673                       powHUerr);
2674    HUdistance = 2 * ulpsFromBoundary(decSig.significandParts(),
2675                                      excessPrecision, isNearest);
2676
2677    /* Are we guaranteed to round correctly if we truncate?  */
2678    if (HUdistance >= HUerr) {
2679      APInt::tcExtract(significandParts(), partCount(), decSig.significandParts(),
2680                       calcSemantics.precision - excessPrecision,
2681                       excessPrecision);
2682      /* Take the exponent of decSig.  If we tcExtract-ed less bits
2683         above we must adjust our exponent to compensate for the
2684         implicit right shift.  */
2685      exponent = (decSig.exponent + semantics->precision
2686                  - (calcSemantics.precision - excessPrecision));
2687      calcLostFraction = lostFractionThroughTruncation(decSig.significandParts(),
2688                                                       decSig.partCount(),
2689                                                       truncatedBits);
2690      return normalize(rounding_mode, calcLostFraction);
2691    }
2692  }
2693}
2694
2695Expected<IEEEFloat::opStatus>
2696IEEEFloat::convertFromDecimalString(StringRef str, roundingMode rounding_mode) {
2697  decimalInfo D;
2698  opStatus fs;
2699
2700  /* Scan the text.  */
2701  StringRef::iterator p = str.begin();
2702  if (Error Err = interpretDecimal(p, str.end(), &D))
2703    return std::move(Err);
2704
2705  /* Handle the quick cases.  First the case of no significant digits,
2706     i.e. zero, and then exponents that are obviously too large or too
2707     small.  Writing L for log 10 / log 2, a number d.ddddd*10^exp
2708     definitely overflows if
2709
2710           (exp - 1) * L >= maxExponent
2711
2712     and definitely underflows to zero where
2713
2714           (exp + 1) * L <= minExponent - precision
2715
2716     With integer arithmetic the tightest bounds for L are
2717
2718           93/28 < L < 196/59            [ numerator <= 256 ]
2719           42039/12655 < L < 28738/8651  [ numerator <= 65536 ]
2720  */
2721
2722  // Test if we have a zero number allowing for strings with no null terminators
2723  // and zero decimals with non-zero exponents.
2724  //
2725  // We computed firstSigDigit by ignoring all zeros and dots. Thus if
2726  // D->firstSigDigit equals str.end(), every digit must be a zero and there can
2727  // be at most one dot. On the other hand, if we have a zero with a non-zero
2728  // exponent, then we know that D.firstSigDigit will be non-numeric.
2729  if (D.firstSigDigit == str.end() || decDigitValue(*D.firstSigDigit) >= 10U) {
2730    category = fcZero;
2731    fs = opOK;
2732
2733  /* Check whether the normalized exponent is high enough to overflow
2734     max during the log-rebasing in the max-exponent check below. */
2735  } else if (D.normalizedExponent - 1 > INT_MAX / 42039) {
2736    fs = handleOverflow(rounding_mode);
2737
2738  /* If it wasn't, then it also wasn't high enough to overflow max
2739     during the log-rebasing in the min-exponent check.  Check that it
2740     won't overflow min in either check, then perform the min-exponent
2741     check. */
2742  } else if (D.normalizedExponent - 1 < INT_MIN / 42039 ||
2743             (D.normalizedExponent + 1) * 28738 <=
2744               8651 * (semantics->minExponent - (int) semantics->precision)) {
2745    /* Underflow to zero and round.  */
2746    category = fcNormal;
2747    zeroSignificand();
2748    fs = normalize(rounding_mode, lfLessThanHalf);
2749
2750  /* We can finally safely perform the max-exponent check. */
2751  } else if ((D.normalizedExponent - 1) * 42039
2752             >= 12655 * semantics->maxExponent) {
2753    /* Overflow and round.  */
2754    fs = handleOverflow(rounding_mode);
2755  } else {
2756    integerPart *decSignificand;
2757    unsigned int partCount;
2758
2759    /* A tight upper bound on number of bits required to hold an
2760       N-digit decimal integer is N * 196 / 59.  Allocate enough space
2761       to hold the full significand, and an extra part required by
2762       tcMultiplyPart.  */
2763    partCount = static_cast<unsigned int>(D.lastSigDigit - D.firstSigDigit) + 1;
2764    partCount = partCountForBits(1 + 196 * partCount / 59);
2765    decSignificand = new integerPart[partCount + 1];
2766    partCount = 0;
2767
2768    /* Convert to binary efficiently - we do almost all multiplication
2769       in an integerPart.  When this would overflow do we do a single
2770       bignum multiplication, and then revert again to multiplication
2771       in an integerPart.  */
2772    do {
2773      integerPart decValue, val, multiplier;
2774
2775      val = 0;
2776      multiplier = 1;
2777
2778      do {
2779        if (*p == '.') {
2780          p++;
2781          if (p == str.end()) {
2782            break;
2783          }
2784        }
2785        decValue = decDigitValue(*p++);
2786        if (decValue >= 10U) {
2787          delete[] decSignificand;
2788          return createError("Invalid character in significand");
2789        }
2790        multiplier *= 10;
2791        val = val * 10 + decValue;
2792        /* The maximum number that can be multiplied by ten with any
2793           digit added without overflowing an integerPart.  */
2794      } while (p <= D.lastSigDigit && multiplier <= (~ (integerPart) 0 - 9) / 10);
2795
2796      /* Multiply out the current part.  */
2797      APInt::tcMultiplyPart(decSignificand, decSignificand, multiplier, val,
2798                            partCount, partCount + 1, false);
2799
2800      /* If we used another part (likely but not guaranteed), increase
2801         the count.  */
2802      if (decSignificand[partCount])
2803        partCount++;
2804    } while (p <= D.lastSigDigit);
2805
2806    category = fcNormal;
2807    fs = roundSignificandWithExponent(decSignificand, partCount,
2808                                      D.exponent, rounding_mode);
2809
2810    delete [] decSignificand;
2811  }
2812
2813  return fs;
2814}
2815
2816bool IEEEFloat::convertFromStringSpecials(StringRef str) {
2817  const size_t MIN_NAME_SIZE = 3;
2818
2819  if (str.size() < MIN_NAME_SIZE)
2820    return false;
2821
2822  if (str.equals("inf") || str.equals("INFINITY") || str.equals("+Inf")) {
2823    makeInf(false);
2824    return true;
2825  }
2826
2827  bool IsNegative = str.front() == '-';
2828  if (IsNegative) {
2829    str = str.drop_front();
2830    if (str.size() < MIN_NAME_SIZE)
2831      return false;
2832
2833    if (str.equals("inf") || str.equals("INFINITY") || str.equals("Inf")) {
2834      makeInf(true);
2835      return true;
2836    }
2837  }
2838
2839  // If we have a 's' (or 'S') prefix, then this is a Signaling NaN.
2840  bool IsSignaling = str.front() == 's' || str.front() == 'S';
2841  if (IsSignaling) {
2842    str = str.drop_front();
2843    if (str.size() < MIN_NAME_SIZE)
2844      return false;
2845  }
2846
2847  if (str.startswith("nan") || str.startswith("NaN")) {
2848    str = str.drop_front(3);
2849
2850    // A NaN without payload.
2851    if (str.empty()) {
2852      makeNaN(IsSignaling, IsNegative);
2853      return true;
2854    }
2855
2856    // Allow the payload to be inside parentheses.
2857    if (str.front() == '(') {
2858      // Parentheses should be balanced (and not empty).
2859      if (str.size() <= 2 || str.back() != ')')
2860        return false;
2861
2862      str = str.slice(1, str.size() - 1);
2863    }
2864
2865    // Determine the payload number's radix.
2866    unsigned Radix = 10;
2867    if (str[0] == '0') {
2868      if (str.size() > 1 && tolower(str[1]) == 'x') {
2869        str = str.drop_front(2);
2870        Radix = 16;
2871      } else
2872        Radix = 8;
2873    }
2874
2875    // Parse the payload and make the NaN.
2876    APInt Payload;
2877    if (!str.getAsInteger(Radix, Payload)) {
2878      makeNaN(IsSignaling, IsNegative, &Payload);
2879      return true;
2880    }
2881  }
2882
2883  return false;
2884}
2885
2886Expected<IEEEFloat::opStatus>
2887IEEEFloat::convertFromString(StringRef str, roundingMode rounding_mode) {
2888  if (str.empty())
2889    return createError("Invalid string length");
2890
2891  // Handle special cases.
2892  if (convertFromStringSpecials(str))
2893    return opOK;
2894
2895  /* Handle a leading minus sign.  */
2896  StringRef::iterator p = str.begin();
2897  size_t slen = str.size();
2898  sign = *p == '-' ? 1 : 0;
2899  if (*p == '-' || *p == '+') {
2900    p++;
2901    slen--;
2902    if (!slen)
2903      return createError("String has no digits");
2904  }
2905
2906  if (slen >= 2 && p[0] == '0' && (p[1] == 'x' || p[1] == 'X')) {
2907    if (slen == 2)
2908      return createError("Invalid string");
2909    return convertFromHexadecimalString(StringRef(p + 2, slen - 2),
2910                                        rounding_mode);
2911  }
2912
2913  return convertFromDecimalString(StringRef(p, slen), rounding_mode);
2914}
2915
2916/* Write out a hexadecimal representation of the floating point value
2917   to DST, which must be of sufficient size, in the C99 form
2918   [-]0xh.hhhhp[+-]d.  Return the number of characters written,
2919   excluding the terminating NUL.
2920
2921   If UPPERCASE, the output is in upper case, otherwise in lower case.
2922
2923   HEXDIGITS digits appear altogether, rounding the value if
2924   necessary.  If HEXDIGITS is 0, the minimal precision to display the
2925   number precisely is used instead.  If nothing would appear after
2926   the decimal point it is suppressed.
2927
2928   The decimal exponent is always printed and has at least one digit.
2929   Zero values display an exponent of zero.  Infinities and NaNs
2930   appear as "infinity" or "nan" respectively.
2931
2932   The above rules are as specified by C99.  There is ambiguity about
2933   what the leading hexadecimal digit should be.  This implementation
2934   uses whatever is necessary so that the exponent is displayed as
2935   stored.  This implies the exponent will fall within the IEEE format
2936   range, and the leading hexadecimal digit will be 0 (for denormals),
2937   1 (normal numbers) or 2 (normal numbers rounded-away-from-zero with
2938   any other digits zero).
2939*/
2940unsigned int IEEEFloat::convertToHexString(char *dst, unsigned int hexDigits,
2941                                           bool upperCase,
2942                                           roundingMode rounding_mode) const {
2943  char *p;
2944
2945  p = dst;
2946  if (sign)
2947    *dst++ = '-';
2948
2949  switch (category) {
2950  case fcInfinity:
2951    memcpy (dst, upperCase ? infinityU: infinityL, sizeof infinityU - 1);
2952    dst += sizeof infinityL - 1;
2953    break;
2954
2955  case fcNaN:
2956    memcpy (dst, upperCase ? NaNU: NaNL, sizeof NaNU - 1);
2957    dst += sizeof NaNU - 1;
2958    break;
2959
2960  case fcZero:
2961    *dst++ = '0';
2962    *dst++ = upperCase ? 'X': 'x';
2963    *dst++ = '0';
2964    if (hexDigits > 1) {
2965      *dst++ = '.';
2966      memset (dst, '0', hexDigits - 1);
2967      dst += hexDigits - 1;
2968    }
2969    *dst++ = upperCase ? 'P': 'p';
2970    *dst++ = '0';
2971    break;
2972
2973  case fcNormal:
2974    dst = convertNormalToHexString (dst, hexDigits, upperCase, rounding_mode);
2975    break;
2976  }
2977
2978  *dst = 0;
2979
2980  return static_cast<unsigned int>(dst - p);
2981}
2982
2983/* Does the hard work of outputting the correctly rounded hexadecimal
2984   form of a normal floating point number with the specified number of
2985   hexadecimal digits.  If HEXDIGITS is zero the minimum number of
2986   digits necessary to print the value precisely is output.  */
2987char *IEEEFloat::convertNormalToHexString(char *dst, unsigned int hexDigits,
2988                                          bool upperCase,
2989                                          roundingMode rounding_mode) const {
2990  unsigned int count, valueBits, shift, partsCount, outputDigits;
2991  const char *hexDigitChars;
2992  const integerPart *significand;
2993  char *p;
2994  bool roundUp;
2995
2996  *dst++ = '0';
2997  *dst++ = upperCase ? 'X': 'x';
2998
2999  roundUp = false;
3000  hexDigitChars = upperCase ? hexDigitsUpper: hexDigitsLower;
3001
3002  significand = significandParts();
3003  partsCount = partCount();
3004
3005  /* +3 because the first digit only uses the single integer bit, so
3006     we have 3 virtual zero most-significant-bits.  */
3007  valueBits = semantics->precision + 3;
3008  shift = integerPartWidth - valueBits % integerPartWidth;
3009
3010  /* The natural number of digits required ignoring trailing
3011     insignificant zeroes.  */
3012  outputDigits = (valueBits - significandLSB () + 3) / 4;
3013
3014  /* hexDigits of zero means use the required number for the
3015     precision.  Otherwise, see if we are truncating.  If we are,
3016     find out if we need to round away from zero.  */
3017  if (hexDigits) {
3018    if (hexDigits < outputDigits) {
3019      /* We are dropping non-zero bits, so need to check how to round.
3020         "bits" is the number of dropped bits.  */
3021      unsigned int bits;
3022      lostFraction fraction;
3023
3024      bits = valueBits - hexDigits * 4;
3025      fraction = lostFractionThroughTruncation (significand, partsCount, bits);
3026      roundUp = roundAwayFromZero(rounding_mode, fraction, bits);
3027    }
3028    outputDigits = hexDigits;
3029  }
3030
3031  /* Write the digits consecutively, and start writing in the location
3032     of the hexadecimal point.  We move the most significant digit
3033     left and add the hexadecimal point later.  */
3034  p = ++dst;
3035
3036  count = (valueBits + integerPartWidth - 1) / integerPartWidth;
3037
3038  while (outputDigits && count) {
3039    integerPart part;
3040
3041    /* Put the most significant integerPartWidth bits in "part".  */
3042    if (--count == partsCount)
3043      part = 0;  /* An imaginary higher zero part.  */
3044    else
3045      part = significand[count] << shift;
3046
3047    if (count && shift)
3048      part |= significand[count - 1] >> (integerPartWidth - shift);
3049
3050    /* Convert as much of "part" to hexdigits as we can.  */
3051    unsigned int curDigits = integerPartWidth / 4;
3052
3053    if (curDigits > outputDigits)
3054      curDigits = outputDigits;
3055    dst += partAsHex (dst, part, curDigits, hexDigitChars);
3056    outputDigits -= curDigits;
3057  }
3058
3059  if (roundUp) {
3060    char *q = dst;
3061
3062    /* Note that hexDigitChars has a trailing '0'.  */
3063    do {
3064      q--;
3065      *q = hexDigitChars[hexDigitValue (*q) + 1];
3066    } while (*q == '0');
3067    assert(q >= p);
3068  } else {
3069    /* Add trailing zeroes.  */
3070    memset (dst, '0', outputDigits);
3071    dst += outputDigits;
3072  }
3073
3074  /* Move the most significant digit to before the point, and if there
3075     is something after the decimal point add it.  This must come
3076     after rounding above.  */
3077  p[-1] = p[0];
3078  if (dst -1 == p)
3079    dst--;
3080  else
3081    p[0] = '.';
3082
3083  /* Finally output the exponent.  */
3084  *dst++ = upperCase ? 'P': 'p';
3085
3086  return writeSignedDecimal (dst, exponent);
3087}
3088
3089hash_code hash_value(const IEEEFloat &Arg) {
3090  if (!Arg.isFiniteNonZero())
3091    return hash_combine((uint8_t)Arg.category,
3092                        // NaN has no sign, fix it at zero.
3093                        Arg.isNaN() ? (uint8_t)0 : (uint8_t)Arg.sign,
3094                        Arg.semantics->precision);
3095
3096  // Normal floats need their exponent and significand hashed.
3097  return hash_combine((uint8_t)Arg.category, (uint8_t)Arg.sign,
3098                      Arg.semantics->precision, Arg.exponent,
3099                      hash_combine_range(
3100                        Arg.significandParts(),
3101                        Arg.significandParts() + Arg.partCount()));
3102}
3103
3104// Conversion from APFloat to/from host float/double.  It may eventually be
3105// possible to eliminate these and have everybody deal with APFloats, but that
3106// will take a while.  This approach will not easily extend to long double.
3107// Current implementation requires integerPartWidth==64, which is correct at
3108// the moment but could be made more general.
3109
3110// Denormals have exponent minExponent in APFloat, but minExponent-1 in
3111// the actual IEEE respresentations.  We compensate for that here.
3112
3113APInt IEEEFloat::convertF80LongDoubleAPFloatToAPInt() const {
3114  assert(semantics == (const llvm::fltSemantics*)&semX87DoubleExtended);
3115  assert(partCount()==2);
3116
3117  uint64_t myexponent, mysignificand;
3118
3119  if (isFiniteNonZero()) {
3120    myexponent = exponent+16383; //bias
3121    mysignificand = significandParts()[0];
3122    if (myexponent==1 && !(mysignificand & 0x8000000000000000ULL))
3123      myexponent = 0;   // denormal
3124  } else if (category==fcZero) {
3125    myexponent = 0;
3126    mysignificand = 0;
3127  } else if (category==fcInfinity) {
3128    myexponent = 0x7fff;
3129    mysignificand = 0x8000000000000000ULL;
3130  } else {
3131    assert(category == fcNaN && "Unknown category");
3132    myexponent = 0x7fff;
3133    mysignificand = significandParts()[0];
3134  }
3135
3136  uint64_t words[2];
3137  words[0] = mysignificand;
3138  words[1] =  ((uint64_t)(sign & 1) << 15) |
3139              (myexponent & 0x7fffLL);
3140  return APInt(80, words);
3141}
3142
3143APInt IEEEFloat::convertPPCDoubleDoubleAPFloatToAPInt() const {
3144  assert(semantics == (const llvm::fltSemantics *)&semPPCDoubleDoubleLegacy);
3145  assert(partCount()==2);
3146
3147  uint64_t words[2];
3148  opStatus fs;
3149  bool losesInfo;
3150
3151  // Convert number to double.  To avoid spurious underflows, we re-
3152  // normalize against the "double" minExponent first, and only *then*
3153  // truncate the mantissa.  The result of that second conversion
3154  // may be inexact, but should never underflow.
3155  // Declare fltSemantics before APFloat that uses it (and
3156  // saves pointer to it) to ensure correct destruction order.
3157  fltSemantics extendedSemantics = *semantics;
3158  extendedSemantics.minExponent = semIEEEdouble.minExponent;
3159  IEEEFloat extended(*this);
3160  fs = extended.convert(extendedSemantics, rmNearestTiesToEven, &losesInfo);
3161  assert(fs == opOK && !losesInfo);
3162  (void)fs;
3163
3164  IEEEFloat u(extended);
3165  fs = u.convert(semIEEEdouble, rmNearestTiesToEven, &losesInfo);
3166  assert(fs == opOK || fs == opInexact);
3167  (void)fs;
3168  words[0] = *u.convertDoubleAPFloatToAPInt().getRawData();
3169
3170  // If conversion was exact or resulted in a special case, we're done;
3171  // just set the second double to zero.  Otherwise, re-convert back to
3172  // the extended format and compute the difference.  This now should
3173  // convert exactly to double.
3174  if (u.isFiniteNonZero() && losesInfo) {
3175    fs = u.convert(extendedSemantics, rmNearestTiesToEven, &losesInfo);
3176    assert(fs == opOK && !losesInfo);
3177    (void)fs;
3178
3179    IEEEFloat v(extended);
3180    v.subtract(u, rmNearestTiesToEven);
3181    fs = v.convert(semIEEEdouble, rmNearestTiesToEven, &losesInfo);
3182    assert(fs == opOK && !losesInfo);
3183    (void)fs;
3184    words[1] = *v.convertDoubleAPFloatToAPInt().getRawData();
3185  } else {
3186    words[1] = 0;
3187  }
3188
3189  return APInt(128, words);
3190}
3191
3192APInt IEEEFloat::convertQuadrupleAPFloatToAPInt() const {
3193  assert(semantics == (const llvm::fltSemantics*)&semIEEEquad);
3194  assert(partCount()==2);
3195
3196  uint64_t myexponent, mysignificand, mysignificand2;
3197
3198  if (isFiniteNonZero()) {
3199    myexponent = exponent+16383; //bias
3200    mysignificand = significandParts()[0];
3201    mysignificand2 = significandParts()[1];
3202    if (myexponent==1 && !(mysignificand2 & 0x1000000000000LL))
3203      myexponent = 0;   // denormal
3204  } else if (category==fcZero) {
3205    myexponent = 0;
3206    mysignificand = mysignificand2 = 0;
3207  } else if (category==fcInfinity) {
3208    myexponent = 0x7fff;
3209    mysignificand = mysignificand2 = 0;
3210  } else {
3211    assert(category == fcNaN && "Unknown category!");
3212    myexponent = 0x7fff;
3213    mysignificand = significandParts()[0];
3214    mysignificand2 = significandParts()[1];
3215  }
3216
3217  uint64_t words[2];
3218  words[0] = mysignificand;
3219  words[1] = ((uint64_t)(sign & 1) << 63) |
3220             ((myexponent & 0x7fff) << 48) |
3221             (mysignificand2 & 0xffffffffffffLL);
3222
3223  return APInt(128, words);
3224}
3225
3226APInt IEEEFloat::convertDoubleAPFloatToAPInt() const {
3227  assert(semantics == (const llvm::fltSemantics*)&semIEEEdouble);
3228  assert(partCount()==1);
3229
3230  uint64_t myexponent, mysignificand;
3231
3232  if (isFiniteNonZero()) {
3233    myexponent = exponent+1023; //bias
3234    mysignificand = *significandParts();
3235    if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
3236      myexponent = 0;   // denormal
3237  } else if (category==fcZero) {
3238    myexponent = 0;
3239    mysignificand = 0;
3240  } else if (category==fcInfinity) {
3241    myexponent = 0x7ff;
3242    mysignificand = 0;
3243  } else {
3244    assert(category == fcNaN && "Unknown category!");
3245    myexponent = 0x7ff;
3246    mysignificand = *significandParts();
3247  }
3248
3249  return APInt(64, ((((uint64_t)(sign & 1) << 63) |
3250                     ((myexponent & 0x7ff) <<  52) |
3251                     (mysignificand & 0xfffffffffffffLL))));
3252}
3253
3254APInt IEEEFloat::convertFloatAPFloatToAPInt() const {
3255  assert(semantics == (const llvm::fltSemantics*)&semIEEEsingle);
3256  assert(partCount()==1);
3257
3258  uint32_t myexponent, mysignificand;
3259
3260  if (isFiniteNonZero()) {
3261    myexponent = exponent+127; //bias
3262    mysignificand = (uint32_t)*significandParts();
3263    if (myexponent == 1 && !(mysignificand & 0x800000))
3264      myexponent = 0;   // denormal
3265  } else if (category==fcZero) {
3266    myexponent = 0;
3267    mysignificand = 0;
3268  } else if (category==fcInfinity) {
3269    myexponent = 0xff;
3270    mysignificand = 0;
3271  } else {
3272    assert(category == fcNaN && "Unknown category!");
3273    myexponent = 0xff;
3274    mysignificand = (uint32_t)*significandParts();
3275  }
3276
3277  return APInt(32, (((sign&1) << 31) | ((myexponent&0xff) << 23) |
3278                    (mysignificand & 0x7fffff)));
3279}
3280
3281APInt IEEEFloat::convertBFloatAPFloatToAPInt() const {
3282  assert(semantics == (const llvm::fltSemantics *)&semBFloat);
3283  assert(partCount() == 1);
3284
3285  uint32_t myexponent, mysignificand;
3286
3287  if (isFiniteNonZero()) {
3288    myexponent = exponent + 127; // bias
3289    mysignificand = (uint32_t)*significandParts();
3290    if (myexponent == 1 && !(mysignificand & 0x80))
3291      myexponent = 0; // denormal
3292  } else if (category == fcZero) {
3293    myexponent = 0;
3294    mysignificand = 0;
3295  } else if (category == fcInfinity) {
3296    myexponent = 0xff;
3297    mysignificand = 0;
3298  } else {
3299    assert(category == fcNaN && "Unknown category!");
3300    myexponent = 0xff;
3301    mysignificand = (uint32_t)*significandParts();
3302  }
3303
3304  return APInt(16, (((sign & 1) << 15) | ((myexponent & 0xff) << 7) |
3305                    (mysignificand & 0x7f)));
3306}
3307
3308APInt IEEEFloat::convertHalfAPFloatToAPInt() const {
3309  assert(semantics == (const llvm::fltSemantics*)&semIEEEhalf);
3310  assert(partCount()==1);
3311
3312  uint32_t myexponent, mysignificand;
3313
3314  if (isFiniteNonZero()) {
3315    myexponent = exponent+15; //bias
3316    mysignificand = (uint32_t)*significandParts();
3317    if (myexponent == 1 && !(mysignificand & 0x400))
3318      myexponent = 0;   // denormal
3319  } else if (category==fcZero) {
3320    myexponent = 0;
3321    mysignificand = 0;
3322  } else if (category==fcInfinity) {
3323    myexponent = 0x1f;
3324    mysignificand = 0;
3325  } else {
3326    assert(category == fcNaN && "Unknown category!");
3327    myexponent = 0x1f;
3328    mysignificand = (uint32_t)*significandParts();
3329  }
3330
3331  return APInt(16, (((sign&1) << 15) | ((myexponent&0x1f) << 10) |
3332                    (mysignificand & 0x3ff)));
3333}
3334
3335// This function creates an APInt that is just a bit map of the floating
3336// point constant as it would appear in memory.  It is not a conversion,
3337// and treating the result as a normal integer is unlikely to be useful.
3338
3339APInt IEEEFloat::bitcastToAPInt() const {
3340  if (semantics == (const llvm::fltSemantics*)&semIEEEhalf)
3341    return convertHalfAPFloatToAPInt();
3342
3343  if (semantics == (const llvm::fltSemantics *)&semBFloat)
3344    return convertBFloatAPFloatToAPInt();
3345
3346  if (semantics == (const llvm::fltSemantics*)&semIEEEsingle)
3347    return convertFloatAPFloatToAPInt();
3348
3349  if (semantics == (const llvm::fltSemantics*)&semIEEEdouble)
3350    return convertDoubleAPFloatToAPInt();
3351
3352  if (semantics == (const llvm::fltSemantics*)&semIEEEquad)
3353    return convertQuadrupleAPFloatToAPInt();
3354
3355  if (semantics == (const llvm::fltSemantics *)&semPPCDoubleDoubleLegacy)
3356    return convertPPCDoubleDoubleAPFloatToAPInt();
3357
3358  assert(semantics == (const llvm::fltSemantics*)&semX87DoubleExtended &&
3359         "unknown format!");
3360  return convertF80LongDoubleAPFloatToAPInt();
3361}
3362
3363float IEEEFloat::convertToFloat() const {
3364  assert(semantics == (const llvm::fltSemantics*)&semIEEEsingle &&
3365         "Float semantics are not IEEEsingle");
3366  APInt api = bitcastToAPInt();
3367  return api.bitsToFloat();
3368}
3369
3370double IEEEFloat::convertToDouble() const {
3371  assert(semantics == (const llvm::fltSemantics*)&semIEEEdouble &&
3372         "Float semantics are not IEEEdouble");
3373  APInt api = bitcastToAPInt();
3374  return api.bitsToDouble();
3375}
3376
3377/// Integer bit is explicit in this format.  Intel hardware (387 and later)
3378/// does not support these bit patterns:
3379///  exponent = all 1's, integer bit 0, significand 0 ("pseudoinfinity")
3380///  exponent = all 1's, integer bit 0, significand nonzero ("pseudoNaN")
3381///  exponent!=0 nor all 1's, integer bit 0 ("unnormal")
3382///  exponent = 0, integer bit 1 ("pseudodenormal")
3383/// At the moment, the first three are treated as NaNs, the last one as Normal.
3384void IEEEFloat::initFromF80LongDoubleAPInt(const APInt &api) {
3385  assert(api.getBitWidth()==80);
3386  uint64_t i1 = api.getRawData()[0];
3387  uint64_t i2 = api.getRawData()[1];
3388  uint64_t myexponent = (i2 & 0x7fff);
3389  uint64_t mysignificand = i1;
3390  uint8_t myintegerbit = mysignificand >> 63;
3391
3392  initialize(&semX87DoubleExtended);
3393  assert(partCount()==2);
3394
3395  sign = static_cast<unsigned int>(i2>>15);
3396  if (myexponent == 0 && mysignificand == 0) {
3397    // exponent, significand meaningless
3398    category = fcZero;
3399  } else if (myexponent==0x7fff && mysignificand==0x8000000000000000ULL) {
3400    // exponent, significand meaningless
3401    category = fcInfinity;
3402  } else if ((myexponent == 0x7fff && mysignificand != 0x8000000000000000ULL) ||
3403             (myexponent != 0x7fff && myexponent != 0 && myintegerbit == 0)) {
3404    // exponent meaningless
3405    category = fcNaN;
3406    significandParts()[0] = mysignificand;
3407    significandParts()[1] = 0;
3408  } else {
3409    category = fcNormal;
3410    exponent = myexponent - 16383;
3411    significandParts()[0] = mysignificand;
3412    significandParts()[1] = 0;
3413    if (myexponent==0)          // denormal
3414      exponent = -16382;
3415  }
3416}
3417
3418void IEEEFloat::initFromPPCDoubleDoubleAPInt(const APInt &api) {
3419  assert(api.getBitWidth()==128);
3420  uint64_t i1 = api.getRawData()[0];
3421  uint64_t i2 = api.getRawData()[1];
3422  opStatus fs;
3423  bool losesInfo;
3424
3425  // Get the first double and convert to our format.
3426  initFromDoubleAPInt(APInt(64, i1));
3427  fs = convert(semPPCDoubleDoubleLegacy, rmNearestTiesToEven, &losesInfo);
3428  assert(fs == opOK && !losesInfo);
3429  (void)fs;
3430
3431  // Unless we have a special case, add in second double.
3432  if (isFiniteNonZero()) {
3433    IEEEFloat v(semIEEEdouble, APInt(64, i2));
3434    fs = v.convert(semPPCDoubleDoubleLegacy, rmNearestTiesToEven, &losesInfo);
3435    assert(fs == opOK && !losesInfo);
3436    (void)fs;
3437
3438    add(v, rmNearestTiesToEven);
3439  }
3440}
3441
3442void IEEEFloat::initFromQuadrupleAPInt(const APInt &api) {
3443  assert(api.getBitWidth()==128);
3444  uint64_t i1 = api.getRawData()[0];
3445  uint64_t i2 = api.getRawData()[1];
3446  uint64_t myexponent = (i2 >> 48) & 0x7fff;
3447  uint64_t mysignificand  = i1;
3448  uint64_t mysignificand2 = i2 & 0xffffffffffffLL;
3449
3450  initialize(&semIEEEquad);
3451  assert(partCount()==2);
3452
3453  sign = static_cast<unsigned int>(i2>>63);
3454  if (myexponent==0 &&
3455      (mysignificand==0 && mysignificand2==0)) {
3456    // exponent, significand meaningless
3457    category = fcZero;
3458  } else if (myexponent==0x7fff &&
3459             (mysignificand==0 && mysignificand2==0)) {
3460    // exponent, significand meaningless
3461    category = fcInfinity;
3462  } else if (myexponent==0x7fff &&
3463             (mysignificand!=0 || mysignificand2 !=0)) {
3464    // exponent meaningless
3465    category = fcNaN;
3466    significandParts()[0] = mysignificand;
3467    significandParts()[1] = mysignificand2;
3468  } else {
3469    category = fcNormal;
3470    exponent = myexponent - 16383;
3471    significandParts()[0] = mysignificand;
3472    significandParts()[1] = mysignificand2;
3473    if (myexponent==0)          // denormal
3474      exponent = -16382;
3475    else
3476      significandParts()[1] |= 0x1000000000000LL;  // integer bit
3477  }
3478}
3479
3480void IEEEFloat::initFromDoubleAPInt(const APInt &api) {
3481  assert(api.getBitWidth()==64);
3482  uint64_t i = *api.getRawData();
3483  uint64_t myexponent = (i >> 52) & 0x7ff;
3484  uint64_t mysignificand = i & 0xfffffffffffffLL;
3485
3486  initialize(&semIEEEdouble);
3487  assert(partCount()==1);
3488
3489  sign = static_cast<unsigned int>(i>>63);
3490  if (myexponent==0 && mysignificand==0) {
3491    // exponent, significand meaningless
3492    category = fcZero;
3493  } else if (myexponent==0x7ff && mysignificand==0) {
3494    // exponent, significand meaningless
3495    category = fcInfinity;
3496  } else if (myexponent==0x7ff && mysignificand!=0) {
3497    // exponent meaningless
3498    category = fcNaN;
3499    *significandParts() = mysignificand;
3500  } else {
3501    category = fcNormal;
3502    exponent = myexponent - 1023;
3503    *significandParts() = mysignificand;
3504    if (myexponent==0)          // denormal
3505      exponent = -1022;
3506    else
3507      *significandParts() |= 0x10000000000000LL;  // integer bit
3508  }
3509}
3510
3511void IEEEFloat::initFromFloatAPInt(const APInt &api) {
3512  assert(api.getBitWidth()==32);
3513  uint32_t i = (uint32_t)*api.getRawData();
3514  uint32_t myexponent = (i >> 23) & 0xff;
3515  uint32_t mysignificand = i & 0x7fffff;
3516
3517  initialize(&semIEEEsingle);
3518  assert(partCount()==1);
3519
3520  sign = i >> 31;
3521  if (myexponent==0 && mysignificand==0) {
3522    // exponent, significand meaningless
3523    category = fcZero;
3524  } else if (myexponent==0xff && mysignificand==0) {
3525    // exponent, significand meaningless
3526    category = fcInfinity;
3527  } else if (myexponent==0xff && mysignificand!=0) {
3528    // sign, exponent, significand meaningless
3529    category = fcNaN;
3530    *significandParts() = mysignificand;
3531  } else {
3532    category = fcNormal;
3533    exponent = myexponent - 127;  //bias
3534    *significandParts() = mysignificand;
3535    if (myexponent==0)    // denormal
3536      exponent = -126;
3537    else
3538      *significandParts() |= 0x800000; // integer bit
3539  }
3540}
3541
3542void IEEEFloat::initFromBFloatAPInt(const APInt &api) {
3543  assert(api.getBitWidth() == 16);
3544  uint32_t i = (uint32_t)*api.getRawData();
3545  uint32_t myexponent = (i >> 7) & 0xff;
3546  uint32_t mysignificand = i & 0x7f;
3547
3548  initialize(&semBFloat);
3549  assert(partCount() == 1);
3550
3551  sign = i >> 15;
3552  if (myexponent == 0 && mysignificand == 0) {
3553    // exponent, significand meaningless
3554    category = fcZero;
3555  } else if (myexponent == 0xff && mysignificand == 0) {
3556    // exponent, significand meaningless
3557    category = fcInfinity;
3558  } else if (myexponent == 0xff && mysignificand != 0) {
3559    // sign, exponent, significand meaningless
3560    category = fcNaN;
3561    *significandParts() = mysignificand;
3562  } else {
3563    category = fcNormal;
3564    exponent = myexponent - 127; // bias
3565    *significandParts() = mysignificand;
3566    if (myexponent == 0) // denormal
3567      exponent = -126;
3568    else
3569      *significandParts() |= 0x80; // integer bit
3570  }
3571}
3572
3573void IEEEFloat::initFromHalfAPInt(const APInt &api) {
3574  assert(api.getBitWidth()==16);
3575  uint32_t i = (uint32_t)*api.getRawData();
3576  uint32_t myexponent = (i >> 10) & 0x1f;
3577  uint32_t mysignificand = i & 0x3ff;
3578
3579  initialize(&semIEEEhalf);
3580  assert(partCount()==1);
3581
3582  sign = i >> 15;
3583  if (myexponent==0 && mysignificand==0) {
3584    // exponent, significand meaningless
3585    category = fcZero;
3586  } else if (myexponent==0x1f && mysignificand==0) {
3587    // exponent, significand meaningless
3588    category = fcInfinity;
3589  } else if (myexponent==0x1f && mysignificand!=0) {
3590    // sign, exponent, significand meaningless
3591    category = fcNaN;
3592    *significandParts() = mysignificand;
3593  } else {
3594    category = fcNormal;
3595    exponent = myexponent - 15;  //bias
3596    *significandParts() = mysignificand;
3597    if (myexponent==0)    // denormal
3598      exponent = -14;
3599    else
3600      *significandParts() |= 0x400; // integer bit
3601  }
3602}
3603
3604/// Treat api as containing the bits of a floating point number.  Currently
3605/// we infer the floating point type from the size of the APInt.  The
3606/// isIEEE argument distinguishes between PPC128 and IEEE128 (not meaningful
3607/// when the size is anything else).
3608void IEEEFloat::initFromAPInt(const fltSemantics *Sem, const APInt &api) {
3609  if (Sem == &semIEEEhalf)
3610    return initFromHalfAPInt(api);
3611  if (Sem == &semBFloat)
3612    return initFromBFloatAPInt(api);
3613  if (Sem == &semIEEEsingle)
3614    return initFromFloatAPInt(api);
3615  if (Sem == &semIEEEdouble)
3616    return initFromDoubleAPInt(api);
3617  if (Sem == &semX87DoubleExtended)
3618    return initFromF80LongDoubleAPInt(api);
3619  if (Sem == &semIEEEquad)
3620    return initFromQuadrupleAPInt(api);
3621  if (Sem == &semPPCDoubleDoubleLegacy)
3622    return initFromPPCDoubleDoubleAPInt(api);
3623
3624  llvm_unreachable(nullptr);
3625}
3626
3627/// Make this number the largest magnitude normal number in the given
3628/// semantics.
3629void IEEEFloat::makeLargest(bool Negative) {
3630  // We want (in interchange format):
3631  //   sign = {Negative}
3632  //   exponent = 1..10
3633  //   significand = 1..1
3634  category = fcNormal;
3635  sign = Negative;
3636  exponent = semantics->maxExponent;
3637
3638  // Use memset to set all but the highest integerPart to all ones.
3639  integerPart *significand = significandParts();
3640  unsigned PartCount = partCount();
3641  memset(significand, 0xFF, sizeof(integerPart)*(PartCount - 1));
3642
3643  // Set the high integerPart especially setting all unused top bits for
3644  // internal consistency.
3645  const unsigned NumUnusedHighBits =
3646    PartCount*integerPartWidth - semantics->precision;
3647  significand[PartCount - 1] = (NumUnusedHighBits < integerPartWidth)
3648                                   ? (~integerPart(0) >> NumUnusedHighBits)
3649                                   : 0;
3650}
3651
3652/// Make this number the smallest magnitude denormal number in the given
3653/// semantics.
3654void IEEEFloat::makeSmallest(bool Negative) {
3655  // We want (in interchange format):
3656  //   sign = {Negative}
3657  //   exponent = 0..0
3658  //   significand = 0..01
3659  category = fcNormal;
3660  sign = Negative;
3661  exponent = semantics->minExponent;
3662  APInt::tcSet(significandParts(), 1, partCount());
3663}
3664
3665void IEEEFloat::makeSmallestNormalized(bool Negative) {
3666  // We want (in interchange format):
3667  //   sign = {Negative}
3668  //   exponent = 0..0
3669  //   significand = 10..0
3670
3671  category = fcNormal;
3672  zeroSignificand();
3673  sign = Negative;
3674  exponent = semantics->minExponent;
3675  significandParts()[partCountForBits(semantics->precision) - 1] |=
3676      (((integerPart)1) << ((semantics->precision - 1) % integerPartWidth));
3677}
3678
3679IEEEFloat::IEEEFloat(const fltSemantics &Sem, const APInt &API) {
3680  initFromAPInt(&Sem, API);
3681}
3682
3683IEEEFloat::IEEEFloat(float f) {
3684  initFromAPInt(&semIEEEsingle, APInt::floatToBits(f));
3685}
3686
3687IEEEFloat::IEEEFloat(double d) {
3688  initFromAPInt(&semIEEEdouble, APInt::doubleToBits(d));
3689}
3690
3691namespace {
3692  void append(SmallVectorImpl<char> &Buffer, StringRef Str) {
3693    Buffer.append(Str.begin(), Str.end());
3694  }
3695
3696  /// Removes data from the given significand until it is no more
3697  /// precise than is required for the desired precision.
3698  void AdjustToPrecision(APInt &significand,
3699                         int &exp, unsigned FormatPrecision) {
3700    unsigned bits = significand.getActiveBits();
3701
3702    // 196/59 is a very slight overestimate of lg_2(10).
3703    unsigned bitsRequired = (FormatPrecision * 196 + 58) / 59;
3704
3705    if (bits <= bitsRequired) return;
3706
3707    unsigned tensRemovable = (bits - bitsRequired) * 59 / 196;
3708    if (!tensRemovable) return;
3709
3710    exp += tensRemovable;
3711
3712    APInt divisor(significand.getBitWidth(), 1);
3713    APInt powten(significand.getBitWidth(), 10);
3714    while (true) {
3715      if (tensRemovable & 1)
3716        divisor *= powten;
3717      tensRemovable >>= 1;
3718      if (!tensRemovable) break;
3719      powten *= powten;
3720    }
3721
3722    significand = significand.udiv(divisor);
3723
3724    // Truncate the significand down to its active bit count.
3725    significand = significand.trunc(significand.getActiveBits());
3726  }
3727
3728
3729  void AdjustToPrecision(SmallVectorImpl<char> &buffer,
3730                         int &exp, unsigned FormatPrecision) {
3731    unsigned N = buffer.size();
3732    if (N <= FormatPrecision) return;
3733
3734    // The most significant figures are the last ones in the buffer.
3735    unsigned FirstSignificant = N - FormatPrecision;
3736
3737    // Round.
3738    // FIXME: this probably shouldn't use 'round half up'.
3739
3740    // Rounding down is just a truncation, except we also want to drop
3741    // trailing zeros from the new result.
3742    if (buffer[FirstSignificant - 1] < '5') {
3743      while (FirstSignificant < N && buffer[FirstSignificant] == '0')
3744        FirstSignificant++;
3745
3746      exp += FirstSignificant;
3747      buffer.erase(&buffer[0], &buffer[FirstSignificant]);
3748      return;
3749    }
3750
3751    // Rounding up requires a decimal add-with-carry.  If we continue
3752    // the carry, the newly-introduced zeros will just be truncated.
3753    for (unsigned I = FirstSignificant; I != N; ++I) {
3754      if (buffer[I] == '9') {
3755        FirstSignificant++;
3756      } else {
3757        buffer[I]++;
3758        break;
3759      }
3760    }
3761
3762    // If we carried through, we have exactly one digit of precision.
3763    if (FirstSignificant == N) {
3764      exp += FirstSignificant;
3765      buffer.clear();
3766      buffer.push_back('1');
3767      return;
3768    }
3769
3770    exp += FirstSignificant;
3771    buffer.erase(&buffer[0], &buffer[FirstSignificant]);
3772  }
3773}
3774
3775void IEEEFloat::toString(SmallVectorImpl<char> &Str, unsigned FormatPrecision,
3776                         unsigned FormatMaxPadding, bool TruncateZero) const {
3777  switch (category) {
3778  case fcInfinity:
3779    if (isNegative())
3780      return append(Str, "-Inf");
3781    else
3782      return append(Str, "+Inf");
3783
3784  case fcNaN: return append(Str, "NaN");
3785
3786  case fcZero:
3787    if (isNegative())
3788      Str.push_back('-');
3789
3790    if (!FormatMaxPadding) {
3791      if (TruncateZero)
3792        append(Str, "0.0E+0");
3793      else {
3794        append(Str, "0.0");
3795        if (FormatPrecision > 1)
3796          Str.append(FormatPrecision - 1, '0');
3797        append(Str, "e+00");
3798      }
3799    } else
3800      Str.push_back('0');
3801    return;
3802
3803  case fcNormal:
3804    break;
3805  }
3806
3807  if (isNegative())
3808    Str.push_back('-');
3809
3810  // Decompose the number into an APInt and an exponent.
3811  int exp = exponent - ((int) semantics->precision - 1);
3812  APInt significand(semantics->precision,
3813                    makeArrayRef(significandParts(),
3814                                 partCountForBits(semantics->precision)));
3815
3816  // Set FormatPrecision if zero.  We want to do this before we
3817  // truncate trailing zeros, as those are part of the precision.
3818  if (!FormatPrecision) {
3819    // We use enough digits so the number can be round-tripped back to an
3820    // APFloat. The formula comes from "How to Print Floating-Point Numbers
3821    // Accurately" by Steele and White.
3822    // FIXME: Using a formula based purely on the precision is conservative;
3823    // we can print fewer digits depending on the actual value being printed.
3824
3825    // FormatPrecision = 2 + floor(significandBits / lg_2(10))
3826    FormatPrecision = 2 + semantics->precision * 59 / 196;
3827  }
3828
3829  // Ignore trailing binary zeros.
3830  int trailingZeros = significand.countTrailingZeros();
3831  exp += trailingZeros;
3832  significand.lshrInPlace(trailingZeros);
3833
3834  // Change the exponent from 2^e to 10^e.
3835  if (exp == 0) {
3836    // Nothing to do.
3837  } else if (exp > 0) {
3838    // Just shift left.
3839    significand = significand.zext(semantics->precision + exp);
3840    significand <<= exp;
3841    exp = 0;
3842  } else { /* exp < 0 */
3843    int texp = -exp;
3844
3845    // We transform this using the identity:
3846    //   (N)(2^-e) == (N)(5^e)(10^-e)
3847    // This means we have to multiply N (the significand) by 5^e.
3848    // To avoid overflow, we have to operate on numbers large
3849    // enough to store N * 5^e:
3850    //   log2(N * 5^e) == log2(N) + e * log2(5)
3851    //                 <= semantics->precision + e * 137 / 59
3852    //   (log_2(5) ~ 2.321928 < 2.322034 ~ 137/59)
3853
3854    unsigned precision = semantics->precision + (137 * texp + 136) / 59;
3855
3856    // Multiply significand by 5^e.
3857    //   N * 5^0101 == N * 5^(1*1) * 5^(0*2) * 5^(1*4) * 5^(0*8)
3858    significand = significand.zext(precision);
3859    APInt five_to_the_i(precision, 5);
3860    while (true) {
3861      if (texp & 1) significand *= five_to_the_i;
3862
3863      texp >>= 1;
3864      if (!texp) break;
3865      five_to_the_i *= five_to_the_i;
3866    }
3867  }
3868
3869  AdjustToPrecision(significand, exp, FormatPrecision);
3870
3871  SmallVector<char, 256> buffer;
3872
3873  // Fill the buffer.
3874  unsigned precision = significand.getBitWidth();
3875  APInt ten(precision, 10);
3876  APInt digit(precision, 0);
3877
3878  bool inTrail = true;
3879  while (significand != 0) {
3880    // digit <- significand % 10
3881    // significand <- significand / 10
3882    APInt::udivrem(significand, ten, significand, digit);
3883
3884    unsigned d = digit.getZExtValue();
3885
3886    // Drop trailing zeros.
3887    if (inTrail && !d) exp++;
3888    else {
3889      buffer.push_back((char) ('0' + d));
3890      inTrail = false;
3891    }
3892  }
3893
3894  assert(!buffer.empty() && "no characters in buffer!");
3895
3896  // Drop down to FormatPrecision.
3897  // TODO: don't do more precise calculations above than are required.
3898  AdjustToPrecision(buffer, exp, FormatPrecision);
3899
3900  unsigned NDigits = buffer.size();
3901
3902  // Check whether we should use scientific notation.
3903  bool FormatScientific;
3904  if (!FormatMaxPadding)
3905    FormatScientific = true;
3906  else {
3907    if (exp >= 0) {
3908      // 765e3 --> 765000
3909      //              ^^^
3910      // But we shouldn't make the number look more precise than it is.
3911      FormatScientific = ((unsigned) exp > FormatMaxPadding ||
3912                          NDigits + (unsigned) exp > FormatPrecision);
3913    } else {
3914      // Power of the most significant digit.
3915      int MSD = exp + (int) (NDigits - 1);
3916      if (MSD >= 0) {
3917        // 765e-2 == 7.65
3918        FormatScientific = false;
3919      } else {
3920        // 765e-5 == 0.00765
3921        //           ^ ^^
3922        FormatScientific = ((unsigned) -MSD) > FormatMaxPadding;
3923      }
3924    }
3925  }
3926
3927  // Scientific formatting is pretty straightforward.
3928  if (FormatScientific) {
3929    exp += (NDigits - 1);
3930
3931    Str.push_back(buffer[NDigits-1]);
3932    Str.push_back('.');
3933    if (NDigits == 1 && TruncateZero)
3934      Str.push_back('0');
3935    else
3936      for (unsigned I = 1; I != NDigits; ++I)
3937        Str.push_back(buffer[NDigits-1-I]);
3938    // Fill with zeros up to FormatPrecision.
3939    if (!TruncateZero && FormatPrecision > NDigits - 1)
3940      Str.append(FormatPrecision - NDigits + 1, '0');
3941    // For !TruncateZero we use lower 'e'.
3942    Str.push_back(TruncateZero ? 'E' : 'e');
3943
3944    Str.push_back(exp >= 0 ? '+' : '-');
3945    if (exp < 0) exp = -exp;
3946    SmallVector<char, 6> expbuf;
3947    do {
3948      expbuf.push_back((char) ('0' + (exp % 10)));
3949      exp /= 10;
3950    } while (exp);
3951    // Exponent always at least two digits if we do not truncate zeros.
3952    if (!TruncateZero && expbuf.size() < 2)
3953      expbuf.push_back('0');
3954    for (unsigned I = 0, E = expbuf.size(); I != E; ++I)
3955      Str.push_back(expbuf[E-1-I]);
3956    return;
3957  }
3958
3959  // Non-scientific, positive exponents.
3960  if (exp >= 0) {
3961    for (unsigned I = 0; I != NDigits; ++I)
3962      Str.push_back(buffer[NDigits-1-I]);
3963    for (unsigned I = 0; I != (unsigned) exp; ++I)
3964      Str.push_back('0');
3965    return;
3966  }
3967
3968  // Non-scientific, negative exponents.
3969
3970  // The number of digits to the left of the decimal point.
3971  int NWholeDigits = exp + (int) NDigits;
3972
3973  unsigned I = 0;
3974  if (NWholeDigits > 0) {
3975    for (; I != (unsigned) NWholeDigits; ++I)
3976      Str.push_back(buffer[NDigits-I-1]);
3977    Str.push_back('.');
3978  } else {
3979    unsigned NZeros = 1 + (unsigned) -NWholeDigits;
3980
3981    Str.push_back('0');
3982    Str.push_back('.');
3983    for (unsigned Z = 1; Z != NZeros; ++Z)
3984      Str.push_back('0');
3985  }
3986
3987  for (; I != NDigits; ++I)
3988    Str.push_back(buffer[NDigits-I-1]);
3989}
3990
3991bool IEEEFloat::getExactInverse(APFloat *inv) const {
3992  // Special floats and denormals have no exact inverse.
3993  if (!isFiniteNonZero())
3994    return false;
3995
3996  // Check that the number is a power of two by making sure that only the
3997  // integer bit is set in the significand.
3998  if (significandLSB() != semantics->precision - 1)
3999    return false;
4000
4001  // Get the inverse.
4002  IEEEFloat reciprocal(*semantics, 1ULL);
4003  if (reciprocal.divide(*this, rmNearestTiesToEven) != opOK)
4004    return false;
4005
4006  // Avoid multiplication with a denormal, it is not safe on all platforms and
4007  // may be slower than a normal division.
4008  if (reciprocal.isDenormal())
4009    return false;
4010
4011  assert(reciprocal.isFiniteNonZero() &&
4012         reciprocal.significandLSB() == reciprocal.semantics->precision - 1);
4013
4014  if (inv)
4015    *inv = APFloat(reciprocal, *semantics);
4016
4017  return true;
4018}
4019
4020bool IEEEFloat::isSignaling() const {
4021  if (!isNaN())
4022    return false;
4023
4024  // IEEE-754R 2008 6.2.1: A signaling NaN bit string should be encoded with the
4025  // first bit of the trailing significand being 0.
4026  return !APInt::tcExtractBit(significandParts(), semantics->precision - 2);
4027}
4028
4029/// IEEE-754R 2008 5.3.1: nextUp/nextDown.
4030///
4031/// *NOTE* since nextDown(x) = -nextUp(-x), we only implement nextUp with
4032/// appropriate sign switching before/after the computation.
4033IEEEFloat::opStatus IEEEFloat::next(bool nextDown) {
4034  // If we are performing nextDown, swap sign so we have -x.
4035  if (nextDown)
4036    changeSign();
4037
4038  // Compute nextUp(x)
4039  opStatus result = opOK;
4040
4041  // Handle each float category separately.
4042  switch (category) {
4043  case fcInfinity:
4044    // nextUp(+inf) = +inf
4045    if (!isNegative())
4046      break;
4047    // nextUp(-inf) = -getLargest()
4048    makeLargest(true);
4049    break;
4050  case fcNaN:
4051    // IEEE-754R 2008 6.2 Par 2: nextUp(sNaN) = qNaN. Set Invalid flag.
4052    // IEEE-754R 2008 6.2: nextUp(qNaN) = qNaN. Must be identity so we do not
4053    //                     change the payload.
4054    if (isSignaling()) {
4055      result = opInvalidOp;
4056      // For consistency, propagate the sign of the sNaN to the qNaN.
4057      makeNaN(false, isNegative(), nullptr);
4058    }
4059    break;
4060  case fcZero:
4061    // nextUp(pm 0) = +getSmallest()
4062    makeSmallest(false);
4063    break;
4064  case fcNormal:
4065    // nextUp(-getSmallest()) = -0
4066    if (isSmallest() && isNegative()) {
4067      APInt::tcSet(significandParts(), 0, partCount());
4068      category = fcZero;
4069      exponent = 0;
4070      break;
4071    }
4072
4073    // nextUp(getLargest()) == INFINITY
4074    if (isLargest() && !isNegative()) {
4075      APInt::tcSet(significandParts(), 0, partCount());
4076      category = fcInfinity;
4077      exponent = semantics->maxExponent + 1;
4078      break;
4079    }
4080
4081    // nextUp(normal) == normal + inc.
4082    if (isNegative()) {
4083      // If we are negative, we need to decrement the significand.
4084
4085      // We only cross a binade boundary that requires adjusting the exponent
4086      // if:
4087      //   1. exponent != semantics->minExponent. This implies we are not in the
4088      //   smallest binade or are dealing with denormals.
4089      //   2. Our significand excluding the integral bit is all zeros.
4090      bool WillCrossBinadeBoundary =
4091        exponent != semantics->minExponent && isSignificandAllZeros();
4092
4093      // Decrement the significand.
4094      //
4095      // We always do this since:
4096      //   1. If we are dealing with a non-binade decrement, by definition we
4097      //   just decrement the significand.
4098      //   2. If we are dealing with a normal -> normal binade decrement, since
4099      //   we have an explicit integral bit the fact that all bits but the
4100      //   integral bit are zero implies that subtracting one will yield a
4101      //   significand with 0 integral bit and 1 in all other spots. Thus we
4102      //   must just adjust the exponent and set the integral bit to 1.
4103      //   3. If we are dealing with a normal -> denormal binade decrement,
4104      //   since we set the integral bit to 0 when we represent denormals, we
4105      //   just decrement the significand.
4106      integerPart *Parts = significandParts();
4107      APInt::tcDecrement(Parts, partCount());
4108
4109      if (WillCrossBinadeBoundary) {
4110        // Our result is a normal number. Do the following:
4111        // 1. Set the integral bit to 1.
4112        // 2. Decrement the exponent.
4113        APInt::tcSetBit(Parts, semantics->precision - 1);
4114        exponent--;
4115      }
4116    } else {
4117      // If we are positive, we need to increment the significand.
4118
4119      // We only cross a binade boundary that requires adjusting the exponent if
4120      // the input is not a denormal and all of said input's significand bits
4121      // are set. If all of said conditions are true: clear the significand, set
4122      // the integral bit to 1, and increment the exponent. If we have a
4123      // denormal always increment since moving denormals and the numbers in the
4124      // smallest normal binade have the same exponent in our representation.
4125      bool WillCrossBinadeBoundary = !isDenormal() && isSignificandAllOnes();
4126
4127      if (WillCrossBinadeBoundary) {
4128        integerPart *Parts = significandParts();
4129        APInt::tcSet(Parts, 0, partCount());
4130        APInt::tcSetBit(Parts, semantics->precision - 1);
4131        assert(exponent != semantics->maxExponent &&
4132               "We can not increment an exponent beyond the maxExponent allowed"
4133               " by the given floating point semantics.");
4134        exponent++;
4135      } else {
4136        incrementSignificand();
4137      }
4138    }
4139    break;
4140  }
4141
4142  // If we are performing nextDown, swap sign so we have -nextUp(-x)
4143  if (nextDown)
4144    changeSign();
4145
4146  return result;
4147}
4148
4149void IEEEFloat::makeInf(bool Negative) {
4150  category = fcInfinity;
4151  sign = Negative;
4152  exponent = semantics->maxExponent + 1;
4153  APInt::tcSet(significandParts(), 0, partCount());
4154}
4155
4156void IEEEFloat::makeZero(bool Negative) {
4157  category = fcZero;
4158  sign = Negative;
4159  exponent = semantics->minExponent-1;
4160  APInt::tcSet(significandParts(), 0, partCount());
4161}
4162
4163void IEEEFloat::makeQuiet() {
4164  assert(isNaN());
4165  APInt::tcSetBit(significandParts(), semantics->precision - 2);
4166}
4167
4168int ilogb(const IEEEFloat &Arg) {
4169  if (Arg.isNaN())
4170    return IEEEFloat::IEK_NaN;
4171  if (Arg.isZero())
4172    return IEEEFloat::IEK_Zero;
4173  if (Arg.isInfinity())
4174    return IEEEFloat::IEK_Inf;
4175  if (!Arg.isDenormal())
4176    return Arg.exponent;
4177
4178  IEEEFloat Normalized(Arg);
4179  int SignificandBits = Arg.getSemantics().precision - 1;
4180
4181  Normalized.exponent += SignificandBits;
4182  Normalized.normalize(IEEEFloat::rmNearestTiesToEven, lfExactlyZero);
4183  return Normalized.exponent - SignificandBits;
4184}
4185
4186IEEEFloat scalbn(IEEEFloat X, int Exp, IEEEFloat::roundingMode RoundingMode) {
4187  auto MaxExp = X.getSemantics().maxExponent;
4188  auto MinExp = X.getSemantics().minExponent;
4189
4190  // If Exp is wildly out-of-scale, simply adding it to X.exponent will
4191  // overflow; clamp it to a safe range before adding, but ensure that the range
4192  // is large enough that the clamp does not change the result. The range we
4193  // need to support is the difference between the largest possible exponent and
4194  // the normalized exponent of half the smallest denormal.
4195
4196  int SignificandBits = X.getSemantics().precision - 1;
4197  int MaxIncrement = MaxExp - (MinExp - SignificandBits) + 1;
4198
4199  // Clamp to one past the range ends to let normalize handle overlflow.
4200  X.exponent += std::min(std::max(Exp, -MaxIncrement - 1), MaxIncrement);
4201  X.normalize(RoundingMode, lfExactlyZero);
4202  if (X.isNaN())
4203    X.makeQuiet();
4204  return X;
4205}
4206
4207IEEEFloat frexp(const IEEEFloat &Val, int &Exp, IEEEFloat::roundingMode RM) {
4208  Exp = ilogb(Val);
4209
4210  // Quiet signalling nans.
4211  if (Exp == IEEEFloat::IEK_NaN) {
4212    IEEEFloat Quiet(Val);
4213    Quiet.makeQuiet();
4214    return Quiet;
4215  }
4216
4217  if (Exp == IEEEFloat::IEK_Inf)
4218    return Val;
4219
4220  // 1 is added because frexp is defined to return a normalized fraction in
4221  // +/-[0.5, 1.0), rather than the usual +/-[1.0, 2.0).
4222  Exp = Exp == IEEEFloat::IEK_Zero ? 0 : Exp + 1;
4223  return scalbn(Val, -Exp, RM);
4224}
4225
4226DoubleAPFloat::DoubleAPFloat(const fltSemantics &S)
4227    : Semantics(&S),
4228      Floats(new APFloat[2]{APFloat(semIEEEdouble), APFloat(semIEEEdouble)}) {
4229  assert(Semantics == &semPPCDoubleDouble);
4230}
4231
4232DoubleAPFloat::DoubleAPFloat(const fltSemantics &S, uninitializedTag)
4233    : Semantics(&S),
4234      Floats(new APFloat[2]{APFloat(semIEEEdouble, uninitialized),
4235                            APFloat(semIEEEdouble, uninitialized)}) {
4236  assert(Semantics == &semPPCDoubleDouble);
4237}
4238
4239DoubleAPFloat::DoubleAPFloat(const fltSemantics &S, integerPart I)
4240    : Semantics(&S), Floats(new APFloat[2]{APFloat(semIEEEdouble, I),
4241                                           APFloat(semIEEEdouble)}) {
4242  assert(Semantics == &semPPCDoubleDouble);
4243}
4244
4245DoubleAPFloat::DoubleAPFloat(const fltSemantics &S, const APInt &I)
4246    : Semantics(&S),
4247      Floats(new APFloat[2]{
4248          APFloat(semIEEEdouble, APInt(64, I.getRawData()[0])),
4249          APFloat(semIEEEdouble, APInt(64, I.getRawData()[1]))}) {
4250  assert(Semantics == &semPPCDoubleDouble);
4251}
4252
4253DoubleAPFloat::DoubleAPFloat(const fltSemantics &S, APFloat &&First,
4254                             APFloat &&Second)
4255    : Semantics(&S),
4256      Floats(new APFloat[2]{std::move(First), std::move(Second)}) {
4257  assert(Semantics == &semPPCDoubleDouble);
4258  assert(&Floats[0].getSemantics() == &semIEEEdouble);
4259  assert(&Floats[1].getSemantics() == &semIEEEdouble);
4260}
4261
4262DoubleAPFloat::DoubleAPFloat(const DoubleAPFloat &RHS)
4263    : Semantics(RHS.Semantics),
4264      Floats(RHS.Floats ? new APFloat[2]{APFloat(RHS.Floats[0]),
4265                                         APFloat(RHS.Floats[1])}
4266                        : nullptr) {
4267  assert(Semantics == &semPPCDoubleDouble);
4268}
4269
4270DoubleAPFloat::DoubleAPFloat(DoubleAPFloat &&RHS)
4271    : Semantics(RHS.Semantics), Floats(std::move(RHS.Floats)) {
4272  RHS.Semantics = &semBogus;
4273  assert(Semantics == &semPPCDoubleDouble);
4274}
4275
4276DoubleAPFloat &DoubleAPFloat::operator=(const DoubleAPFloat &RHS) {
4277  if (Semantics == RHS.Semantics && RHS.Floats) {
4278    Floats[0] = RHS.Floats[0];
4279    Floats[1] = RHS.Floats[1];
4280  } else if (this != &RHS) {
4281    this->~DoubleAPFloat();
4282    new (this) DoubleAPFloat(RHS);
4283  }
4284  return *this;
4285}
4286
4287// Implement addition, subtraction, multiplication and division based on:
4288// "Software for Doubled-Precision Floating-Point Computations",
4289// by Seppo Linnainmaa, ACM TOMS vol 7 no 3, September 1981, pages 272-283.
4290APFloat::opStatus DoubleAPFloat::addImpl(const APFloat &a, const APFloat &aa,
4291                                         const APFloat &c, const APFloat &cc,
4292                                         roundingMode RM) {
4293  int Status = opOK;
4294  APFloat z = a;
4295  Status |= z.add(c, RM);
4296  if (!z.isFinite()) {
4297    if (!z.isInfinity()) {
4298      Floats[0] = std::move(z);
4299      Floats[1].makeZero(/* Neg = */ false);
4300      return (opStatus)Status;
4301    }
4302    Status = opOK;
4303    auto AComparedToC = a.compareAbsoluteValue(c);
4304    z = cc;
4305    Status |= z.add(aa, RM);
4306    if (AComparedToC == APFloat::cmpGreaterThan) {
4307      // z = cc + aa + c + a;
4308      Status |= z.add(c, RM);
4309      Status |= z.add(a, RM);
4310    } else {
4311      // z = cc + aa + a + c;
4312      Status |= z.add(a, RM);
4313      Status |= z.add(c, RM);
4314    }
4315    if (!z.isFinite()) {
4316      Floats[0] = std::move(z);
4317      Floats[1].makeZero(/* Neg = */ false);
4318      return (opStatus)Status;
4319    }
4320    Floats[0] = z;
4321    APFloat zz = aa;
4322    Status |= zz.add(cc, RM);
4323    if (AComparedToC == APFloat::cmpGreaterThan) {
4324      // Floats[1] = a - z + c + zz;
4325      Floats[1] = a;
4326      Status |= Floats[1].subtract(z, RM);
4327      Status |= Floats[1].add(c, RM);
4328      Status |= Floats[1].add(zz, RM);
4329    } else {
4330      // Floats[1] = c - z + a + zz;
4331      Floats[1] = c;
4332      Status |= Floats[1].subtract(z, RM);
4333      Status |= Floats[1].add(a, RM);
4334      Status |= Floats[1].add(zz, RM);
4335    }
4336  } else {
4337    // q = a - z;
4338    APFloat q = a;
4339    Status |= q.subtract(z, RM);
4340
4341    // zz = q + c + (a - (q + z)) + aa + cc;
4342    // Compute a - (q + z) as -((q + z) - a) to avoid temporary copies.
4343    auto zz = q;
4344    Status |= zz.add(c, RM);
4345    Status |= q.add(z, RM);
4346    Status |= q.subtract(a, RM);
4347    q.changeSign();
4348    Status |= zz.add(q, RM);
4349    Status |= zz.add(aa, RM);
4350    Status |= zz.add(cc, RM);
4351    if (zz.isZero() && !zz.isNegative()) {
4352      Floats[0] = std::move(z);
4353      Floats[1].makeZero(/* Neg = */ false);
4354      return opOK;
4355    }
4356    Floats[0] = z;
4357    Status |= Floats[0].add(zz, RM);
4358    if (!Floats[0].isFinite()) {
4359      Floats[1].makeZero(/* Neg = */ false);
4360      return (opStatus)Status;
4361    }
4362    Floats[1] = std::move(z);
4363    Status |= Floats[1].subtract(Floats[0], RM);
4364    Status |= Floats[1].add(zz, RM);
4365  }
4366  return (opStatus)Status;
4367}
4368
4369APFloat::opStatus DoubleAPFloat::addWithSpecial(const DoubleAPFloat &LHS,
4370                                                const DoubleAPFloat &RHS,
4371                                                DoubleAPFloat &Out,
4372                                                roundingMode RM) {
4373  if (LHS.getCategory() == fcNaN) {
4374    Out = LHS;
4375    return opOK;
4376  }
4377  if (RHS.getCategory() == fcNaN) {
4378    Out = RHS;
4379    return opOK;
4380  }
4381  if (LHS.getCategory() == fcZero) {
4382    Out = RHS;
4383    return opOK;
4384  }
4385  if (RHS.getCategory() == fcZero) {
4386    Out = LHS;
4387    return opOK;
4388  }
4389  if (LHS.getCategory() == fcInfinity && RHS.getCategory() == fcInfinity &&
4390      LHS.isNegative() != RHS.isNegative()) {
4391    Out.makeNaN(false, Out.isNegative(), nullptr);
4392    return opInvalidOp;
4393  }
4394  if (LHS.getCategory() == fcInfinity) {
4395    Out = LHS;
4396    return opOK;
4397  }
4398  if (RHS.getCategory() == fcInfinity) {
4399    Out = RHS;
4400    return opOK;
4401  }
4402  assert(LHS.getCategory() == fcNormal && RHS.getCategory() == fcNormal);
4403
4404  APFloat A(LHS.Floats[0]), AA(LHS.Floats[1]), C(RHS.Floats[0]),
4405      CC(RHS.Floats[1]);
4406  assert(&A.getSemantics() == &semIEEEdouble);
4407  assert(&AA.getSemantics() == &semIEEEdouble);
4408  assert(&C.getSemantics() == &semIEEEdouble);
4409  assert(&CC.getSemantics() == &semIEEEdouble);
4410  assert(&Out.Floats[0].getSemantics() == &semIEEEdouble);
4411  assert(&Out.Floats[1].getSemantics() == &semIEEEdouble);
4412  return Out.addImpl(A, AA, C, CC, RM);
4413}
4414
4415APFloat::opStatus DoubleAPFloat::add(const DoubleAPFloat &RHS,
4416                                     roundingMode RM) {
4417  return addWithSpecial(*this, RHS, *this, RM);
4418}
4419
4420APFloat::opStatus DoubleAPFloat::subtract(const DoubleAPFloat &RHS,
4421                                          roundingMode RM) {
4422  changeSign();
4423  auto Ret = add(RHS, RM);
4424  changeSign();
4425  return Ret;
4426}
4427
4428APFloat::opStatus DoubleAPFloat::multiply(const DoubleAPFloat &RHS,
4429                                          APFloat::roundingMode RM) {
4430  const auto &LHS = *this;
4431  auto &Out = *this;
4432  /* Interesting observation: For special categories, finding the lowest
4433     common ancestor of the following layered graph gives the correct
4434     return category:
4435
4436        NaN
4437       /   \
4438     Zero  Inf
4439       \   /
4440       Normal
4441
4442     e.g. NaN * NaN = NaN
4443          Zero * Inf = NaN
4444          Normal * Zero = Zero
4445          Normal * Inf = Inf
4446  */
4447  if (LHS.getCategory() == fcNaN) {
4448    Out = LHS;
4449    return opOK;
4450  }
4451  if (RHS.getCategory() == fcNaN) {
4452    Out = RHS;
4453    return opOK;
4454  }
4455  if ((LHS.getCategory() == fcZero && RHS.getCategory() == fcInfinity) ||
4456      (LHS.getCategory() == fcInfinity && RHS.getCategory() == fcZero)) {
4457    Out.makeNaN(false, false, nullptr);
4458    return opOK;
4459  }
4460  if (LHS.getCategory() == fcZero || LHS.getCategory() == fcInfinity) {
4461    Out = LHS;
4462    return opOK;
4463  }
4464  if (RHS.getCategory() == fcZero || RHS.getCategory() == fcInfinity) {
4465    Out = RHS;
4466    return opOK;
4467  }
4468  assert(LHS.getCategory() == fcNormal && RHS.getCategory() == fcNormal &&
4469         "Special cases not handled exhaustively");
4470
4471  int Status = opOK;
4472  APFloat A = Floats[0], B = Floats[1], C = RHS.Floats[0], D = RHS.Floats[1];
4473  // t = a * c
4474  APFloat T = A;
4475  Status |= T.multiply(C, RM);
4476  if (!T.isFiniteNonZero()) {
4477    Floats[0] = T;
4478    Floats[1].makeZero(/* Neg = */ false);
4479    return (opStatus)Status;
4480  }
4481
4482  // tau = fmsub(a, c, t), that is -fmadd(-a, c, t).
4483  APFloat Tau = A;
4484  T.changeSign();
4485  Status |= Tau.fusedMultiplyAdd(C, T, RM);
4486  T.changeSign();
4487  {
4488    // v = a * d
4489    APFloat V = A;
4490    Status |= V.multiply(D, RM);
4491    // w = b * c
4492    APFloat W = B;
4493    Status |= W.multiply(C, RM);
4494    Status |= V.add(W, RM);
4495    // tau += v + w
4496    Status |= Tau.add(V, RM);
4497  }
4498  // u = t + tau
4499  APFloat U = T;
4500  Status |= U.add(Tau, RM);
4501
4502  Floats[0] = U;
4503  if (!U.isFinite()) {
4504    Floats[1].makeZero(/* Neg = */ false);
4505  } else {
4506    // Floats[1] = (t - u) + tau
4507    Status |= T.subtract(U, RM);
4508    Status |= T.add(Tau, RM);
4509    Floats[1] = T;
4510  }
4511  return (opStatus)Status;
4512}
4513
4514APFloat::opStatus DoubleAPFloat::divide(const DoubleAPFloat &RHS,
4515                                        APFloat::roundingMode RM) {
4516  assert(Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
4517  APFloat Tmp(semPPCDoubleDoubleLegacy, bitcastToAPInt());
4518  auto Ret =
4519      Tmp.divide(APFloat(semPPCDoubleDoubleLegacy, RHS.bitcastToAPInt()), RM);
4520  *this = DoubleAPFloat(semPPCDoubleDouble, Tmp.bitcastToAPInt());
4521  return Ret;
4522}
4523
4524APFloat::opStatus DoubleAPFloat::remainder(const DoubleAPFloat &RHS) {
4525  assert(Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
4526  APFloat Tmp(semPPCDoubleDoubleLegacy, bitcastToAPInt());
4527  auto Ret =
4528      Tmp.remainder(APFloat(semPPCDoubleDoubleLegacy, RHS.bitcastToAPInt()));
4529  *this = DoubleAPFloat(semPPCDoubleDouble, Tmp.bitcastToAPInt());
4530  return Ret;
4531}
4532
4533APFloat::opStatus DoubleAPFloat::mod(const DoubleAPFloat &RHS) {
4534  assert(Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
4535  APFloat Tmp(semPPCDoubleDoubleLegacy, bitcastToAPInt());
4536  auto Ret = Tmp.mod(APFloat(semPPCDoubleDoubleLegacy, RHS.bitcastToAPInt()));
4537  *this = DoubleAPFloat(semPPCDoubleDouble, Tmp.bitcastToAPInt());
4538  return Ret;
4539}
4540
4541APFloat::opStatus
4542DoubleAPFloat::fusedMultiplyAdd(const DoubleAPFloat &Multiplicand,
4543                                const DoubleAPFloat &Addend,
4544                                APFloat::roundingMode RM) {
4545  assert(Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
4546  APFloat Tmp(semPPCDoubleDoubleLegacy, bitcastToAPInt());
4547  auto Ret = Tmp.fusedMultiplyAdd(
4548      APFloat(semPPCDoubleDoubleLegacy, Multiplicand.bitcastToAPInt()),
4549      APFloat(semPPCDoubleDoubleLegacy, Addend.bitcastToAPInt()), RM);
4550  *this = DoubleAPFloat(semPPCDoubleDouble, Tmp.bitcastToAPInt());
4551  return Ret;
4552}
4553
4554APFloat::opStatus DoubleAPFloat::roundToIntegral(APFloat::roundingMode RM) {
4555  assert(Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
4556  APFloat Tmp(semPPCDoubleDoubleLegacy, bitcastToAPInt());
4557  auto Ret = Tmp.roundToIntegral(RM);
4558  *this = DoubleAPFloat(semPPCDoubleDouble, Tmp.bitcastToAPInt());
4559  return Ret;
4560}
4561
4562void DoubleAPFloat::changeSign() {
4563  Floats[0].changeSign();
4564  Floats[1].changeSign();
4565}
4566
4567APFloat::cmpResult
4568DoubleAPFloat::compareAbsoluteValue(const DoubleAPFloat &RHS) const {
4569  auto Result = Floats[0].compareAbsoluteValue(RHS.Floats[0]);
4570  if (Result != cmpEqual)
4571    return Result;
4572  Result = Floats[1].compareAbsoluteValue(RHS.Floats[1]);
4573  if (Result == cmpLessThan || Result == cmpGreaterThan) {
4574    auto Against = Floats[0].isNegative() ^ Floats[1].isNegative();
4575    auto RHSAgainst = RHS.Floats[0].isNegative() ^ RHS.Floats[1].isNegative();
4576    if (Against && !RHSAgainst)
4577      return cmpLessThan;
4578    if (!Against && RHSAgainst)
4579      return cmpGreaterThan;
4580    if (!Against && !RHSAgainst)
4581      return Result;
4582    if (Against && RHSAgainst)
4583      return (cmpResult)(cmpLessThan + cmpGreaterThan - Result);
4584  }
4585  return Result;
4586}
4587
4588APFloat::fltCategory DoubleAPFloat::getCategory() const {
4589  return Floats[0].getCategory();
4590}
4591
4592bool DoubleAPFloat::isNegative() const { return Floats[0].isNegative(); }
4593
4594void DoubleAPFloat::makeInf(bool Neg) {
4595  Floats[0].makeInf(Neg);
4596  Floats[1].makeZero(/* Neg = */ false);
4597}
4598
4599void DoubleAPFloat::makeZero(bool Neg) {
4600  Floats[0].makeZero(Neg);
4601  Floats[1].makeZero(/* Neg = */ false);
4602}
4603
4604void DoubleAPFloat::makeLargest(bool Neg) {
4605  assert(Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
4606  Floats[0] = APFloat(semIEEEdouble, APInt(64, 0x7fefffffffffffffull));
4607  Floats[1] = APFloat(semIEEEdouble, APInt(64, 0x7c8ffffffffffffeull));
4608  if (Neg)
4609    changeSign();
4610}
4611
4612void DoubleAPFloat::makeSmallest(bool Neg) {
4613  assert(Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
4614  Floats[0].makeSmallest(Neg);
4615  Floats[1].makeZero(/* Neg = */ false);
4616}
4617
4618void DoubleAPFloat::makeSmallestNormalized(bool Neg) {
4619  assert(Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
4620  Floats[0] = APFloat(semIEEEdouble, APInt(64, 0x0360000000000000ull));
4621  if (Neg)
4622    Floats[0].changeSign();
4623  Floats[1].makeZero(/* Neg = */ false);
4624}
4625
4626void DoubleAPFloat::makeNaN(bool SNaN, bool Neg, const APInt *fill) {
4627  Floats[0].makeNaN(SNaN, Neg, fill);
4628  Floats[1].makeZero(/* Neg = */ false);
4629}
4630
4631APFloat::cmpResult DoubleAPFloat::compare(const DoubleAPFloat &RHS) const {
4632  auto Result = Floats[0].compare(RHS.Floats[0]);
4633  // |Float[0]| > |Float[1]|
4634  if (Result == APFloat::cmpEqual)
4635    return Floats[1].compare(RHS.Floats[1]);
4636  return Result;
4637}
4638
4639bool DoubleAPFloat::bitwiseIsEqual(const DoubleAPFloat &RHS) const {
4640  return Floats[0].bitwiseIsEqual(RHS.Floats[0]) &&
4641         Floats[1].bitwiseIsEqual(RHS.Floats[1]);
4642}
4643
4644hash_code hash_value(const DoubleAPFloat &Arg) {
4645  if (Arg.Floats)
4646    return hash_combine(hash_value(Arg.Floats[0]), hash_value(Arg.Floats[1]));
4647  return hash_combine(Arg.Semantics);
4648}
4649
4650APInt DoubleAPFloat::bitcastToAPInt() const {
4651  assert(Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
4652  uint64_t Data[] = {
4653      Floats[0].bitcastToAPInt().getRawData()[0],
4654      Floats[1].bitcastToAPInt().getRawData()[0],
4655  };
4656  return APInt(128, 2, Data);
4657}
4658
4659Expected<APFloat::opStatus> DoubleAPFloat::convertFromString(StringRef S,
4660                                                             roundingMode RM) {
4661  assert(Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
4662  APFloat Tmp(semPPCDoubleDoubleLegacy);
4663  auto Ret = Tmp.convertFromString(S, RM);
4664  *this = DoubleAPFloat(semPPCDoubleDouble, Tmp.bitcastToAPInt());
4665  return Ret;
4666}
4667
4668APFloat::opStatus DoubleAPFloat::next(bool nextDown) {
4669  assert(Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
4670  APFloat Tmp(semPPCDoubleDoubleLegacy, bitcastToAPInt());
4671  auto Ret = Tmp.next(nextDown);
4672  *this = DoubleAPFloat(semPPCDoubleDouble, Tmp.bitcastToAPInt());
4673  return Ret;
4674}
4675
4676APFloat::opStatus
4677DoubleAPFloat::convertToInteger(MutableArrayRef<integerPart> Input,
4678                                unsigned int Width, bool IsSigned,
4679                                roundingMode RM, bool *IsExact) const {
4680  assert(Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
4681  return APFloat(semPPCDoubleDoubleLegacy, bitcastToAPInt())
4682      .convertToInteger(Input, Width, IsSigned, RM, IsExact);
4683}
4684
4685APFloat::opStatus DoubleAPFloat::convertFromAPInt(const APInt &Input,
4686                                                  bool IsSigned,
4687                                                  roundingMode RM) {
4688  assert(Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
4689  APFloat Tmp(semPPCDoubleDoubleLegacy);
4690  auto Ret = Tmp.convertFromAPInt(Input, IsSigned, RM);
4691  *this = DoubleAPFloat(semPPCDoubleDouble, Tmp.bitcastToAPInt());
4692  return Ret;
4693}
4694
4695APFloat::opStatus
4696DoubleAPFloat::convertFromSignExtendedInteger(const integerPart *Input,
4697                                              unsigned int InputSize,
4698                                              bool IsSigned, roundingMode RM) {
4699  assert(Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
4700  APFloat Tmp(semPPCDoubleDoubleLegacy);
4701  auto Ret = Tmp.convertFromSignExtendedInteger(Input, InputSize, IsSigned, RM);
4702  *this = DoubleAPFloat(semPPCDoubleDouble, Tmp.bitcastToAPInt());
4703  return Ret;
4704}
4705
4706APFloat::opStatus
4707DoubleAPFloat::convertFromZeroExtendedInteger(const integerPart *Input,
4708                                              unsigned int InputSize,
4709                                              bool IsSigned, roundingMode RM) {
4710  assert(Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
4711  APFloat Tmp(semPPCDoubleDoubleLegacy);
4712  auto Ret = Tmp.convertFromZeroExtendedInteger(Input, InputSize, IsSigned, RM);
4713  *this = DoubleAPFloat(semPPCDoubleDouble, Tmp.bitcastToAPInt());
4714  return Ret;
4715}
4716
4717unsigned int DoubleAPFloat::convertToHexString(char *DST,
4718                                               unsigned int HexDigits,
4719                                               bool UpperCase,
4720                                               roundingMode RM) const {
4721  assert(Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
4722  return APFloat(semPPCDoubleDoubleLegacy, bitcastToAPInt())
4723      .convertToHexString(DST, HexDigits, UpperCase, RM);
4724}
4725
4726bool DoubleAPFloat::isDenormal() const {
4727  return getCategory() == fcNormal &&
4728         (Floats[0].isDenormal() || Floats[1].isDenormal() ||
4729          // (double)(Hi + Lo) == Hi defines a normal number.
4730          Floats[0] != Floats[0] + Floats[1]);
4731}
4732
4733bool DoubleAPFloat::isSmallest() const {
4734  if (getCategory() != fcNormal)
4735    return false;
4736  DoubleAPFloat Tmp(*this);
4737  Tmp.makeSmallest(this->isNegative());
4738  return Tmp.compare(*this) == cmpEqual;
4739}
4740
4741bool DoubleAPFloat::isLargest() const {
4742  if (getCategory() != fcNormal)
4743    return false;
4744  DoubleAPFloat Tmp(*this);
4745  Tmp.makeLargest(this->isNegative());
4746  return Tmp.compare(*this) == cmpEqual;
4747}
4748
4749bool DoubleAPFloat::isInteger() const {
4750  assert(Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
4751  return Floats[0].isInteger() && Floats[1].isInteger();
4752}
4753
4754void DoubleAPFloat::toString(SmallVectorImpl<char> &Str,
4755                             unsigned FormatPrecision,
4756                             unsigned FormatMaxPadding,
4757                             bool TruncateZero) const {
4758  assert(Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
4759  APFloat(semPPCDoubleDoubleLegacy, bitcastToAPInt())
4760      .toString(Str, FormatPrecision, FormatMaxPadding, TruncateZero);
4761}
4762
4763bool DoubleAPFloat::getExactInverse(APFloat *inv) const {
4764  assert(Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
4765  APFloat Tmp(semPPCDoubleDoubleLegacy, bitcastToAPInt());
4766  if (!inv)
4767    return Tmp.getExactInverse(nullptr);
4768  APFloat Inv(semPPCDoubleDoubleLegacy);
4769  auto Ret = Tmp.getExactInverse(&Inv);
4770  *inv = APFloat(semPPCDoubleDouble, Inv.bitcastToAPInt());
4771  return Ret;
4772}
4773
4774DoubleAPFloat scalbn(DoubleAPFloat Arg, int Exp, APFloat::roundingMode RM) {
4775  assert(Arg.Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
4776  return DoubleAPFloat(semPPCDoubleDouble, scalbn(Arg.Floats[0], Exp, RM),
4777                       scalbn(Arg.Floats[1], Exp, RM));
4778}
4779
4780DoubleAPFloat frexp(const DoubleAPFloat &Arg, int &Exp,
4781                    APFloat::roundingMode RM) {
4782  assert(Arg.Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
4783  APFloat First = frexp(Arg.Floats[0], Exp, RM);
4784  APFloat Second = Arg.Floats[1];
4785  if (Arg.getCategory() == APFloat::fcNormal)
4786    Second = scalbn(Second, -Exp, RM);
4787  return DoubleAPFloat(semPPCDoubleDouble, std::move(First), std::move(Second));
4788}
4789
4790} // End detail namespace
4791
4792APFloat::Storage::Storage(IEEEFloat F, const fltSemantics &Semantics) {
4793  if (usesLayout<IEEEFloat>(Semantics)) {
4794    new (&IEEE) IEEEFloat(std::move(F));
4795    return;
4796  }
4797  if (usesLayout<DoubleAPFloat>(Semantics)) {
4798    const fltSemantics& S = F.getSemantics();
4799    new (&Double)
4800        DoubleAPFloat(Semantics, APFloat(std::move(F), S),
4801                      APFloat(semIEEEdouble));
4802    return;
4803  }
4804  llvm_unreachable("Unexpected semantics");
4805}
4806
4807Expected<APFloat::opStatus> APFloat::convertFromString(StringRef Str,
4808                                                       roundingMode RM) {
4809  APFLOAT_DISPATCH_ON_SEMANTICS(convertFromString(Str, RM));
4810}
4811
4812hash_code hash_value(const APFloat &Arg) {
4813  if (APFloat::usesLayout<detail::IEEEFloat>(Arg.getSemantics()))
4814    return hash_value(Arg.U.IEEE);
4815  if (APFloat::usesLayout<detail::DoubleAPFloat>(Arg.getSemantics()))
4816    return hash_value(Arg.U.Double);
4817  llvm_unreachable("Unexpected semantics");
4818}
4819
4820APFloat::APFloat(const fltSemantics &Semantics, StringRef S)
4821    : APFloat(Semantics) {
4822  auto StatusOrErr = convertFromString(S, rmNearestTiesToEven);
4823  assert(StatusOrErr && "Invalid floating point representation");
4824  consumeError(StatusOrErr.takeError());
4825}
4826
4827APFloat::opStatus APFloat::convert(const fltSemantics &ToSemantics,
4828                                   roundingMode RM, bool *losesInfo) {
4829  if (&getSemantics() == &ToSemantics) {
4830    *losesInfo = false;
4831    return opOK;
4832  }
4833  if (usesLayout<IEEEFloat>(getSemantics()) &&
4834      usesLayout<IEEEFloat>(ToSemantics))
4835    return U.IEEE.convert(ToSemantics, RM, losesInfo);
4836  if (usesLayout<IEEEFloat>(getSemantics()) &&
4837      usesLayout<DoubleAPFloat>(ToSemantics)) {
4838    assert(&ToSemantics == &semPPCDoubleDouble);
4839    auto Ret = U.IEEE.convert(semPPCDoubleDoubleLegacy, RM, losesInfo);
4840    *this = APFloat(ToSemantics, U.IEEE.bitcastToAPInt());
4841    return Ret;
4842  }
4843  if (usesLayout<DoubleAPFloat>(getSemantics()) &&
4844      usesLayout<IEEEFloat>(ToSemantics)) {
4845    auto Ret = getIEEE().convert(ToSemantics, RM, losesInfo);
4846    *this = APFloat(std::move(getIEEE()), ToSemantics);
4847    return Ret;
4848  }
4849  llvm_unreachable("Unexpected semantics");
4850}
4851
4852APFloat APFloat::getAllOnesValue(const fltSemantics &Semantics,
4853                                 unsigned BitWidth) {
4854  return APFloat(Semantics, APInt::getAllOnesValue(BitWidth));
4855}
4856
4857void APFloat::print(raw_ostream &OS) const {
4858  SmallVector<char, 16> Buffer;
4859  toString(Buffer);
4860  OS << Buffer << "\n";
4861}
4862
4863#if !defined(NDEBUG) || defined(LLVM_ENABLE_DUMP)
4864LLVM_DUMP_METHOD void APFloat::dump() const { print(dbgs()); }
4865#endif
4866
4867void APFloat::Profile(FoldingSetNodeID &NID) const {
4868  NID.Add(bitcastToAPInt());
4869}
4870
4871/* Same as convertToInteger(integerPart*, ...), except the result is returned in
4872   an APSInt, whose initial bit-width and signed-ness are used to determine the
4873   precision of the conversion.
4874 */
4875APFloat::opStatus APFloat::convertToInteger(APSInt &result,
4876                                            roundingMode rounding_mode,
4877                                            bool *isExact) const {
4878  unsigned bitWidth = result.getBitWidth();
4879  SmallVector<uint64_t, 4> parts(result.getNumWords());
4880  opStatus status = convertToInteger(parts, bitWidth, result.isSigned(),
4881                                     rounding_mode, isExact);
4882  // Keeps the original signed-ness.
4883  result = APInt(bitWidth, parts);
4884  return status;
4885}
4886
4887} // End llvm namespace
4888
4889#undef APFLOAT_DISPATCH_ON_SEMANTICS
4890