1// Copyright 2018 Ulf Adams 2// 3// The contents of this file may be used under the terms of the Apache License, 4// Version 2.0. 5// 6// (See accompanying file LICENSE-Apache or copy at 7// http://www.apache.org/licenses/LICENSE-2.0) 8// 9// Alternatively, the contents of this file may be used under the terms of 10// the Boost Software License, Version 1.0. 11// (See accompanying file LICENSE-Boost or copy at 12// https://www.boost.org/LICENSE_1_0.txt) 13// 14// Unless required by applicable law or agreed to in writing, this software 15// is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY 16// KIND, either express or implied. 17 18// Runtime compiler options: 19// -DRYU_DEBUG Generate verbose debugging output to stdout. 20 21 22 23 24#ifdef RYU_DEBUG 25static char* s(uint128_t v) { 26 int len = decimalLength(v); 27 char* b = (char*) malloc((len + 1) * sizeof(char)); 28 for (int i = 0; i < len; i++) { 29 const uint32_t c = (uint32_t) (v % 10); 30 v /= 10; 31 b[len - 1 - i] = (char) ('0' + c); 32 } 33 b[len] = 0; 34 return b; 35} 36#endif 37 38#define ONE ((uint128_t) 1) 39 40struct floating_decimal_128 generic_binary_to_decimal( 41 const uint128_t ieeeMantissa, const uint32_t ieeeExponent, const bool ieeeSign, 42 const uint32_t mantissaBits, const uint32_t exponentBits, const bool explicitLeadingBit) { 43#ifdef RYU_DEBUG 44 printf("IN="); 45 for (int32_t bit = 127; bit >= 0; --bit) { 46 printf("%u", (uint32_t) ((bits >> bit) & 1)); 47 } 48 printf("\n"); 49#endif 50 51 const uint32_t bias = (1u << (exponentBits - 1)) - 1; 52 53 if (ieeeExponent == 0 && ieeeMantissa == 0) { 54 struct floating_decimal_128 fd; 55 fd.mantissa = 0; 56 fd.exponent = 0; 57 fd.sign = ieeeSign; 58 return fd; 59 } 60 if (ieeeExponent == ((1u << exponentBits) - 1u)) { 61 struct floating_decimal_128 fd; 62 fd.mantissa = explicitLeadingBit ? ieeeMantissa & ((ONE << (mantissaBits - 1)) - 1) : ieeeMantissa; 63 fd.exponent = FD128_EXCEPTIONAL_EXPONENT; 64 fd.sign = ieeeSign; 65 return fd; 66 } 67 68 int32_t e2; 69 uint128_t m2; 70 // We subtract 2 in all cases so that the bounds computation has 2 additional bits. 71 if (explicitLeadingBit) { 72 // mantissaBits includes the explicit leading bit, so we need to correct for that here. 73 if (ieeeExponent == 0) { 74 e2 = 1 - bias - mantissaBits + 1 - 2; 75 } else { 76 e2 = ieeeExponent - bias - mantissaBits + 1 - 2; 77 } 78 m2 = ieeeMantissa; 79 } else { 80 if (ieeeExponent == 0) { 81 e2 = 1 - bias - mantissaBits - 2; 82 m2 = ieeeMantissa; 83 } else { 84 e2 = ieeeExponent - bias - mantissaBits - 2; 85 m2 = (ONE << mantissaBits) | ieeeMantissa; 86 } 87 } 88 const bool even = (m2 & 1) == 0; 89 const bool acceptBounds = even; 90 91#ifdef RYU_DEBUG 92 printf("-> %s %s * 2^%d\n", ieeeSign ? "-" : "+", s(m2), e2 + 2); 93#endif 94 95 // Step 2: Determine the interval of legal decimal representations. 96 const uint128_t mv = 4 * m2; 97 // Implicit bool -> int conversion. True is 1, false is 0. 98 const uint32_t mmShift = 99 (ieeeMantissa != (explicitLeadingBit ? ONE << (mantissaBits - 1) : 0)) 100 || (ieeeExponent == 0); 101 102 // Step 3: Convert to a decimal power base using 128-bit arithmetic. 103 uint128_t vr, vp, vm; 104 int32_t e10; 105 bool vmIsTrailingZeros = false; 106 bool vrIsTrailingZeros = false; 107 if (e2 >= 0) { 108 // I tried special-casing q == 0, but there was no effect on performance. 109 // This expression is slightly faster than max(0, log10Pow2(e2) - 1). 110 const uint32_t q = log10Pow2(e2) - (e2 > 3); 111 e10 = q; 112 const int32_t k = FLOAT_128_POW5_INV_BITCOUNT + pow5bits(q) - 1; 113 const int32_t i = -e2 + q + k; 114 uint64_t pow5[4]; 115 generic_computeInvPow5(q, pow5); 116 vr = mulShift(4 * m2, pow5, i); 117 vp = mulShift(4 * m2 + 2, pow5, i); 118 vm = mulShift(4 * m2 - 1 - mmShift, pow5, i); 119#ifdef RYU_DEBUG 120 printf("%s * 2^%d / 10^%d\n", s(mv), e2, q); 121 printf("V+=%s\nV =%s\nV-=%s\n", s(vp), s(vr), s(vm)); 122#endif 123 // floor(log_5(2^128)) = 55, this is very conservative 124 if (q <= 55) { 125 // Only one of mp, mv, and mm can be a multiple of 5, if any. 126 if (mv % 5 == 0) { 127 vrIsTrailingZeros = multipleOfPowerOf5(mv, q - 1); 128 } else if (acceptBounds) { 129 // Same as min(e2 + (~mm & 1), pow5Factor(mm)) >= q 130 // <=> e2 + (~mm & 1) >= q && pow5Factor(mm) >= q 131 // <=> true && pow5Factor(mm) >= q, since e2 >= q. 132 vmIsTrailingZeros = multipleOfPowerOf5(mv - 1 - mmShift, q); 133 } else { 134 // Same as min(e2 + 1, pow5Factor(mp)) >= q. 135 vp -= multipleOfPowerOf5(mv + 2, q); 136 } 137 } 138 } else { 139 // This expression is slightly faster than max(0, log10Pow5(-e2) - 1). 140 const uint32_t q = log10Pow5(-e2) - (-e2 > 1); 141 e10 = q + e2; 142 const int32_t i = -e2 - q; 143 const int32_t k = pow5bits(i) - FLOAT_128_POW5_BITCOUNT; 144 const int32_t j = q - k; 145 uint64_t pow5[4]; 146 generic_computePow5(i, pow5); 147 vr = mulShift(4 * m2, pow5, j); 148 vp = mulShift(4 * m2 + 2, pow5, j); 149 vm = mulShift(4 * m2 - 1 - mmShift, pow5, j); 150#ifdef RYU_DEBUG 151 printf("%s * 5^%d / 10^%d\n", s(mv), -e2, q); 152 printf("%d %d %d %d\n", q, i, k, j); 153 printf("V+=%s\nV =%s\nV-=%s\n", s(vp), s(vr), s(vm)); 154#endif 155 if (q <= 1) { 156 // {vr,vp,vm} is trailing zeros if {mv,mp,mm} has at least q trailing 0 bits. 157 // mv = 4 m2, so it always has at least two trailing 0 bits. 158 vrIsTrailingZeros = true; 159 if (acceptBounds) { 160 // mm = mv - 1 - mmShift, so it has 1 trailing 0 bit iff mmShift == 1. 161 vmIsTrailingZeros = mmShift == 1; 162 } else { 163 // mp = mv + 2, so it always has at least one trailing 0 bit. 164 --vp; 165 } 166 } else if (q < 127) { // TODO(ulfjack): Use a tighter bound here. 167 // We need to compute min(ntz(mv), pow5Factor(mv) - e2) >= q-1 168 // <=> ntz(mv) >= q-1 && pow5Factor(mv) - e2 >= q-1 169 // <=> ntz(mv) >= q-1 (e2 is negative and -e2 >= q) 170 // <=> (mv & ((1 << (q-1)) - 1)) == 0 171 // We also need to make sure that the left shift does not overflow. 172 vrIsTrailingZeros = multipleOfPowerOf2(mv, q - 1); 173#ifdef RYU_DEBUG 174 printf("vr is trailing zeros=%s\n", vrIsTrailingZeros ? "true" : "false"); 175#endif 176 } 177 } 178#ifdef RYU_DEBUG 179 printf("e10=%d\n", e10); 180 printf("V+=%s\nV =%s\nV-=%s\n", s(vp), s(vr), s(vm)); 181 printf("vm is trailing zeros=%s\n", vmIsTrailingZeros ? "true" : "false"); 182 printf("vr is trailing zeros=%s\n", vrIsTrailingZeros ? "true" : "false"); 183#endif 184 185 // Step 4: Find the shortest decimal representation in the interval of legal representations. 186 uint32_t removed = 0; 187 uint8_t lastRemovedDigit = 0; 188 uint128_t output; 189 190 while (vp / 10 > vm / 10) { 191 vmIsTrailingZeros &= vm % 10 == 0; 192 vrIsTrailingZeros &= lastRemovedDigit == 0; 193 lastRemovedDigit = (uint8_t) (vr % 10); 194 vr /= 10; 195 vp /= 10; 196 vm /= 10; 197 ++removed; 198 } 199#ifdef RYU_DEBUG 200 printf("V+=%s\nV =%s\nV-=%s\n", s(vp), s(vr), s(vm)); 201 printf("d-10=%s\n", vmIsTrailingZeros ? "true" : "false"); 202#endif 203 if (vmIsTrailingZeros) { 204 while (vm % 10 == 0) { 205 vrIsTrailingZeros &= lastRemovedDigit == 0; 206 lastRemovedDigit = (uint8_t) (vr % 10); 207 vr /= 10; 208 vp /= 10; 209 vm /= 10; 210 ++removed; 211 } 212 } 213#ifdef RYU_DEBUG 214 printf("%s %d\n", s(vr), lastRemovedDigit); 215 printf("vr is trailing zeros=%s\n", vrIsTrailingZeros ? "true" : "false"); 216#endif 217 if (vrIsTrailingZeros && (lastRemovedDigit == 5) && (vr % 2 == 0)) { 218 // Round even if the exact numbers is .....50..0. 219 lastRemovedDigit = 4; 220 } 221 // We need to take vr+1 if vr is outside bounds or we need to round up. 222 output = vr + 223 ((vr == vm && (!acceptBounds || !vmIsTrailingZeros)) || (lastRemovedDigit >= 5)); 224 const int32_t exp = e10 + removed; 225 226#ifdef RYU_DEBUG 227 printf("V+=%s\nV =%s\nV-=%s\n", s(vp), s(vr), s(vm)); 228 printf("O=%s\n", s(output)); 229 printf("EXP=%d\n", exp); 230#endif 231 232 struct floating_decimal_128 fd; 233 fd.mantissa = output; 234 fd.exponent = exp; 235 fd.sign = ieeeSign; 236 return fd; 237} 238 239static inline int copy_special_str(char * const result, const struct floating_decimal_128 fd) { 240 if (fd.mantissa) { 241 memcpy(result, "NaN", 3); 242 return 3; 243 } 244 if (fd.sign) { 245 result[0] = '-'; 246 } 247 memcpy(result + fd.sign, "Infinity", 8); 248 return fd.sign + 8; 249} 250 251int generic_to_chars(const struct floating_decimal_128 v, char* const result) { 252 if (v.exponent == FD128_EXCEPTIONAL_EXPONENT) { 253 return copy_special_str(result, v); 254 } 255 256 // Step 5: Print the decimal representation. 257 int index = 0; 258 if (v.sign) { 259 result[index++] = '-'; 260 } 261 262 uint128_t output = v.mantissa; 263 const uint32_t olength = decimalLength(output); 264 265#ifdef RYU_DEBUG 266 printf("DIGITS=%s\n", s(v.mantissa)); 267 printf("OLEN=%u\n", olength); 268 printf("EXP=%u\n", v.exponent + olength); 269#endif 270 271 for (uint32_t i = 0; i < olength - 1; ++i) { 272 const uint32_t c = (uint32_t) (output % 10); 273 output /= 10; 274 result[index + olength - i] = (char) ('0' + c); 275 } 276 result[index] = '0' + (uint32_t) (output % 10); // output should be < 10 by now. 277 278 // Print decimal point if needed. 279 if (olength > 1) { 280 result[index + 1] = '.'; 281 index += olength + 1; 282 } else { 283 ++index; 284 } 285 286 // Print the exponent. 287 result[index++] = 'e'; 288 int32_t exp = v.exponent + olength - 1; 289 if (exp < 0) { 290 result[index++] = '-'; 291 exp = -exp; 292 } else 293 result[index++] = '+'; 294 295 uint32_t elength = decimalLength(exp); 296 if (elength == 1) 297 ++elength; 298 for (uint32_t i = 0; i < elength; ++i) { 299 const uint32_t c = exp % 10; 300 exp /= 10; 301 result[index + elength - 1 - i] = (char) ('0' + c); 302 } 303 index += elength; 304 return index; 305} 306