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