nostrdb

an unfairly fast embedded nostr database backed by lmdb
git clone git://jb55.com/nostrdb
Log | Files | Refs | Submodules | README | LICENSE

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 */