ctfefloat.d

   1 // stb_sprintf - v1.05 - public domain snprintf() implementation
   2 // originally by Jeff Roberts / RAD Game Tools, 2015/10/20
   3 // http://github.com/nothings/stb
   4 //
   5 // allowed types:  sc uidBboXx p AaGgEef n
   6 // lengths      :  h ll j z t I64 I32 I
   7 //
   8 // Contributors:
   9 //    Fabian "ryg" Giesen (reformatting)
  10 //
  11 // Contributors (bugfixes):
  12 //    github:d26435
  13 //    github:trex78
  14 //    Jari Komppa (SI suffixes)
  15 //    Rohit Nirmal
  16 //    Marcin Wojdyr
  17 //    Leonard Ritter
  18 //
  19 // LICENSE:
  20 // ------------------------------------------------------------------------------
  21 // This software is available under 2 licenses -- choose whichever you prefer.
  22 // ------------------------------------------------------------------------------
  23 // ALTERNATIVE A - MIT License
  24 // Copyright (c) 2017 Sean Barrett
  25 // Permission is hereby granted, free of charge, to any person obtaining a copy of
  26 // this software and associated documentation files (the "Software"), to deal in
  27 // the Software without restriction, including without limitation the rights to
  28 // use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
  29 // of the Software, and to permit persons to whom the Software is furnished to do
  30 // so, subject to the following conditions:
  31 // The above copyright notice and this permission notice shall be included in all
  32 // copies or substantial portions of the Software.
  33 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  34 // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  35 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  36 // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  37 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  38 // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  39 // SOFTWARE.
  40 // ------------------------------------------------------------------------------
  41 // This program is free software: you can redistribute it and/or modify
  42 // it under the terms of the GNU General Public License as published by
  43 // the Free Software Foundation, version 3 of the License ONLY.
  44 //
  45 // This program is distributed in the hope that it will be useful,
  46 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  47 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  48 // GNU General Public License for more details.
  49 //
  50 // You should have received a copy of the GNU General Public License
  51 // along with this program. If not, see <http://www.gnu.org/licenses/>.
  52 module iv.ctfefloat;
  53
  54
  55 // ////////////////////////////////////////////////////////////////////////// //
  56 //pragma(msg, double2str(66.66));
  57
  58
  59 // ////////////////////////////////////////////////////////////////////////// //
  60 public string double2str (double d, int fracdigs=-1) nothrow @trusted {
  61   char[STBSP__NUMSZ] num;
  62   uint fdc = (fracdigs >= 0 ? cast(uint)fracdigs : 6);
  63   uint olen;
  64   int dpos;
  65   const(char)* start;
  66   int sign = stbsp__real_to_str(&start, &olen, num.ptr, &dpos, d, fdc);
  67   if (dpos == STBSP__SPECIAL) {
  68     return start[0..olen].idup;
  69   } else {
  70     while (olen > dpos && start[olen-1] == '0') --olen;
  71     uint reslen = (sign ? 1 : 0)+olen+(dpos == 0 ? 2 : 1);
  72     char[] res = new char[](reslen+(olen-dpos < fdc ? fdc-(olen-dpos) : 0));
  73     res[] = '0'; // why not?
  74     uint rpos = 0;
  75     if (sign) res.ptr[rpos++] = '-';
  76     if (dpos > 0) { res.ptr[rpos..rpos+dpos] = start[0..dpos]; rpos += dpos; }
  77     if (fdc != 0) {
  78       res.ptr[rpos++] = '.';
  79       res.ptr[rpos..rpos+(olen-dpos)] = start[dpos..olen];
  80       rpos += olen-dpos;
  81       if (fracdigs > 0 && olen-dpos < fdc) rpos += fdc-(olen-dpos);
  82     }
  83     assert(rpos <= res.length);
  84     return cast(string)res.ptr[0..rpos]; // it is safe to cast here
  85   }
  86 }
  87
  88
  89 // ////////////////////////////////////////////////////////////////////////// //
  90 private:
  91
  92 public enum STBSP__SPECIAL = 0x7000;
  93
  94 public enum STBSP__NUMSZ = 512; // big enough for e308 (with commas) or e-307
  95 //char[STBSP__NUMSZ] num;
  96
  97
  98 // ////////////////////////////////////////////////////////////////////////// //
  99 static immutable double[23] stbsp__bot = [
 100   0x1p+0,0x1.4p+3,0x1.9p+6,0x1.f4p+9,0x1.388p+13,0x1.86ap+16,0x1.e848p+19,0x1.312dp+23,
 101   0x1.7d784p+26,0x1.dcd65p+29,0x1.2a05f2p+33,0x1.74876e8p+36,0x1.d1a94a2p+39,0x1.2309ce54p+43,
 102   0x1.6bcc41e9p+46,0x1.c6bf52634p+49,0x1.1c37937e08p+53,0x1.6345785d8ap+56,0x1.bc16d674ec8p+59,
 103   0x1.158e460913dp+63,0x1.5af1d78b58c4p+66,0x1.b1ae4d6e2ef5p+69,0x1.0f0cf064dd592p+73,
 104 ];
 105
 106 static immutable double[22] stbsp__negbot = [
 107   0x1.999999999999ap-4,0x1.47ae147ae147bp-7,0x1.0624dd2f1a9fcp-10,0x1.a36e2eb1c432dp-14,
 108   0x1.4f8b588e368f1p-17,0x1.0c6f7a0b5ed8dp-20,0x1.ad7f29abcaf48p-24,0x1.5798ee2308c3ap-27,
 109   0x1.12e0be826d695p-30,0x1.b7cdfd9d7bdbbp-34,0x1.5fd7fe1796495p-37,0x1.19799812dea11p-40,
 110   0x1.c25c268497682p-44,0x1.6849b86a12b9bp-47,0x1.203af9ee75616p-50,0x1.cd2b297d889bcp-54,
 111   0x1.70ef54646d497p-57,0x1.2725dd1d243acp-60,0x1.d83c94fb6d2acp-64,0x1.79ca10c924223p-67,
 112   0x1.2e3b40a0e9b4fp-70,0x1.e392010175ee6p-74,
 113 ];
 114
 115 static immutable double[22] stbsp__negboterr = [
 116   -0x1.999999999999ap-58,-0x1.eb851eb851eb8p-63,-0x1.89374bc6a7efap-66,-0x1.6a161e4f765fep-68,
 117   -0x1.ee78183f91e64p-71,0x1.b5a63f9a49c2cp-75,0x1.5e1e99483b023p-78,-0x1.03023df2d4c94p-82,
 118   -0x1.34674bfabb83bp-84,-0x1.20a5465df8d2cp-88,0x1.7f7bc7b4d28aap-91,0x1.97f27f0f6e886p-96,
 119   -0x1.ecd79a5a0df95p-99,0x1.ea70909833de7p-107,-0x1.937831647f5ap-104,0x1.5b4c2ebe68799p-109,
 120   -0x1.db7b2080a3029p-111,-0x1.7c628066e8ceep-114,0x1.a52b31e9e3d07p-119,0x1.75447a5d8e536p-121,
 121   0x1.f769fb7e0b75ep-124,-0x1.a7566d9cba769p-128,
 122 ];
 123
 124 static immutable double[13] stbsp__top = [
 125   0x1.52d02c7e14af6p+76,0x1.c06a5ec5433c6p+152,0x1.28bc8abe49f64p+229,0x1.88ba3bf284e24p+305,
 126   0x1.03e29f5c2b18cp+382,0x1.57f48bb41db7cp+458,0x1.c73892ecbfbf4p+534,0x1.2d3d6f88f0b3dp+611,
 127   0x1.8eb0138858d0ap+687,0x1.07d457124123dp+764,0x1.5d2ce55747a18p+840,0x1.ce2137f743382p+916,
 128   0x1.31cfd3999f7bp+993,
 129 ];
 130
 131 static immutable double[13] stbsp__negtop = [
 132   0x1.82db34012b251p-77,0x1.244ce242c5561p-153,0x1.b9b6364f30304p-230,0x1.4dbf7b3f71cb7p-306,
 133   0x1.f8587e7083e3p-383,0x1.7d12a4670c123p-459,0x1.1fee341fc585dp-535,0x1.b31bb5dc320d2p-612,
 134   0x1.48c22ca71a1bdp-688,0x1.f0ce4839198dbp-765,0x1.77603725064a8p-841,0x1.1ba03f5b21p-917,
 135   0x1.ac9a7b3b7302fp-994,
 136 ];
 137
 138 static immutable double[13] stbsp__toperr = [
 139   0x1p+23,0x1.bb542c80deb4p+95,-0x1.83b80b9aab60ap+175,-0x1.32e22d17a166cp+251,
 140   -0x1.23606902e180ep+326,-0x1.96fb782462e87p+403,-0x1.358952c0bd011p+480,-0x1.78c1376a34b6cp+555,
 141   -0x1.17569fc243adfp+633,-0x1.d9365a897aaa6p+710,0x1.9050c256123ap+786,-0x1.b1799d76cc7a6p+862,
 142   -0x1.213fe39571a38p+939,
 143 ];
 144
 145 static immutable double[13] stbsp__negtoperr = [
 146   0x1.13badb829e079p-131,-0x1.e46a98d3d9f64p-209,0x1.227c7218a2b65p-284,0x1.1d96999aa01e9p-362,
 147   -0x1.cc2229efc3962p-437,-0x1.cd04a2263407ap-513,-0x1.23b80f187a157p-590,-0x1.c4e22914ed912p-666,
 148   0x1.bc296cdf42f82p-742,-0x1.9f9e7f4e16fe1p-819,-0x1.aeb0a72a8902ap-895,-0x1.e228e12c13408p-971,
 149   0x0.0000000fa1259p-1022,
 150 ];
 151
 152
 153 static immutable ulong[20] stbsp__powten = [
 154    1UL,
 155    10UL,
 156    100UL,
 157    1000UL,
 158    10000UL,
 159    100000UL,
 160    1000000UL,
 161    10000000UL,
 162    100000000UL,
 163    1000000000UL,
 164    10000000000UL,
 165    100000000000UL,
 166    1000000000000UL,
 167    10000000000000UL,
 168    100000000000000UL,
 169    1000000000000000UL,
 170    10000000000000000UL,
 171    100000000000000000UL,
 172    1000000000000000000UL,
 173    10000000000000000000UL
 174 ];
 175 enum stbsp__tento19th = 1000000000000000000UL;
 176
 177 static immutable string stbsp__digitpair = "00010203040506070809101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899";
 178
 179
 180 // get float info; returns `-1` for negatives (including negative zero), and `0` for positives
 181 public int stbsp__real_to_parts (long* bits, int* expo, double value) nothrow @trusted @nogc {
 182   // load value and round at the frac_digits
 183   double d = value;
 184   long b = *(cast(const(long)*)&d);
 185   //STBSP__COPYFP(b, d);
 186
 187   *bits = b&(((cast(ulong)1)<<52)-1);
 188   *expo = cast(int)(((b>>52)&2047)-1023);
 189
 190   return cast(int)(b>>63);
 191 }
 192
 193
 194 void stbsp__ddmulthi (ref double oh, ref double ol, in double xh, in double yh) nothrow @trusted @nogc {
 195   double ahi = 0, alo, bhi = 0, blo;
 196   long bt;
 197   oh = xh * yh;
 198   //STBSP__COPYFP(bt, xh);
 199   bt = *(cast(const(long)*)&xh);
 200   bt &= ((~cast(ulong)0) << 27);
 201   //STBSP__COPYFP(ahi, bt);
 202   ahi = *(cast(const(double)*)&bt);
 203   alo = xh - ahi;
 204   //STBSP__COPYFP(bt, yh);
 205   bt = *(cast(const(long)*)&yh);
 206   bt &= ((~cast(ulong)0) << 27);
 207   //STBSP__COPYFP(bhi, bt);
 208   bhi = *(cast(const(double)*)&bt);
 209   blo = yh - bhi;
 210   ol = ((ahi * bhi - oh) + ahi * blo + alo * bhi) + alo * blo;
 211 }
 212
 213 void stbsp__ddtoS64 (ref long ob, in double xh, in double xl, in double ph) nothrow @trusted @nogc {
 214   double ahi = 0, alo, vh, t;
 215   ob = cast(long)ph;
 216   vh = cast(double)ob;
 217   ahi = (xh - vh);
 218   t = (ahi - xh);
 219   alo = (xh - (ahi - t)) - (vh + t);
 220   ob += cast(long)(ahi + alo + xl);
 221 }
 222
 223 void stbsp__ddrenorm (ref double oh, ref double ol) nothrow @trusted @nogc {
 224   double s = oh + ol;
 225   ol = ol - (s - oh);
 226   oh = s;
 227 }
 228
 229 void stbsp__ddmultlo (in double oh, ref double ol, in double xh, in double xl, in double yh, in double yl) nothrow @trusted @nogc { ol = ol + (xh * yl + xl * yh); }
 230
 231 void stbsp__ddmultlos (in double oh, ref double ol, in double xh, in double yl) nothrow @trusted @nogc { ol = ol + (xh * yl); }
 232
 233 // power can be -323 to +350
 234 void stbsp__raise_to_power10 (double* ohi, double* olo, double d, int power) nothrow @trusted @nogc {
 235   double ph, pl;
 236   if (power >= 0 && power <= 22) {
 237     stbsp__ddmulthi(ph, pl, d, stbsp__bot[power]);
 238   } else {
 239     int e, et, eb;
 240     double p2h, p2l;
 241
 242     e = power;
 243     if (power < 0) e = -e;
 244     et = (e * 0x2c9) >> 14; /* %23 */
 245     if (et > 13) et = 13;
 246     eb = e - (et * 23);
 247
 248     ph = d;
 249     pl = 0.0;
 250     if (power < 0) {
 251       if (eb) {
 252         --eb;
 253         stbsp__ddmulthi(ph, pl, d, stbsp__negbot[eb]);
 254         stbsp__ddmultlos(ph, pl, d, stbsp__negboterr[eb]);
 255       }
 256       if (et) {
 257         stbsp__ddrenorm(ph, pl);
 258         --et;
 259         stbsp__ddmulthi(p2h, p2l, ph, stbsp__negtop[et]);
 260         stbsp__ddmultlo(p2h, p2l, ph, pl, stbsp__negtop[et], stbsp__negtoperr[et]);
 261         ph = p2h;
 262         pl = p2l;
 263       }
 264     } else {
 265       if (eb) {
 266         e = eb;
 267         if (eb > 22) eb = 22;
 268         e -= eb;
 269         stbsp__ddmulthi(ph, pl, d, stbsp__bot[eb]);
 270         if (e) {
 271           stbsp__ddrenorm(ph, pl);
 272           stbsp__ddmulthi(p2h, p2l, ph, stbsp__bot[e]);
 273           stbsp__ddmultlos(p2h, p2l, stbsp__bot[e], pl);
 274           ph = p2h;
 275           pl = p2l;
 276         }
 277       }
 278       if (et) {
 279         stbsp__ddrenorm(ph, pl);
 280         --et;
 281         stbsp__ddmulthi(p2h, p2l, ph, stbsp__top[et]);
 282         stbsp__ddmultlo(p2h, p2l, ph, pl, stbsp__top[et], stbsp__toperr[et]);
 283         ph = p2h;
 284         pl = p2l;
 285       }
 286     }
 287   }
 288   stbsp__ddrenorm(ph, pl);
 289   *ohi = ph;
 290   *olo = pl;
 291 }
 292
 293
 294 // given a float value, returns the significant bits in bits, and the position of the decimal point in decimal_pos.
 295 // +INF/-INF and NAN are specified by special values returned in the decimal_pos parameter.
 296 // frac_digits is absolute normally, but if you want from first significant digits (got %g and %e), or in 0x80000000.
 297 // returns `-1` for negatives (including negative zero), and `0` for positives.
 298 public int stbsp__real_to_str (const(char)** start, uint* len, char* outbuf, int* decimal_pos, double value, uint frac_digits=6) nothrow @trusted @nogc {
 299   double d;
 300   long bits = 0;
 301   int expo, e, ng, tens;
 302
 303   d = value;
 304   //STBSP__COPYFP(bits, d);
 305   bits = *(cast(const(long)*)&d);
 306   expo = cast(int)((bits >> 52) & 2047);
 307   ng = cast(int)(bits >> 63);
 308   if (ng) d = -d;
 309
 310   // is nan or inf?
 311   if (expo == 2047) {
 312     *start = (bits & (((cast(ulong)1) << 52) - 1)) ? "NaN" : "Inf";
 313     *decimal_pos = STBSP__SPECIAL;
 314     *len = 3;
 315     return ng;
 316   }
 317
 318   // is zero or denormal?
 319   if (expo == 0) {
 320     if ((bits << 1) == 0) {
 321       // do zero
 322       *decimal_pos = 1;
 323       *start = outbuf;
 324       outbuf[0] = '0';
 325       *len = 1;
 326       return ng;
 327     }
 328     // find the right expo for denormals
 329     {
 330       long v = (cast(ulong)1) << 51;
 331       while ((bits & v) == 0) {
 332         --expo;
 333         v >>= 1;
 334       }
 335     }
 336   }
 337
 338   // find the decimal exponent as well as the decimal bits of the value
 339   {
 340     double ph, pl;
 341
 342     // log10 estimate - very specifically tweaked to hit or undershoot by no more than 1 of log10 of all expos 1..2046
 343     tens = expo - 1023;
 344     tens = (tens < 0) ? ((tens * 617) / 2048) : (((tens * 1233) / 4096) + 1);
 345
 346     // move the significant bits into position and stick them into an int
 347     stbsp__raise_to_power10(&ph, &pl, d, 18 - tens);
 348
 349     // get full as much precision from double-double as possible
 350     stbsp__ddtoS64(bits, ph, pl, ph);
 351
 352     // check if we undershot
 353     if ((cast(ulong)bits) >= stbsp__tento19th) ++tens;
 354   }
 355
 356   // now do the rounding in integer land
 357   frac_digits = (frac_digits & 0x80000000) ? ((frac_digits & 0x7ffffff) + 1) : (tens + frac_digits);
 358   if ((frac_digits < 24)) {
 359     uint dg = 1;
 360     if (cast(ulong)bits >= stbsp__powten[9]) dg = 10;
 361     while (cast(ulong)bits >= stbsp__powten[dg]) {
 362       ++dg;
 363       if (dg == 20) goto noround;
 364     }
 365     if (frac_digits < dg) {
 366       ulong r;
 367       // add 0.5 at the right position and round
 368       e = dg - frac_digits;
 369       if (cast(uint)e >= 24) goto noround;
 370       r = stbsp__powten[e];
 371       bits = bits + (r / 2);
 372       if (cast(ulong)bits >= stbsp__powten[dg])
 373          ++tens;
 374       bits /= r;
 375     }
 376   noround:;
 377   }
 378
 379   // kill long trailing runs of zeros
 380   if (bits) {
 381     uint n;
 382     for (;;) {
 383       if (bits <= 0xffffffff) break;
 384       if (bits % 1000) goto donez;
 385       bits /= 1000;
 386     }
 387     n = cast(uint)bits;
 388     while ((n % 1000) == 0) n /= 1000;
 389     bits = n;
 390   donez:;
 391   }
 392
 393   // convert to string
 394   outbuf += 64;
 395   e = 0;
 396   for (;;) {
 397     uint n;
 398     char* o = outbuf-8;
 399     // do the conversion in chunks of U32s (avoid most 64-bit divides, worth it, constant denomiators be damned)
 400     if (bits >= 100000000) {
 401       n = cast(uint)(bits % 100000000);
 402       bits /= 100000000;
 403     } else {
 404       n = cast(uint)bits;
 405       bits = 0;
 406     }
 407     while (n) {
 408       outbuf -= 2;
 409       if (__ctfe) {
 410         outbuf[0..2] = stbsp__digitpair[(n%100)*2..(n%100)*2+2];
 411       } else {
 412         *cast(ushort*)outbuf = *cast(ushort*)&stbsp__digitpair[(n%100)*2];
 413       }
 414       n /= 100;
 415       e += 2;
 416     }
 417     if (bits == 0) {
 418       if (e && outbuf[0] == '0') {
 419         ++outbuf;
 420         --e;
 421       }
 422       break;
 423     }
 424     while (outbuf != o) {
 425       *--outbuf = '0';
 426       ++e;
 427     }
 428   }
 429
 430   *decimal_pos = tens;
 431   *start = outbuf;
 432   *len = e;
 433   return ng;
 434 }