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