code / tomo

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