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) ((x).hi)
25 #define fpword1(x) ((x).lo)
27 /* Ten_pmax = floor(P*log(2)/log(5)) */
28 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
29 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
30 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
34 #define Exp_msk1 0x100000
35 #define Exp_msk11 0x100000
36 #define Exp_mask 0x7ff00000
40 #define Exp_1 0x3ff00000
41 #define Exp_11 0x3ff00000
43 #define Frac_mask 0xfffff
44 #define Frac_mask1 0xfffff
47 #define Bndry_mask 0xfffff
48 #define Bndry_mask1 0xfffff
49 #define Sign_bit 0x80000000
56 #define rounded_product(a,b) a *= b
57 #define rounded_quotient(a,b) a /= b
59 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
60 #define Big1 0xffffffff
62 #define FFFFFFFF 0xffffffffUL
71 int k, maxwds, sign, wds;
75 typedef struct Bigint Bigint;
77 static Bigint *freelist[Kmax+1];
86 assert(k < nelem(freelist));
89 if (rv = freelist[k]) {
90 freelist[k] = rv->next;
93 len = (sizeof(Bigint) + (x - 1) * sizeof(unsigned int) + sizeof(double) -1)
95 if (pmem_next - private_mem + len <= PRIVATE_mem) {
96 rv = (Bigint * )pmem_next;
99 rv = (Bigint * )malloc(len * sizeof(double));
104 rv->sign = rv->wds = 0;
112 ACQUIRE_DTOA_LOCK(0);
113 v->next = freelist[v->k];
119 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
120 y->wds*sizeof(int) + 2*sizeof(int))
123 multadd(Bigint *b, int m, int a) /* multiply by m and add a */
126 unsigned int carry, *x, y;
136 y = (xi & 0xffff) * m + carry;
137 z = (xi >> 16) * m + (y >> 16);
139 *x++ = (z << 16) + (y & 0xffff);
142 if (wds >= b->maxwds) {
143 b1 = Balloc(b->k + 1);
155 s2b(const char *s, int nd0, int nd, unsigned int y9)
162 for (k = 0, y = 1; x > y; y <<= 1, k++)
172 b = multadd(b, 10, *s++ - '0');
178 b = multadd(b, 10, *s++ - '0');
183 hi0bits(register unsigned int x)
187 if (!(x & 0xffff0000)) {
191 if (!(x & 0xff000000)) {
195 if (!(x & 0xf0000000)) {
199 if (!(x & 0xc0000000)) {
203 if (!(x & 0x80000000)) {
205 if (!(x & 0x40000000))
212 lo0bits(unsigned int *y)
215 register unsigned int x = *y;
266 mult(Bigint *a, Bigint *b)
270 unsigned int * x, *xa, *xae, *xb, *xbe, *xc, *xc0;
272 unsigned int carry, z;
275 if (a->wds < b->wds) {
287 for (x = c->x, xa = x + wc; x < xa; x++)
294 for (; xb < xbe; xb++, xc0++) {
295 if (y = *xb & 0xffff) {
300 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
302 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
314 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
317 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
323 for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc)
332 pow5mult(Bigint *b, int k)
334 Bigint * b1, *p5, *p51;
336 static int p05[3] = {
340 b = multadd(b, p05[i-1], 0);
346 ACQUIRE_DTOA_LOCK(1);
361 if (!(p51 = p5->next)) {
362 ACQUIRE_DTOA_LOCK(1);
363 if (!(p51 = p5->next)) {
364 p51 = p5->next = mult(p5, p5);
375 lshift(Bigint *b, int k)
379 unsigned int * x, *x1, *xe, z;
384 for (i = b->maxwds; n1 > i; i <<= 1)
388 for (i = 0; i < n; i++)
411 cmp(Bigint *a, Bigint *b)
413 unsigned int * xa, *xa0, *xb, *xb0;
426 return * xa < *xb ? -1 : 1;
434 diff(Bigint *a, Bigint *b)
438 unsigned int * xa, *xae, *xb, *xbe, *xc;
439 unsigned int borrow, y;
467 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
468 borrow = (y & 0x10000) >> 16;
469 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
470 borrow = (z & 0x10000) >> 16;
474 y = (*xa & 0xffff) - borrow;
475 borrow = (y & 0x10000) >> 16;
476 z = (*xa++ >> 16) - borrow;
477 borrow = (z & 0x10000) >> 16;
492 L = (fpword0(x) & Exp_mask) - (P - 1) * Exp_msk1;
499 b2d(Bigint *a, int *e)
501 unsigned int * xa, *xa0, w, y, z;
503 #define d0 fpword0(d)
504 #define d1 fpword1(d)
513 d0 = Exp_1 | y >> Ebits - k;
514 w = xa > xa0 ? *--xa : 0;
515 d1 = y << (32 - Ebits) + k | w >> Ebits - k;
518 z = xa > xa0 ? *--xa : 0;
520 d0 = Exp_1 | y << k | z >> 32 - k;
521 y = xa > xa0 ? *--xa : 0;
522 d1 = z << k | y >> 32 - k;
534 d2b(FPdbleword d, int *e, int *bits)
538 unsigned int * x, y, z;
546 d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
547 de = (int)(d0 >> Exp_shift);
550 if (k = lo0bits(&y)) {
551 x[0] = y | z << 32 - k;
555 b->wds = (x[1] = z) ? 2 : 1;
562 *e = de - Bias - (P - 1) + k;
571 ratio(Bigint *a, Bigint *b)
578 k = ka - kb + 32 * (a->wds - b->wds);
580 fpword0(da) += k * Exp_msk1;
583 fpword0(db) += k * Exp_msk1;
590 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
591 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
597 1e16, 1e32, 1e64, 1e128, 1e256 };
599 static const double tinytens[] = {
600 1e-16, 1e-32, 1e-64, 1e-128,
601 9007199254740992.e-256
604 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
605 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
606 #define Scale_Bit 0x10
609 #define NAN_WORD0 0x7ff80000
614 match(const char **sp, char *t)
617 const char * s = *sp;
620 if ((c = *++s) >= 'A' && c <= 'Z')
630 gethex(FPdbleword *rvp, const char **sp)
632 unsigned int c, x[2];
634 int havedig, udx0, xshift;
637 havedig = xshift = 0;
640 while (c = *(const unsigned char * )++s) {
641 if (c >= '0' && c <= '9')
643 else if (c >= 'a' && c <= 'f')
645 else if (c >= 'A' && c <= 'F')
648 if (udx0 && havedig) {
653 } else if (/*(*/ c == ')') {
657 return; /* invalid form: don't change *sp */
665 x[0] = (x[0] << 4) | (x[1] >> 28);
666 x[1] = (x[1] << 4) | c;
668 if ((x[0] &= 0xfffff) || x[1]) {
669 fpword0(*rvp) = Exp_mask | x[0];
670 fpword1(*rvp) = x[1];
675 quorem(Bigint *b, Bigint *S)
678 unsigned int * bx, *bxe, q, *sx, *sxe;
679 unsigned int borrow, carry, y, ys;
680 unsigned int si, z, zs;
689 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
695 ys = (si & 0xffff) * q + carry;
696 zs = (si >> 16) * q + (ys >> 16);
698 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
699 borrow = (y & 0x10000) >> 16;
700 z = (*bx >> 16) - (zs & 0xffff) - borrow;
701 borrow = (z & 0x10000) >> 16;
706 while (--bxe > bx && !*bxe)
711 if (cmp(b, S) >= 0) {
719 ys = (si & 0xffff) + carry;
720 zs = (si >> 16) + (ys >> 16);
722 y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
723 borrow = (y & 0x10000) >> 16;
724 z = (*bx >> 16) - (zs & 0xffff) - borrow;
725 borrow = (z & 0x10000) >> 16;
731 while (--bxe > bx && !*bxe)
744 j = sizeof(unsigned int);
746 sizeof(Bigint) - sizeof(unsigned int) - sizeof(int) + j <= i;
749 r = (int * )Balloc(k);
756 nrv_alloc(char *s, char **rve, int n)
760 t = rv = rv_alloc(n);
768 /* freedtoa(s) must be used to free values s returned by dtoa
769 * when MULTIPLE_THREADS is #defined. It should be used in all cases,
770 * but for consistency with earlier versions of dtoa, it is optional
771 * when MULTIPLE_THREADS is not defined.
777 Bigint * b = (Bigint * )((int *)s - 1);
778 b->maxwds = 1 << (b->k = *(int * )b);
782 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
784 * Inspired by "How to Print Floating-Point Numbers Accurately" by
785 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
788 * 1. Rather than iterating, we use a simple numeric overestimate
789 * to determine k = floor(log10(d)). We scale relevant
790 * quantities using O(log2(k)) rather than O(k) multiplications.
791 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
792 * try to generate digits strictly left to right. Instead, we
793 * compute with fewer bits and propagate the carry if necessary
794 * when rounding the final digit up. This is often faster.
795 * 3. Under the assumption that input will be rounded nearest,
796 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
797 * That is, we allow equality in stopping tests when the
798 * round-nearest rule will give the same floating-point value
799 * as would satisfaction of the stopping test with strict
801 * 4. We remove common factors of powers of 2 from relevant
803 * 5. When converting floating-point integers less than 1e16,
804 * we use floating-point arithmetic rather than resorting
805 * to multiple-precision integers.
806 * 6. When asked to produce fewer than 15 digits, we first try
807 * to get by with floating-point arithmetic; we resort to
808 * multiple-precision integer arithmetic only if we cannot
809 * guarantee that the floating-point calculation has given
810 * the correctly rounded result. For k requested digits and
811 * "uniformly" distributed input, the probability is
812 * something like 10^(k-15) that we must resort to the int
817 dtoa(double _d, int mode, int ndigits, int *decpt, int *sign, char **rve)
819 /* Arguments ndigits, decpt, sign are similar to those
820 of ecvt and fcvt; trailing zeros are suppressed from
821 the returned string. If not null, *rve is set to point
822 to the end of the return value. If d is +-Infinity or NaN,
823 then *decpt is set to 9999.
826 0 ==> shortest string that yields d when read in
827 and rounded to nearest.
828 1 ==> like 0, but with Steele & White stopping rule;
829 e.g. with IEEE P754 arithmetic , mode 0 gives
830 1e23 whereas mode 1 gives 9.999999999999999e22.
831 2 ==> max(1,ndigits) significant digits. This gives a
832 return value similar to that of ecvt, except
833 that trailing zeros are suppressed.
834 3 ==> through ndigits past the decimal point. This
835 gives a return value similar to that from fcvt,
836 except that trailing zeros are suppressed, and
837 ndigits can be negative.
838 4-9 should give the same return values as 2-3, i.e.,
839 4 <= mode <= 9 ==> same return as mode
840 2 + (mode & 1). These modes are mainly for
841 debugging; often they run slower but sometimes
842 faster than modes 2-3.
843 4,5,8,9 ==> left-to-right digit generation.
844 6-9 ==> don't try fast floating-point estimate
847 Values of mode other than 0-9 are treated as mode 0.
849 Sufficient space is allocated to the return value
850 to hold the suppressed trailing zeros.
853 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
854 j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
855 spec_case, try_quick;
857 Bigint * b, *b1, *delta, *mlo=nil, *mhi, *S;
859 FPdbleword d, d2, eps;
863 if (fpword0(d) & Sign_bit) {
864 /* set sign for everything, including 0's and NaNs */
866 fpword0(d) &= ~Sign_bit; /* clear sign bit */
870 if ((fpword0(d) & Exp_mask) == Exp_mask) {
871 /* Infinity or NaN */
873 if (!fpword1(d) && !(fpword0(d) & 0xfffff))
874 return nrv_alloc("Infinity", rve, 8);
875 return nrv_alloc("NaN", rve, 3);
879 return nrv_alloc("0", rve, 1);
882 b = d2b(d, &be, &bbits);
883 i = (int)(fpword0(d) >> Exp_shift1 & (Exp_mask >> Exp_shift1));
885 fpword0(d2) &= Frac_mask1;
886 fpword0(d2) |= Exp_11;
888 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
889 * log10(x) = log(x) / log(10)
890 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
891 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
893 * This suggests computing an approximation k to log10(d) by
895 * k = (i - Bias)*0.301029995663981
896 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
898 * We want k to be too large rather than too small.
899 * The error in the first-order Taylor series approximation
900 * is in our favor, so we just round up the constant enough
901 * to compensate for any error in the multiplication of
902 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
903 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
904 * adding 1e-13 to the constant term more than suffices.
905 * Hence we adjust the constant term to 0.1760912590558.
906 * (We could get a more accurate k by invoking log10,
907 * but this is probably not worthwhile.)
911 ds = (d2.x - 1.5) * 0.289529654602168 + 0.1760912590558 + i * 0.301029995663981;
913 if (ds < 0. && ds != k)
914 k--; /* want k = floor(ds) */
916 if (k >= 0 && k <= Ten_pmax) {
939 if (mode < 0 || mode > 9)
961 ilim = ilim1 = i = ndigits;
973 s = s0 = rv_alloc(i);
975 if (ilim >= 0 && ilim <= Quick_max && try_quick) {
977 /* Try to get by with floating-point arithmetic. */
983 ieps = 2; /* conservative */
988 /* prevent overflows */
990 d.x /= bigtens[n_bigtens-1];
993 for (; j; j >>= 1, i++)
999 } else if (j1 = -k) {
1000 d.x *= tens[j1 & 0xf];
1001 for (j = j1 >> 4; j; j >>= 1, i++)
1007 if (k_check && d.x < 1. && ilim > 0) {
1015 eps.x = ieps * d.x + 7.;
1016 fpword0(eps) -= (P - 1) * Exp_msk1;
1026 /* Generate ilim digits, then fix them up. */
1027 eps.x *= tens[ilim-1];
1028 for (i = 1; ; i++, d.x *= 10.) {
1031 *s++ = '0' + (int)L;
1033 if (d.x > 0.5 + eps.x)
1035 else if (d.x < 0.5 - eps.x) {
1051 /* Do we have a "small" integer? */
1053 if (be >= 0 && k <= Int_max) {
1056 if (ndigits < 0 && ilim <= 0) {
1058 if (ilim < 0 || d.x <= 5 * ds)
1062 for (i = 1; ; i++) {
1065 *s++ = '0' + (int)L;
1068 if (d.x > ds || d.x == ds && L & 1) {
1102 if ((i = ilim) < 0) {
1111 if (m2 > 0 && s2 > 0) {
1112 i = m2 < s2 ? m2 : s2;
1120 mhi = pow5mult(mhi, m5);
1128 b = pow5mult(b, b5);
1132 S = pow5mult(S, s5);
1134 /* Check for special case that d is a normalized power of 2. */
1138 if (!fpword1(d) && !(fpword0(d) & Bndry_mask)
1140 /* The special case */
1147 /* Arrange for convenient computation of quotients:
1148 * shift left if necessary so divisor has 4 leading 0 bits.
1150 * Perhaps we should just compute leading 28 bits of S once
1151 * and for all and pass them and a shift to quorem, so it
1152 * can do shifts and ors to compute the numerator for q.
1154 if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)
1172 if (cmp(b, S) < 0) {
1174 b = multadd(b, 10, 0); /* we botched the k estimate */
1176 mhi = multadd(mhi, 10, 0);
1180 if (ilim <= 0 && mode > 2) {
1181 if (ilim < 0 || cmp(b, S = multadd(S, 5, 0)) <= 0) {
1182 /* no digits, fcvt style */
1194 mhi = lshift(mhi, m2);
1196 /* Compute mlo -- check for special case
1197 * that d is a normalized power of 2.
1202 mhi = Balloc(mhi->k);
1204 mhi = lshift(mhi, Log2P);
1207 for (i = 1; ; i++) {
1208 dig = quorem(b, S) + '0';
1209 /* Do we yet have the shortest decimal string
1210 * that will round to d?
1213 delta = diff(S, mhi);
1214 j1 = delta->sign ? 1 : cmp(b, delta);
1216 if (j1 == 0 && !mode && !(fpword1(d) & 1)) {
1224 if (j < 0 || j == 0 && !mode
1225 && !(fpword1(d) & 1)
1230 if ((j1 > 0 || j1 == 0 && dig & 1)
1238 if (dig == '9') { /* possible if i == 1 */
1249 b = multadd(b, 10, 0);
1251 mlo = mhi = multadd(mhi, 10, 0);
1253 mlo = multadd(mlo, 10, 0);
1254 mhi = multadd(mhi, 10, 0);
1258 for (i = 1; ; i++) {
1259 *s++ = dig = quorem(b, S) + '0';
1262 b = multadd(b, 10, 0);
1265 /* Round off last digit */
1269 if (j > 0 || j == 0 && dig & 1) {
1286 if (mlo && mlo != mhi)