src/misc/fpconv.c

   1 #include <stdbool.h>
   2 #include <string.h>
   3
   4 #include "fpconv.h"
   5 #include "powers.h"
   6
   7 #define fracmask  0x000FFFFFFFFFFFFFU
   8 #define expmask   0x7FF0000000000000U
   9 #define hiddenbit 0x0010000000000000U
  10 #define signmask  0x8000000000000000U
  11 #define expbias   (1023 + 52)
  12
  13 #define absv(n) ((n) < 0 ? -(n) : (n))
  14 #define minv(a, b) ((a) < (b) ? (a) : (b))
  15
  16 static uint64_t tens[] = {
  17     10000000000000000000U, 1000000000000000000U, 100000000000000000U,
  18     10000000000000000U, 1000000000000000U, 100000000000000U,
  19     10000000000000U, 1000000000000U, 100000000000U,
  20     10000000000U, 1000000000U, 100000000U,
  21     10000000U, 1000000U, 100000U,
  22     10000U, 1000U, 100U,
  23     10U, 1U
  24 };
  25
  26 static inline uint64_t get_dbits(double d)
  27 {
  28     union {
  29         double   dbl;
  30         uint64_t i;
  31     } dbl_bits = { d };
  32
  33     return dbl_bits.i;
  34 }
  35
  36 static Fp build_fp(double d)
  37 {
  38     uint64_t bits = get_dbits(d);
  39
  40     Fp fp;
  41     fp.frac = bits & fracmask;
  42     fp.exp = (bits & expmask) >> 52;
  43
  44     if(fp.exp) {
  45         fp.frac += hiddenbit;
  46         fp.exp -= expbias;
  47
  48     } else {
  49         fp.exp = -expbias + 1;
  50     }
  51
  52     return fp;
  53 }
  54
  55 static void normalize(Fp* fp)
  56 {
  57     while ((fp->frac & hiddenbit) == 0) {
  58         fp->frac <<= 1;
  59         fp->exp--;
  60     }
  61
  62     int shift = 64 - 52 - 1;
  63     fp->frac <<= shift;
  64     fp->exp -= shift;
  65 }
  66
  67 static void get_normalized_boundaries(Fp* fp, Fp* lower, Fp* upper)
  68 {
  69     upper->frac = (fp->frac << 1) + 1;
  70     upper->exp  = fp->exp - 1;
  71
  72     while ((upper->frac & (hiddenbit << 1)) == 0) {
  73         upper->frac <<= 1;
  74         upper->exp--;
  75     }
  76
  77     int u_shift = 64 - 52 - 2;
  78
  79     upper->frac <<= u_shift;
  80     upper->exp = upper->exp - u_shift;
  81
  82
  83     int l_shift = fp->frac == hiddenbit ? 2 : 1;
  84
  85     lower->frac = (fp->frac << l_shift) - 1;
  86     lower->exp = fp->exp - l_shift;
  87
  88
  89     lower->frac <<= lower->exp - upper->exp;
  90     lower->exp = upper->exp;
  91 }
  92
  93 static Fp multiply(Fp* a, Fp* b)
  94 {
  95     const uint64_t lomask = 0x00000000FFFFFFFF;
  96
  97     uint64_t ah_bl = (a->frac >> 32)    * (b->frac & lomask);
  98     uint64_t al_bh = (a->frac & lomask) * (b->frac >> 32);
  99     uint64_t al_bl = (a->frac & lomask) * (b->frac & lomask);
 100     uint64_t ah_bh = (a->frac >> 32)    * (b->frac >> 32);
 101
 102     uint64_t tmp = (ah_bl & lomask) + (al_bh & lomask) + (al_bl >> 32);
 103     /* round up */
 104     tmp += 1U << 31;
 105
 106     Fp fp = {
 107         ah_bh + (ah_bl >> 32) + (al_bh >> 32) + (tmp >> 32),
 108         a->exp + b->exp + 64
 109     };
 110
 111     return fp;
 112 }
 113
 114 static void round_digit(char* digits, int ndigits, uint64_t delta, uint64_t rem, uint64_t kappa, uint64_t frac)
 115 {
 116     while (rem < frac && delta - rem >= kappa &&
 117            (rem + kappa < frac || frac - rem > rem + kappa - frac)) {
 118
 119         digits[ndigits - 1]--;
 120         rem += kappa;
 121     }
 122 }
 123
 124 static int generate_digits(Fp* fp, Fp* upper, Fp* lower, char* digits, int* K)
 125 {
 126     uint64_t wfrac = upper->frac - fp->frac;
 127     uint64_t delta = upper->frac - lower->frac;
 128
 129     Fp one;
 130     one.frac = 1ULL << -upper->exp;
 131     one.exp  = upper->exp;
 132
 133     uint64_t part1 = upper->frac >> -one.exp;
 134     uint64_t part2 = upper->frac & (one.frac - 1);
 135
 136     int idx = 0, kappa = 10;
 137     uint64_t* divp;
 138     /* 1000000000 */
 139     for(divp = tens + 10; kappa > 0; divp++) {
 140
 141         uint64_t div = *divp;
 142         unsigned digit = part1 / div;
 143
 144         if (digit || idx) {
 145             digits[idx++] = digit + '0';
 146         }
 147
 148         part1 -= digit * div;
 149         kappa--;
 150
 151         uint64_t tmp = (part1 <<-one.exp) + part2;
 152         if (tmp <= delta) {
 153             *K += kappa;
 154             round_digit(digits, idx, delta, tmp, div << -one.exp, wfrac);
 155
 156             return idx;
 157         }
 158     }
 159
 160     /* 10 */
 161     uint64_t* unit = tens + 18;
 162
 163     while(true) {
 164         part2 *= 10;
 165         delta *= 10;
 166         kappa--;
 167
 168         unsigned digit = part2 >> -one.exp;
 169         if (digit || idx) {
 170             digits[idx++] = digit + '0';
 171         }
 172
 173         part2 &= one.frac - 1;
 174         if (part2 < delta) {
 175             *K += kappa;
 176             round_digit(digits, idx, delta, part2, one.frac, wfrac * *unit);
 177
 178             return idx;
 179         }
 180
 181         unit--;
 182     }
 183 }
 184
 185 static int grisu2(double d, char* digits, int* K)
 186 {
 187     Fp w = build_fp(d);
 188
 189     Fp lower, upper;
 190     get_normalized_boundaries(&w, &lower, &upper);
 191
 192     normalize(&w);
 193
 194     int k;
 195     Fp cp = find_cachedpow10(upper.exp, &k);
 196
 197     w     = multiply(&w,     &cp);
 198     upper = multiply(&upper, &cp);
 199     lower = multiply(&lower, &cp);
 200
 201     lower.frac++;
 202     upper.frac--;
 203
 204     *K = -k;
 205
 206     return generate_digits(&w, &upper, &lower, digits, K);
 207 }
 208
 209 static int emit_digits(char* digits, int ndigits, char* dest, int K, bool neg)
 210 {
 211     int exp = absv(K + ndigits - 1);
 212
 213     /* write plain integer */
 214     if(K >= 0 && (exp < (ndigits + 7))) {
 215         memcpy(dest, digits, ndigits);
 216         memset(dest + ndigits, '0', K);
 217
 218         return ndigits + K;
 219     }
 220
 221     /* write decimal w/o scientific notation */
 222     if(K < 0 && (K > -7 || exp < 4)) {
 223         int offset = ndigits - absv(K);
 224         /* fp < 1.0 -> write leading zero */
 225         if(offset <= 0) {
 226             offset = -offset;
 227             dest[0] = '0';
 228             dest[1] = '.';
 229             memset(dest + 2, '0', offset);
 230             memcpy(dest + offset + 2, digits, ndigits);
 231
 232             return ndigits + 2 + offset;
 233
 234         /* fp > 1.0 */
 235         } else {
 236             memcpy(dest, digits, offset);
 237             dest[offset] = '.';
 238             memcpy(dest + offset + 1, digits + offset, ndigits - offset);
 239
 240             return ndigits + 1;
 241         }
 242     }
 243
 244     /* write decimal w/ scientific notation */
 245     ndigits = minv(ndigits, 18 - neg);
 246
 247     int idx = 0;
 248     dest[idx++] = digits[0];
 249
 250     if(ndigits > 1) {
 251         dest[idx++] = '.';
 252         memcpy(dest + idx, digits + 1, ndigits - 1);
 253         idx += ndigits - 1;
 254     }
 255
 256     dest[idx++] = 'e';
 257
 258     char sign = K + ndigits - 1 < 0 ? '-' : '+';
 259     dest[idx++] = sign;
 260
 261     int cent = 0;
 262
 263     if(exp > 99) {
 264         cent = exp / 100;
 265         dest[idx++] = cent + '0';
 266         exp -= cent * 100;
 267     }
 268     if(exp > 9) {
 269         int dec = exp / 10;
 270         dest[idx++] = dec + '0';
 271         exp -= dec * 10;
 272
 273     } else if(cent) {
 274         dest[idx++] = '0';
 275     }
 276
 277     dest[idx++] = exp % 10 + '0';
 278
 279     return idx;
 280 }
 281
 282 static int filter_special(double fp, char* dest)
 283 {
 284     if(fp == 0.0) {
 285         dest[0] = '0';
 286         return 1;
 287     }
 288
 289     uint64_t bits = get_dbits(fp);
 290
 291     bool nan = (bits & expmask) == expmask;
 292
 293     if(!nan) {
 294         return 0;
 295     }
 296
 297     if(bits & fracmask) {
 298         dest[0] = 'n'; dest[1] = 'a'; dest[2] = 'n';
 299
 300     } else {
 301         dest[0] = 'i'; dest[1] = 'n'; dest[2] = 'f';
 302     }
 303
 304     return 3;
 305 }
 306
 307 int fpconv_dtoa(double d, char dest[24])
 308 {
 309     char digits[18];
 310
 311     int str_len = 0;
 312     bool neg = false;
 313
 314     if(get_dbits(d) & signmask) {
 315         dest[0] = '-';
 316         str_len++;
 317         neg = true;
 318     }
 319
 320     int spec = filter_special(d, dest + str_len);
 321
 322     if(spec) {
 323         return str_len + spec;
 324     }
 325
 326     int K = 0;
 327     int ndigits = grisu2(d, digits, &K);
 328
 329     str_len += emit_digits(digits, ndigits, dest + str_len, K, neg);
 330
 331     return str_len;
 332 }