1 /* derived from /netlib/fp/dtoa.c assuming IEEE, Standard C */
2 /* kudos to dmg@bell-labs.com, gripes to ehg@bell-labs.com */
4 /* Let x be the exact mathematical number defined by a decimal
5 * string s. Then atof(s) is the round-nearest-even IEEE
6 * floating point value.
7 * Let y be an IEEE floating point value and let s be the string
8 * printed as %.17g. Then atof(s) is exactly y.
13 static Lock _dtoalk[2];
14 #define ACQUIRE_DTOA_LOCK(n) lock(&_dtoalk[n])
15 #define FREE_DTOA_LOCK(n) unlock(&_dtoalk[n])
16 #define PRIVATE_mem ((2000+sizeof(double)-1)/sizeof(double))
17 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
20 #define DBL_MAX_10_EXP 308
21 #define DBL_MAX_EXP 1024
23 #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
24 #define fpword0(x) ((FPdbleword*)&x)->hi
25 #define fpword1(x) ((FPdbleword*)&x)->lo
26 /* Ten_pmax = floor(P*log(2)/log(5)) */
27 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
28 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
29 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
33 #define Exp_msk1 0x100000
34 #define Exp_msk11 0x100000
35 #define Exp_mask 0x7ff00000
39 #define Exp_1 0x3ff00000
40 #define Exp_11 0x3ff00000
42 #define Frac_mask 0xfffff
43 #define Frac_mask1 0xfffff
46 #define Bndry_mask 0xfffff
47 #define Bndry_mask1 0xfffff
49 #define Sign_bit 0x80000000
55 #define Avoid_Underflow
57 #define rounded_product(a,b) a *= b
58 #define rounded_quotient(a,b) a /= b
60 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
61 #define Big1 0xffffffff
63 #define FFFFFFFF 0xffffffffUL
72 int k, maxwds, sign, wds;
76 typedef struct Bigint Bigint;
78 static Bigint *freelist[Kmax+1];
88 if (rv = freelist[k]) {
89 freelist[k] = rv->next;
92 len = (sizeof(Bigint) + (x - 1) * sizeof(unsigned int) + sizeof(double) -1)
94 if (pmem_next - private_mem + len <= PRIVATE_mem) {
95 rv = (Bigint * )pmem_next;
98 rv = (Bigint * )malloc(len * sizeof(double));
103 rv->sign = rv->wds = 0;
111 ACQUIRE_DTOA_LOCK(0);
112 v->next = freelist[v->k];
118 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
119 y->wds*sizeof(int) + 2*sizeof(int))
122 multadd(Bigint *b, int m, int a) /* multiply by m and add a */
125 unsigned int carry, *x, y;
135 y = (xi & 0xffff) * m + carry;
136 z = (xi >> 16) * m + (y >> 16);
138 *x++ = (z << 16) + (y & 0xffff);
141 if (wds >= b->maxwds) {
142 b1 = Balloc(b->k + 1);
154 s2b(const char *s, int nd0, int nd, unsigned int y9)
161 for (k = 0, y = 1; x > y; y <<= 1, k++)
171 b = multadd(b, 10, *s++ - '0');
177 b = multadd(b, 10, *s++ - '0');
182 hi0bits(register unsigned int x)
186 if (!(x & 0xffff0000)) {
190 if (!(x & 0xff000000)) {
194 if (!(x & 0xf0000000)) {
198 if (!(x & 0xc0000000)) {
202 if (!(x & 0x80000000)) {
204 if (!(x & 0x40000000))
211 lo0bits(unsigned int *y)
214 register unsigned int x = *y;
265 mult(Bigint *a, Bigint *b)
269 unsigned int * x, *xa, *xae, *xb, *xbe, *xc, *xc0;
271 unsigned int carry, z;
274 if (a->wds < b->wds) {
286 for (x = c->x, xa = x + wc; x < xa; x++)
293 for (; xb < xbe; xb++, xc0++) {
294 if (y = *xb & 0xffff) {
299 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
301 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
313 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
316 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
322 for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc)
331 pow5mult(Bigint *b, int k)
333 Bigint * b1, *p5, *p51;
335 static int p05[3] = {
339 b = multadd(b, p05[i-1], 0);
345 ACQUIRE_DTOA_LOCK(1);
360 if (!(p51 = p5->next)) {
361 ACQUIRE_DTOA_LOCK(1);
362 if (!(p51 = p5->next)) {
363 p51 = p5->next = mult(p5, p5);
374 lshift(Bigint *b, int k)
378 unsigned int * x, *x1, *xe, z;
383 for (i = b->maxwds; n1 > i; i <<= 1)
387 for (i = 0; i < n; i++)
410 cmp(Bigint *a, Bigint *b)
412 unsigned int * xa, *xa0, *xb, *xb0;
425 return * xa < *xb ? -1 : 1;
433 diff(Bigint *a, Bigint *b)
437 unsigned int * xa, *xae, *xb, *xbe, *xc;
438 unsigned int borrow, y;
466 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
467 borrow = (y & 0x10000) >> 16;
468 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
469 borrow = (z & 0x10000) >> 16;
473 y = (*xa & 0xffff) - borrow;
474 borrow = (y & 0x10000) >> 16;
475 z = (*xa++ >> 16) - borrow;
476 borrow = (z & 0x10000) >> 16;
491 L = (fpword0(x) & Exp_mask) - (P - 1) * Exp_msk1;
498 b2d(Bigint *a, int *e)
500 unsigned int * xa, *xa0, w, y, z;
503 #define d0 fpword0(d)
504 #define d1 fpword1(d)
512 d0 = Exp_1 | y >> Ebits - k;
513 w = xa > xa0 ? *--xa : 0;
514 d1 = y << (32 - Ebits) + k | w >> Ebits - k;
517 z = xa > xa0 ? *--xa : 0;
519 d0 = Exp_1 | y << k | z >> 32 - k;
520 y = xa > xa0 ? *--xa : 0;
521 d1 = z << k | y >> 32 - k;
533 d2b(double d, int *e, int *bits)
537 unsigned int * x, y, z;
538 #define d0 fpword0(d)
539 #define d1 fpword1(d)
545 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
546 de = (int)(d0 >> Exp_shift);
549 if (k = lo0bits(&y)) {
550 x[0] = y | z << 32 - k;
554 i = b->wds = (x[1] = z) ? 2 : 1;
561 *e = de - Bias - (P - 1) + k;
570 ratio(Bigint *a, Bigint *b)
577 k = ka - kb + 32 * (a->wds - b->wds);
579 fpword0(da) += k * Exp_msk1;
582 fpword0(db) += k * Exp_msk1;
589 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
590 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
596 1e16, 1e32, 1e64, 1e128, 1e256 };
598 static const double tinytens[] = {
599 1e-16, 1e-32, 1e-64, 1e-128,
600 9007199254740992.e-256
603 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
604 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
605 #define Scale_Bit 0x10
608 #define NAN_WORD0 0x7ff80000
613 match(const char **sp, char *t)
616 const char * s = *sp;
619 if ((c = *++s) >= 'A' && c <= 'Z')
629 gethex(double *rvp, const char **sp)
631 unsigned int c, x[2];
633 int havedig, udx0, xshift;
636 havedig = xshift = 0;
639 while (c = *(const unsigned char * )++s) {
640 if (c >= '0' && c <= '9')
642 else if (c >= 'a' && c <= 'f')
644 else if (c >= 'A' && c <= 'F')
647 if (udx0 && havedig) {
652 } else if (/*(*/ c == ')') {
656 return; /* invalid form: don't change *sp */
664 x[0] = (x[0] << 4) | (x[1] >> 28);
665 x[1] = (x[1] << 4) | c;
667 if ((x[0] &= 0xfffff) || x[1]) {
668 fpword0(*rvp) = Exp_mask | x[0];
669 fpword1(*rvp) = x[1];
674 quorem(Bigint *b, Bigint *S)
677 unsigned int * bx, *bxe, q, *sx, *sxe;
678 unsigned int borrow, carry, y, ys;
679 unsigned int si, z, zs;
688 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
694 ys = (si & 0xffff) * q + carry;
695 zs = (si >> 16) * q + (ys >> 16);
697 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
698 borrow = (y & 0x10000) >> 16;
699 z = (*bx >> 16) - (zs & 0xffff) - borrow;
700 borrow = (z & 0x10000) >> 16;
705 while (--bxe > bx && !*bxe)
710 if (cmp(b, S) >= 0) {
718 ys = (si & 0xffff) + carry;
719 zs = (si >> 16) + (ys >> 16);
721 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
722 borrow = (y & 0x10000) >> 16;
723 z = (*bx >> 16) - (zs & 0xffff) - borrow;
724 borrow = (z & 0x10000) >> 16;
730 while (--bxe > bx && !*bxe)
743 j = sizeof(unsigned int);
745 sizeof(Bigint) - sizeof(unsigned int) - sizeof(int) + j <= i;
748 r = (int * )Balloc(k);
755 nrv_alloc(char *s, char **rve, int n)
759 t = rv = rv_alloc(n);
767 /* freedtoa(s) must be used to free values s returned by dtoa
768 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
769 * but for consistency with earlier versions of dtoa, it is optional
770 * when MULTIPLE_THREADS is not defined.
776 Bigint * b = (Bigint * )((int *)s - 1);
777 b->maxwds = 1 << (b->k = *(int * )b);
781 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
783 * Inspired by "How to Print Floating-Point Numbers Accurately" by
784 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
787 * 1. Rather than iterating, we use a simple numeric overestimate
788 * to determine k = floor(log10(d)). We scale relevant
789 * quantities using O(log2(k)) rather than O(k) multiplications.
790 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
791 * try to generate digits strictly left to right. Instead, we
792 * compute with fewer bits and propagate the carry if necessary
793 * when rounding the final digit up. This is often faster.
794 * 3. Under the assumption that input will be rounded nearest,
795 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
796 * That is, we allow equality in stopping tests when the
797 * round-nearest rule will give the same floating-point value
798 * as would satisfaction of the stopping test with strict
800 * 4. We remove common factors of powers of 2 from relevant
802 * 5. When converting floating-point integers less than 1e16,
803 * we use floating-point arithmetic rather than resorting
804 * to multiple-precision integers.
805 * 6. When asked to produce fewer than 15 digits, we first try
806 * to get by with floating-point arithmetic; we resort to
807 * multiple-precision integer arithmetic only if we cannot
808 * guarantee that the floating-point calculation has given
809 * the correctly rounded result. For k requested digits and
810 * "uniformly" distributed input, the probability is
811 * something like 10^(k-15) that we must resort to the int
816 dtoa(double d, int mode, int ndigits, int *decpt, int *sign, char **rve)
818 /* Arguments ndigits, decpt, sign are similar to those
819 of ecvt and fcvt; trailing zeros are suppressed from
820 the returned string. If not null, *rve is set to point
821 to the end of the return value. If d is +-Infinity or NaN,
822 then *decpt is set to 9999.
825 0 ==> shortest string that yields d when read in
826 and rounded to nearest.
827 1 ==> like 0, but with Steele & White stopping rule;
828 e.g. with IEEE P754 arithmetic , mode 0 gives
829 1e23 whereas mode 1 gives 9.999999999999999e22.
830 2 ==> max(1,ndigits) significant digits. This gives a
831 return value similar to that of ecvt, except
832 that trailing zeros are suppressed.
833 3 ==> through ndigits past the decimal point. This
834 gives a return value similar to that from fcvt,
835 except that trailing zeros are suppressed, and
836 ndigits can be negative.
837 4-9 should give the same return values as 2-3, i.e.,
838 4 <= mode <= 9 ==> same return as mode
839 2 + (mode & 1). These modes are mainly for
840 debugging; often they run slower but sometimes
841 faster than modes 2-3.
842 4,5,8,9 ==> left-to-right digit generation.
843 6-9 ==> don't try fast floating-point estimate
846 Values of mode other than 0-9 are treated as mode 0.
848 Sufficient space is allocated to the return value
849 to hold the suppressed trailing zeros.
852 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
853 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
854 spec_case, try_quick;
856 Bigint * b, *b1, *delta, *mlo=nil, *mhi, *S;
860 if (fpword0(d) & Sign_bit) {
861 /* set sign for everything, including 0's and NaNs */
863 fpword0(d) &= ~Sign_bit; /* clear sign bit */
867 if ((fpword0(d) & Exp_mask) == Exp_mask) {
868 /* Infinity or NaN */
870 if (!fpword1(d) && !(fpword0(d) & 0xfffff))
871 return nrv_alloc("Infinity", rve, 8);
872 return nrv_alloc("NaN", rve, 3);
876 return nrv_alloc("0", rve, 1);
879 b = d2b(d, &be, &bbits);
880 i = (int)(fpword0(d) >> Exp_shift1 & (Exp_mask >> Exp_shift1));
882 fpword0(d2) &= Frac_mask1;
883 fpword0(d2) |= Exp_11;
885 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
886 * log10(x) = log(x) / log(10)
887 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
888 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
890 * This suggests computing an approximation k to log10(d) by
892 * k = (i - Bias)*0.301029995663981
893 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
895 * We want k to be too large rather than too small.
896 * The error in the first-order Taylor series approximation
897 * is in our favor, so we just round up the constant enough
898 * to compensate for any error in the multiplication of
899 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
900 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
901 * adding 1e-13 to the constant term more than suffices.
902 * Hence we adjust the constant term to 0.1760912590558.
903 * (We could get a more accurate k by invoking log10,
904 * but this is probably not worthwhile.)
908 ds = (d2 - 1.5) * 0.289529654602168 + 0.1760912590558 + i * 0.301029995663981;
910 if (ds < 0. && ds != k)
911 k--; /* want k = floor(ds) */
913 if (k >= 0 && k <= Ten_pmax) {
935 if (mode < 0 || mode > 9)
956 ilim = ilim1 = i = ndigits;
968 s = s0 = rv_alloc(i);
970 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
972 /* Try to get by with floating-point arithmetic. */
978 ieps = 2; /* conservative */
983 /* prevent overflows */
985 d /= bigtens[n_bigtens-1];
988 for (; j; j >>= 1, i++)
994 } else if (j1 = -k) {
996 for (j = j1 >> 4; j; j >>= 1, i++)
1002 if (k_check && d < 1. && ilim > 0) {
1010 eps = ieps * d + 7.;
1011 fpword0(eps) -= (P - 1) * Exp_msk1;
1021 /* Generate ilim digits, then fix them up. */
1022 eps *= tens[ilim-1];
1023 for (i = 1; ; i++, d *= 10.) {
1026 *s++ = '0' + (int)L;
1030 else if (d < 0.5 - eps) {
1046 /* Do we have a "small" integer? */
1048 if (be >= 0 && k <= Int_max) {
1051 if (ndigits < 0 && ilim <= 0) {
1053 if (ilim < 0 || d <= 5 * ds)
1057 for (i = 1; ; i++) {
1060 *s++ = '0' + (int)L;
1063 if (d > ds || d == ds && L & 1) {
1097 if ((i = ilim) < 0) {
1106 if (m2 > 0 && s2 > 0) {
1107 i = m2 < s2 ? m2 : s2;
1115 mhi = pow5mult(mhi, m5);
1123 b = pow5mult(b, b5);
1127 S = pow5mult(S, s5);
1129 /* Check for special case that d is a normalized power of 2. */
1133 if (!fpword1(d) && !(fpword0(d) & Bndry_mask)
1135 /* The special case */
1142 /* Arrange for convenient computation of quotients:
1143 * shift left if necessary so divisor has 4 leading 0 bits.
1145 * Perhaps we should just compute leading 28 bits of S once
1146 * and for all and pass them and a shift to quorem, so it
1147 * can do shifts and ors to compute the numerator for q.
1149 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)
1167 if (cmp(b, S) < 0) {
1169 b = multadd(b, 10, 0); /* we botched the k estimate */
1171 mhi = multadd(mhi, 10, 0);
1175 if (ilim <= 0 && mode > 2) {
1176 if (ilim < 0 || cmp(b, S = multadd(S, 5, 0)) <= 0) {
1177 /* no digits, fcvt style */
1189 mhi = lshift(mhi, m2);
1191 /* Compute mlo -- check for special case
1192 * that d is a normalized power of 2.
1197 mhi = Balloc(mhi->k);
1199 mhi = lshift(mhi, Log2P);
1202 for (i = 1; ; i++) {
1203 dig = quorem(b, S) + '0';
1204 /* Do we yet have the shortest decimal string
1205 * that will round to d?
1208 delta = diff(S, mhi);
1209 j1 = delta->sign ? 1 : cmp(b, delta);
1211 if (j1 == 0 && !mode && !(fpword1(d) & 1)) {
1219 if (j < 0 || j == 0 && !mode
1220 && !(fpword1(d) & 1)
1225 if ((j1 > 0 || j1 == 0 && dig & 1)
1233 if (dig == '9') { /* possible if i == 1 */
1244 b = multadd(b, 10, 0);
1246 mlo = mhi = multadd(mhi, 10, 0);
1248 mlo = multadd(mlo, 10, 0);
1249 mhi = multadd(mhi, 10, 0);
1253 for (i = 1; ; i++) {
1254 *s++ = dig = quorem(b, S) + '0';
1257 b = multadd(b, 10, 0);
1260 /* Round off last digit */
1264 if (j > 0 || j == 0 && dig & 1) {
1281 if (mlo && mlo != mhi)