1 /****************************************************************
3 * The author of this software is David M. Gay.
5 * Copyright (c) 1991, 2006 by AT&T.
7 * Permission to use, copy, modify, and distribute this software for any
8 * purpose without fee is hereby granted, provided that this entire notice
9 * is included in all copies of any software which is or includes a copy
10 * or modification of this software and in all copies of the supporting
11 * documentation for such software.
13 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
14 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
15 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
16 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
18 ***************************************************************/
20 /* Please send bug reports to
22 AT&T Bell Laboratories, Room 2C-463
24 Murray Hill, NJ 07974-2070
26 dmg@research.att.com or research!dmg
35 _Jv_Bigint * b _AND _Jv_Bigint * S)
39 unsigned long carry, q, ys;
40 unsigned long *bx, *bxe, *sx, *sxe;
48 /*debug*/ if (b->_wds > n)
49 /*debug*/ Bug ("oversize b in quorem");
57 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
60 /*debug*/ Bug ("oversized quotient in quorem");
70 ys = (si & 0xffff) * q + carry;
71 zs = (si >> 16) * q + (ys >> 16);
73 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
75 Sign_Extend (borrow, y);
76 z = (*bx >> 16) - (zs & 0xffff) + borrow;
78 Sign_Extend (borrow, z);
81 ys = *sx++ * q + carry;
83 y = *bx - (ys & 0xffff) + borrow;
85 Sign_Extend (borrow, y);
93 while (--bxe > bx && !*bxe)
109 ys = (si & 0xffff) + carry;
110 zs = (si >> 16) + (ys >> 16);
112 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
114 Sign_Extend (borrow, y);
115 z = (*bx >> 16) - (zs & 0xffff) + borrow;
117 Sign_Extend (borrow, z);
122 y = *bx - (ys & 0xffff) + borrow;
124 Sign_Extend (borrow, y);
133 while (--bxe > bx && !*bxe)
145 print (_Jv_Bigint * b)
155 fprintf (stderr, "%08x", *x);
158 fprintf (stderr, "\n");
162 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
164 * Inspired by "How to Print Floating-Point Numbers Accurately" by
165 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
168 * 1. Rather than iterating, we use a simple numeric overestimate
169 * to determine k = floor(log10(d)). We scale relevant
170 * quantities using O(log2(k)) rather than O(k) multiplications.
171 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
172 * try to generate digits strictly left to right. Instead, we
173 * compute with fewer bits and propagate the carry if necessary
174 * when rounding the final digit up. This is often faster.
175 * 3. Under the assumption that input will be rounded nearest,
176 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
177 * That is, we allow equality in stopping tests when the
178 * round-nearest rule will give the same floating-point value
179 * as would satisfaction of the stopping test with strict
181 * 4. We remove common factors of powers of 2 from relevant
183 * 5. When converting floating-point integers less than 1e16,
184 * we use floating-point arithmetic rather than resorting
185 * to multiple-precision integers.
186 * 6. When asked to produce fewer than 15 digits, we first try
187 * to get by with floating-point arithmetic; we resort to
188 * multiple-precision integer arithmetic only if we cannot
189 * guarantee that the floating-point calculation has given
190 * the correctly rounded result. For k requested digits and
191 * "uniformly" distributed input, the probability is
192 * something like 10^(k-15) that we must resort to the long
199 (ptr, _d, mode, ndigits, decpt, sign, rve, float_type),
200 struct _Jv_reent *ptr _AND
210 float_type == 0 for double precision, 1 for float.
212 Arguments ndigits, decpt, sign are similar to those
213 of ecvt and fcvt; trailing zeros are suppressed from
214 the returned string. If not null, *rve is set to point
215 to the end of the return value. If d is +-Infinity or NaN,
216 then *decpt is set to 9999.
219 0 ==> shortest string that yields d when read in
220 and rounded to nearest.
221 1 ==> like 0, but with Steele & White stopping rule;
222 e.g. with IEEE P754 arithmetic , mode 0 gives
223 1e23 whereas mode 1 gives 9.999999999999999e22.
224 2 ==> max(1,ndigits) significant digits. This gives a
225 return value similar to that of ecvt, except
226 that trailing zeros are suppressed.
227 3 ==> through ndigits past the decimal point. This
228 gives a return value similar to that from fcvt,
229 except that trailing zeros are suppressed, and
230 ndigits can be negative.
231 4-9 should give the same return values as 2-3, i.e.,
232 4 <= mode <= 9 ==> same return as mode
233 2 + (mode & 1). These modes are mainly for
234 debugging; often they run slower but sometimes
235 faster than modes 2-3.
236 4,5,8,9 ==> left-to-right digit generation.
237 6-9 ==> don't try fast floating-point estimate
240 > 16 ==> Floating-point arg is treated as single precision.
242 Values of mode other than 0-9 are treated as mode 0.
244 Sufficient space is allocated to the return value
245 to hold the suppressed trailing zeros.
248 int bbits, b2, b5, be, dig, i, ieps, ilim0, j, j1, k, k0,
249 k_check, leftright, m2, m5, s2, s5, try_quick;
250 int ilim = 0, ilim1 = 0, spec_case = 0;
251 union double_union d, d2, eps;
253 #ifndef Sudden_Underflow
257 _Jv_Bigint *b, *b1, *delta, *mlo = NULL, *mhi, *S;
265 ptr->_result->_k = ptr->_result_k;
266 ptr->_result->_maxwds = 1 << ptr->_result_k;
267 Bfree (ptr, ptr->_result);
271 if (word0 (d) & Sign_bit)
273 /* set sign for everything, including 0's and NaNs */
275 word0 (d) &= ~Sign_bit; /* clear sign bit */
280 #if defined(IEEE_Arith) + defined(VAX)
282 if ((word0 (d) & Exp_mask) == Exp_mask)
284 if (word0 (d) == 0x8000)
287 /* Infinity or NaN */
291 !word1 (d) && !(word0 (d) & 0xfffff) ? "Infinity" :
304 d.d += 0; /* normalize */
315 b = d2b (ptr, d.d, &be, &bbits);
316 #ifdef Sudden_Underflow
317 i = (int) (word0 (d) >> Exp_shift1 & (Exp_mask >> Exp_shift1));
319 if ((i = (int) (word0 (d) >> Exp_shift1 & (Exp_mask >> Exp_shift1))))
323 word0 (d2) &= Frac_mask1;
324 word0 (d2) |= Exp_11;
326 if (j = 11 - hi0bits (word0 (d2) & Frac_mask))
330 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
331 * log10(x) = log(x) / log(10)
332 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
333 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
335 * This suggests computing an approximation k to log10(d) by
337 * k = (i - Bias)*0.301029995663981
338 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
340 * We want k to be too large rather than too small.
341 * The error in the first-order Taylor series approximation
342 * is in our favor, so we just round up the constant enough
343 * to compensate for any error in the multiplication of
344 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
345 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
346 * adding 1e-13 to the constant term more than suffices.
347 * Hence we adjust the constant term to 0.1760912590558.
348 * (We could get a more accurate k by invoking log10,
349 * but this is probably not worthwhile.)
357 #ifndef Sudden_Underflow
362 /* d is denormalized */
364 i = bbits + be + (Bias + (P - 1) - 1);
365 x = i > 32 ? word0 (d) << (64 - i) | word1 (d) >> (i - 32)
366 : word1 (d) << (32 - i);
368 word0 (d2) -= 31 * Exp_msk1; /* adjust exponent */
369 i -= (Bias + (P - 1) - 1) + 1;
373 ds = (d2.d - 1.5) * 0.289529654602168 + 0.1760912590558 + i * 0.301029995663981;
375 if (ds < 0. && ds != k)
376 k--; /* want k = floor(ds) */
378 if (k >= 0 && k <= Ten_pmax)
407 if (mode < 0 || mode > 9)
430 ilim = ilim1 = i = ndigits;
442 j = sizeof (unsigned long);
443 for (ptr->_result_k = 0; (int) (sizeof (_Jv_Bigint) - sizeof (unsigned long)) + j <= i;
446 ptr->_result = Balloc (ptr, ptr->_result_k);
447 s = s0 = (char *) ptr->_result;
449 if (ilim >= 0 && ilim <= Quick_max && try_quick)
451 /* Try to get by with floating-point arithmetic. */
457 ieps = 2; /* conservative */
464 /* prevent overflows */
466 d.d /= bigtens[n_bigtens - 1];
469 for (; j; j >>= 1, i++)
479 d.d *= tens[j1 & 0xf];
480 for (j = j1 >> 4; j; j >>= 1, i++)
487 if (k_check && d.d < 1. && ilim > 0)
496 eps.d = ieps * d.d + 7.;
497 word0 (eps) -= (P - 1) * Exp_msk1;
511 /* Use Steele & White method of only
512 * generating digits needed.
514 eps.d = 0.5 / tens[ilim - 1] - eps.d;
519 *s++ = '0' + (int) L;
522 if (1. - d.d < eps.d)
533 /* Generate ilim digits, then fix them up. */
534 eps.d *= tens[ilim - 1];
535 for (i = 1;; i++, d.d *= 10.)
539 *s++ = '0' + (int) L;
542 if (d.d > 0.5 + eps.d)
544 else if (d.d < 0.5 - eps.d)
563 /* Do we have a "small" integer? */
565 if (be >= 0 && k <= Int_max)
569 if (ndigits < 0 && ilim <= 0)
572 if (ilim < 0 || d.d <= 5 * ds)
580 #ifdef Check_FLT_ROUNDS
581 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
588 *s++ = '0' + (int) L;
592 if (d.d > ds || (d.d == ds && L & 1))
620 #ifndef Sudden_Underflow
621 denorm ? be + (Bias + (P - 1) - 1 + 1) :
624 1 + 4 * P - 3 - bbits + ((bbits + be - 1) & 3);
650 if (m2 > 0 && s2 > 0)
652 i = m2 < s2 ? m2 : s2;
663 mhi = pow5mult (ptr, mhi, m5);
664 b1 = mult (ptr, mhi, b);
669 b = pow5mult (ptr, b, j);
672 b = pow5mult (ptr, b, b5);
676 S = pow5mult (ptr, S, s5);
678 /* Check for special case that d is a normalized power of 2. */
682 if (!word1 (d) && !(word0 (d) & Bndry_mask)
683 #ifndef Sudden_Underflow
684 && word0(d) & Exp_mask
688 /* The special case */
697 /* Arrange for convenient computation of quotients:
698 * shift left if necessary so divisor has 4 leading 0 bits.
700 * Perhaps we should just compute leading 28 bits of S once
701 * and for all and pass them and a shift to quorem, so it
702 * can do shifts and ors to compute the numerator for q.
706 if ((i = ((s5 ? 32 - hi0bits (S->_x[S->_wds - 1]) : 1) + s2) & 0x1f))
709 if ((i = ((s5 ? 32 - hi0bits (S->_x[S->_wds - 1]) : 1) + s2) & 0xf))
727 b = lshift (ptr, b, b2);
729 S = lshift (ptr, S, s2);
735 b = multadd (ptr, b, 10, 0); /* we botched the k estimate */
737 mhi = multadd (ptr, mhi, 10, 0);
741 if (ilim <= 0 && mode > 2)
743 if (ilim < 0 || cmp (b, S = multadd (ptr, S, 5, 0)) <= 0)
745 /* no digits, fcvt style */
758 mhi = lshift (ptr, mhi, m2);
760 /* Single precision case, */
762 mhi = lshift (ptr, mhi, 29);
764 /* Compute mlo -- check for special case
765 * that d is a normalized power of 2.
771 mhi = Balloc (ptr, mhi->_k);
773 mhi = lshift (ptr, mhi, Log2P);
778 dig = quorem (b, S) + '0';
779 /* Do we yet have the shortest decimal string
780 * that will round to d?
783 delta = diff (ptr, S, mhi);
784 j1 = delta->_sign ? 1 : cmp (b, delta);
787 if (j1 == 0 && !mode && !(word1 (d) & 1))
797 if (j < 0 || (j == 0 && !mode
805 b = lshift (ptr, b, 1);
807 if ((j1 > 0 || (j1 == 0 && dig & 1))
817 { /* possible if i == 1 */
828 b = multadd (ptr, b, 10, 0);
830 mlo = mhi = multadd (ptr, mhi, 10, 0);
833 mlo = multadd (ptr, mlo, 10, 0);
834 mhi = multadd (ptr, mhi, 10, 0);
841 *s++ = dig = quorem (b, S) + '0';
844 b = multadd (ptr, b, 10, 0);
847 /* Round off last digit */
849 b = lshift (ptr, b, 1);
851 if (j > 0 || (j == 0 && dig & 1))
872 if (mlo && mlo != mhi)
888 (_d, mode, ndigits, decpt, sign, rve, buf, float_type),
898 struct _Jv_reent reent;
902 memset (&reent, 0, sizeof reent);
904 p = _dtoa_r (&reent, _d, mode, ndigits, decpt, sign, rve, float_type);
907 for (i = 0; i < reent._result_k; ++i)
909 struct _Jv_Bigint *l = reent._freelist[i];
912 struct _Jv_Bigint *next = l->_next;
918 free (reent._freelist);