grisu3_math.h (11693B)
1 /* 2 * Copyright (c) 2016 Mikkel F. Jørgensen, dvide.com 3 * Copyright author of MathGeoLib (https://github.com/juj) 4 * 5 * Licensed under the Apache License, Version 2.0 (the "License"); 6 * you may not use this file except in compliance with the License. 7 * You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. http://www.apache.org/licenses/LICENSE-2.0 16 */ 17 18 /* 2016-02-02: Updated by mikkelfj 19 * 20 * Extracted from MatGeoLib grisu3.c, Apache 2.0 license, and extended. 21 * 22 * This file is usually include via grisu3_print.h or grisu3_parse.h. 23 * 24 * The original MatGeoLib dtoa_grisu3 implementation is largely 25 * unchanged except for the uint64 to double cast. The remaining changes 26 * are file structure, name changes, and new additions for parsing: 27 * 28 * - Split into header files only: 29 * grisu3_math.h, grisu3_print.h, (added grisu3_parse.h) 30 * 31 * - names prefixed with grisu3_, grisu3_diy_fp_, GRISU3_. 32 * - added static to all functions. 33 * - disabled clang unused function warnings. 34 * - guarded <stdint.h> to allow for alternative impl. 35 * - added extra numeric constants needed for parsing. 36 * - added dec_pow, cast_double_from_diy_fp. 37 * - changed some function names for consistency. 38 * - moved printing specific grisu3 functions to grisu3_print.h. 39 * - changed double to uint64 cast to avoid aliasing. 40 * - added new grisu3_parse.h for parsing doubles. 41 * - grisu3_print_double (dtoa_grisu3) format .1 as 0.1 needed for valid JSON output 42 * and grisu3_parse_double wouldn't consume it. 43 * - grsu3_print_double changed formatting to prefer 0.012 over 1.2e-2. 44 * 45 * These changes make it possible to include the files as headers only 46 * in other software libraries without risking name conflicts, and to 47 * extend the implementation with a port of Googles Double Conversion 48 * strtod functionality for parsing doubles. 49 * 50 * Extracted from: rev. 915501a / Dec 22, 2015 51 * <https://github.com/juj/MathGeoLib/blob/master/src/Math/grisu3.c> 52 * MathGeoLib License: http://www.apache.org/licenses/LICENSE-2.0.html 53 */ 54 55 #ifndef GRISU3_MATH_H 56 #define GRISU3_MATH_H 57 58 #ifdef __cplusplus 59 extern "C" { 60 #endif 61 62 /* Guarded to allow inclusion of pstdint.h first, if stdint.h is not supported. */ 63 #ifndef UINT8_MAX 64 #include <stdint.h> /* uint64_t etc. */ 65 #endif 66 67 #ifdef GRISU3_NO_ASSERT 68 #undef GRISU3_ASSERT 69 #define GRISU3_ASSERT(x) ((void)0) 70 #endif 71 72 #ifndef GRISU3_ASSERT 73 #include <assert.h> /* assert */ 74 #define GRISU3_ASSERT(x) assert(x) 75 #endif 76 77 #ifdef _MSC_VER 78 #pragma warning(disable : 4204) /* nonstandard extension used : non-constant aggregate initializer */ 79 #endif 80 81 #define GRISU3_D64_SIGN 0x8000000000000000ULL 82 #define GRISU3_D64_EXP_MASK 0x7FF0000000000000ULL 83 #define GRISU3_D64_FRACT_MASK 0x000FFFFFFFFFFFFFULL 84 #define GRISU3_D64_IMPLICIT_ONE 0x0010000000000000ULL 85 #define GRISU3_D64_EXP_POS 52 86 #define GRISU3_D64_EXP_BIAS 1075 87 #define GRISU3_D64_DENORM_EXP (-GRISU3_D64_EXP_BIAS + 1) 88 #define GRISU3_DIY_FP_FRACT_SIZE 64 89 #define GRISU3_D_1_LOG2_10 0.30102999566398114 /* 1 / lg(10) */ 90 #define GRISU3_MIN_TARGET_EXP -60 91 #define GRISU3_MASK32 0xFFFFFFFFULL 92 #define GRISU3_MIN_CACHED_EXP -348 93 #define GRISU3_MAX_CACHED_EXP 340 94 #define GRISU3_CACHED_EXP_STEP 8 95 #define GRISU3_D64_MAX_DEC_EXP 309 96 #define GRISU3_D64_MIN_DEC_EXP -324 97 #define GRISU3_D64_INF GRISU3_D64_EXP_MASK 98 99 #define GRISU3_MIN(x,y) ((x) <= (y) ? (x) : (y)) 100 #define GRISU3_MAX(x,y) ((x) >= (y) ? (x) : (y)) 101 102 103 typedef struct grisu3_diy_fp 104 { 105 uint64_t f; 106 int e; 107 } grisu3_diy_fp_t; 108 109 typedef struct grisu3_diy_fp_power 110 { 111 uint64_t fract; 112 int16_t b_exp, d_exp; 113 } grisu3_diy_fp_power_t; 114 115 typedef union { 116 uint64_t u64; 117 double d64; 118 } grisu3_cast_double_t; 119 120 static uint64_t grisu3_cast_uint64_from_double(double d) 121 { 122 grisu3_cast_double_t cd; 123 cd.d64 = d; 124 return cd.u64; 125 } 126 127 static double grisu3_cast_double_from_uint64(uint64_t u) 128 { 129 grisu3_cast_double_t cd; 130 cd.u64 = u; 131 return cd.d64; 132 } 133 134 #define grisu3_double_infinity grisu3_cast_double_from_uint64(GRISU3_D64_INF) 135 #define grisu3_double_nan grisu3_cast_double_from_uint64(GRISU3_D64_INF + 1) 136 137 static const grisu3_diy_fp_power_t grisu3_diy_fp_pow_cache[] = 138 { 139 { 0xfa8fd5a0081c0288ULL, -1220, -348 }, 140 { 0xbaaee17fa23ebf76ULL, -1193, -340 }, 141 { 0x8b16fb203055ac76ULL, -1166, -332 }, 142 { 0xcf42894a5dce35eaULL, -1140, -324 }, 143 { 0x9a6bb0aa55653b2dULL, -1113, -316 }, 144 { 0xe61acf033d1a45dfULL, -1087, -308 }, 145 { 0xab70fe17c79ac6caULL, -1060, -300 }, 146 { 0xff77b1fcbebcdc4fULL, -1034, -292 }, 147 { 0xbe5691ef416bd60cULL, -1007, -284 }, 148 { 0x8dd01fad907ffc3cULL, -980, -276 }, 149 { 0xd3515c2831559a83ULL, -954, -268 }, 150 { 0x9d71ac8fada6c9b5ULL, -927, -260 }, 151 { 0xea9c227723ee8bcbULL, -901, -252 }, 152 { 0xaecc49914078536dULL, -874, -244 }, 153 { 0x823c12795db6ce57ULL, -847, -236 }, 154 { 0xc21094364dfb5637ULL, -821, -228 }, 155 { 0x9096ea6f3848984fULL, -794, -220 }, 156 { 0xd77485cb25823ac7ULL, -768, -212 }, 157 { 0xa086cfcd97bf97f4ULL, -741, -204 }, 158 { 0xef340a98172aace5ULL, -715, -196 }, 159 { 0xb23867fb2a35b28eULL, -688, -188 }, 160 { 0x84c8d4dfd2c63f3bULL, -661, -180 }, 161 { 0xc5dd44271ad3cdbaULL, -635, -172 }, 162 { 0x936b9fcebb25c996ULL, -608, -164 }, 163 { 0xdbac6c247d62a584ULL, -582, -156 }, 164 { 0xa3ab66580d5fdaf6ULL, -555, -148 }, 165 { 0xf3e2f893dec3f126ULL, -529, -140 }, 166 { 0xb5b5ada8aaff80b8ULL, -502, -132 }, 167 { 0x87625f056c7c4a8bULL, -475, -124 }, 168 { 0xc9bcff6034c13053ULL, -449, -116 }, 169 { 0x964e858c91ba2655ULL, -422, -108 }, 170 { 0xdff9772470297ebdULL, -396, -100 }, 171 { 0xa6dfbd9fb8e5b88fULL, -369, -92 }, 172 { 0xf8a95fcf88747d94ULL, -343, -84 }, 173 { 0xb94470938fa89bcfULL, -316, -76 }, 174 { 0x8a08f0f8bf0f156bULL, -289, -68 }, 175 { 0xcdb02555653131b6ULL, -263, -60 }, 176 { 0x993fe2c6d07b7facULL, -236, -52 }, 177 { 0xe45c10c42a2b3b06ULL, -210, -44 }, 178 { 0xaa242499697392d3ULL, -183, -36 }, 179 { 0xfd87b5f28300ca0eULL, -157, -28 }, 180 { 0xbce5086492111aebULL, -130, -20 }, 181 { 0x8cbccc096f5088ccULL, -103, -12 }, 182 { 0xd1b71758e219652cULL, -77, -4 }, 183 { 0x9c40000000000000ULL, -50, 4 }, 184 { 0xe8d4a51000000000ULL, -24, 12 }, 185 { 0xad78ebc5ac620000ULL, 3, 20 }, 186 { 0x813f3978f8940984ULL, 30, 28 }, 187 { 0xc097ce7bc90715b3ULL, 56, 36 }, 188 { 0x8f7e32ce7bea5c70ULL, 83, 44 }, 189 { 0xd5d238a4abe98068ULL, 109, 52 }, 190 { 0x9f4f2726179a2245ULL, 136, 60 }, 191 { 0xed63a231d4c4fb27ULL, 162, 68 }, 192 { 0xb0de65388cc8ada8ULL, 189, 76 }, 193 { 0x83c7088e1aab65dbULL, 216, 84 }, 194 { 0xc45d1df942711d9aULL, 242, 92 }, 195 { 0x924d692ca61be758ULL, 269, 100 }, 196 { 0xda01ee641a708deaULL, 295, 108 }, 197 { 0xa26da3999aef774aULL, 322, 116 }, 198 { 0xf209787bb47d6b85ULL, 348, 124 }, 199 { 0xb454e4a179dd1877ULL, 375, 132 }, 200 { 0x865b86925b9bc5c2ULL, 402, 140 }, 201 { 0xc83553c5c8965d3dULL, 428, 148 }, 202 { 0x952ab45cfa97a0b3ULL, 455, 156 }, 203 { 0xde469fbd99a05fe3ULL, 481, 164 }, 204 { 0xa59bc234db398c25ULL, 508, 172 }, 205 { 0xf6c69a72a3989f5cULL, 534, 180 }, 206 { 0xb7dcbf5354e9beceULL, 561, 188 }, 207 { 0x88fcf317f22241e2ULL, 588, 196 }, 208 { 0xcc20ce9bd35c78a5ULL, 614, 204 }, 209 { 0x98165af37b2153dfULL, 641, 212 }, 210 { 0xe2a0b5dc971f303aULL, 667, 220 }, 211 { 0xa8d9d1535ce3b396ULL, 694, 228 }, 212 { 0xfb9b7cd9a4a7443cULL, 720, 236 }, 213 { 0xbb764c4ca7a44410ULL, 747, 244 }, 214 { 0x8bab8eefb6409c1aULL, 774, 252 }, 215 { 0xd01fef10a657842cULL, 800, 260 }, 216 { 0x9b10a4e5e9913129ULL, 827, 268 }, 217 { 0xe7109bfba19c0c9dULL, 853, 276 }, 218 { 0xac2820d9623bf429ULL, 880, 284 }, 219 { 0x80444b5e7aa7cf85ULL, 907, 292 }, 220 { 0xbf21e44003acdd2dULL, 933, 300 }, 221 { 0x8e679c2f5e44ff8fULL, 960, 308 }, 222 { 0xd433179d9c8cb841ULL, 986, 316 }, 223 { 0x9e19db92b4e31ba9ULL, 1013, 324 }, 224 { 0xeb96bf6ebadf77d9ULL, 1039, 332 }, 225 { 0xaf87023b9bf0ee6bULL, 1066, 340 } 226 }; 227 228 /* Avoid dependence on lib math to get (int)ceil(v) */ 229 static int grisu3_iceil(double v) 230 { 231 int k = (int)v; 232 if (v < 0) return k; 233 return v - k == 0 ? k : k + 1; 234 } 235 236 static int grisu3_diy_fp_cached_pow(int exp, grisu3_diy_fp_t *p) 237 { 238 int k = grisu3_iceil((exp+GRISU3_DIY_FP_FRACT_SIZE-1) * GRISU3_D_1_LOG2_10); 239 int i = (k-GRISU3_MIN_CACHED_EXP-1) / GRISU3_CACHED_EXP_STEP + 1; 240 p->f = grisu3_diy_fp_pow_cache[i].fract; 241 p->e = grisu3_diy_fp_pow_cache[i].b_exp; 242 return grisu3_diy_fp_pow_cache[i].d_exp; 243 } 244 245 static grisu3_diy_fp_t grisu3_diy_fp_minus(grisu3_diy_fp_t x, grisu3_diy_fp_t y) 246 { 247 grisu3_diy_fp_t d; d.f = x.f - y.f; d.e = x.e; 248 GRISU3_ASSERT(x.e == y.e && x.f >= y.f); 249 return d; 250 } 251 252 static grisu3_diy_fp_t grisu3_diy_fp_multiply(grisu3_diy_fp_t x, grisu3_diy_fp_t y) 253 { 254 uint64_t a, b, c, d, ac, bc, ad, bd, tmp; 255 grisu3_diy_fp_t r; 256 a = x.f >> 32; b = x.f & GRISU3_MASK32; 257 c = y.f >> 32; d = y.f & GRISU3_MASK32; 258 ac = a*c; bc = b*c; 259 ad = a*d; bd = b*d; 260 tmp = (bd >> 32) + (ad & GRISU3_MASK32) + (bc & GRISU3_MASK32); 261 tmp += 1U << 31; /* round */ 262 r.f = ac + (ad >> 32) + (bc >> 32) + (tmp >> 32); 263 r.e = x.e + y.e + 64; 264 return r; 265 } 266 267 static grisu3_diy_fp_t grisu3_diy_fp_normalize(grisu3_diy_fp_t n) 268 { 269 GRISU3_ASSERT(n.f != 0); 270 while(!(n.f & 0xFFC0000000000000ULL)) { n.f <<= 10; n.e -= 10; } 271 while(!(n.f & GRISU3_D64_SIGN)) { n.f <<= 1; --n.e; } 272 return n; 273 } 274 275 static grisu3_diy_fp_t grisu3_cast_diy_fp_from_double(double d) 276 { 277 grisu3_diy_fp_t fp; 278 uint64_t u64 = grisu3_cast_uint64_from_double(d); 279 if (!(u64 & GRISU3_D64_EXP_MASK)) { fp.f = u64 & GRISU3_D64_FRACT_MASK; fp.e = 1 - GRISU3_D64_EXP_BIAS; } 280 else { fp.f = (u64 & GRISU3_D64_FRACT_MASK) + GRISU3_D64_IMPLICIT_ONE; fp.e = (int)((u64 & GRISU3_D64_EXP_MASK) >> GRISU3_D64_EXP_POS) - GRISU3_D64_EXP_BIAS; } 281 return fp; 282 } 283 284 static double grisu3_cast_double_from_diy_fp(grisu3_diy_fp_t n) 285 { 286 const uint64_t hidden_bit = GRISU3_D64_IMPLICIT_ONE; 287 const uint64_t frac_mask = GRISU3_D64_FRACT_MASK; 288 const int denorm_exp = GRISU3_D64_DENORM_EXP; 289 const int exp_bias = GRISU3_D64_EXP_BIAS; 290 const int exp_pos = GRISU3_D64_EXP_POS; 291 292 grisu3_diy_fp_t v = n; 293 uint64_t e_biased; 294 295 while (v.f > hidden_bit + frac_mask) { 296 v.f >>= 1; 297 ++v.e; 298 } 299 if (v.e < denorm_exp) { 300 return 0.0; 301 } 302 while (v.e > denorm_exp && (v.f & hidden_bit) == 0) { 303 v.f <<= 1; 304 --v.e; 305 } 306 if (v.e == denorm_exp && (v.f & hidden_bit) == 0) { 307 e_biased = 0; 308 } else { 309 e_biased = (uint64_t)(v.e + exp_bias); 310 } 311 return grisu3_cast_double_from_uint64((v.f & frac_mask) | (e_biased << exp_pos)); 312 } 313 314 /* pow10_cache[i] = 10^(i-1) */ 315 static const unsigned int grisu3_pow10_cache[] = { 0, 1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000 }; 316 317 static int grisu3_largest_pow10(uint32_t n, int n_bits, uint32_t *power) 318 { 319 int guess = ((n_bits + 1) * 1233 >> 12) + 1/*skip first entry*/; 320 if (n < grisu3_pow10_cache[guess]) --guess; /* We don't have any guarantees that 2^n_bits <= n. */ 321 *power = grisu3_pow10_cache[guess]; 322 return guess; 323 } 324 325 #ifdef __cplusplus 326 } 327 #endif 328 329 #endif /* GRISU3_MATH_H */