1 /****************************************************************
3 * The author of this software is David M. Gay.
5 * Copyright (c) 1991 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
34 _Jv_Bigint * b _AND _Jv_Bigint * S)
38 unsigned long carry, q, ys;
39 unsigned long *bx, *bxe, *sx, *sxe;
47 /*debug*/ if (b->_wds > n)
48 /*debug*/ Bug ("oversize b in quorem");
56 q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
59 /*debug*/ Bug ("oversized quotient in quorem");
69 ys = (si & 0xffff) * q + carry;
70 zs = (si >> 16) * q + (ys >> 16);
72 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
74 Sign_Extend (borrow, y);
75 z = (*bx >> 16) - (zs & 0xffff) + borrow;
77 Sign_Extend (borrow, z);
80 ys = *sx++ * q + carry;
82 y = *bx - (ys & 0xffff) + borrow;
84 Sign_Extend (borrow, y);
92 while (--bxe > bx && !*bxe)
108 ys = (si & 0xffff) + carry;
109 zs = (si >> 16) + (ys >> 16);
111 y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
113 Sign_Extend (borrow, y);
114 z = (*bx >> 16) - (zs & 0xffff) + borrow;
116 Sign_Extend (borrow, z);
121 y = *bx - (ys & 0xffff) + borrow;
123 Sign_Extend (borrow, y);
132 while (--bxe > bx && !*bxe)
144 print (_Jv_Bigint * b)
154 fprintf (stderr, "%08x", *x);
157 fprintf (stderr, "\n");
161 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
163 * Inspired by "How to Print Floating-Point Numbers Accurately" by
164 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
167 * 1. Rather than iterating, we use a simple numeric overestimate
168 * to determine k = floor(log10(d)). We scale relevant
169 * quantities using O(log2(k)) rather than O(k) multiplications.
170 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
171 * try to generate digits strictly left to right. Instead, we
172 * compute with fewer bits and propagate the carry if necessary
173 * when rounding the final digit up. This is often faster.
174 * 3. Under the assumption that input will be rounded nearest,
175 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
176 * That is, we allow equality in stopping tests when the
177 * round-nearest rule will give the same floating-point value
178 * as would satisfaction of the stopping test with strict
180 * 4. We remove common factors of powers of 2 from relevant
182 * 5. When converting floating-point integers less than 1e16,
183 * we use floating-point arithmetic rather than resorting
184 * to multiple-precision integers.
185 * 6. When asked to produce fewer than 15 digits, we first try
186 * to get by with floating-point arithmetic; we resort to
187 * multiple-precision integer arithmetic only if we cannot
188 * guarantee that the floating-point calculation has given
189 * the correctly rounded result. For k requested digits and
190 * "uniformly" distributed input, the probability is
191 * something like 10^(k-15) that we must resort to the long
198 (ptr, _d, mode, ndigits, decpt, sign, rve, float_type),
199 struct _Jv_reent *ptr _AND
209 float_type == 0 for double precision, 1 for float.
211 Arguments ndigits, decpt, sign are similar to those
212 of ecvt and fcvt; trailing zeros are suppressed from
213 the returned string. If not null, *rve is set to point
214 to the end of the return value. If d is +-Infinity or NaN,
215 then *decpt is set to 9999.
218 0 ==> shortest string that yields d when read in
219 and rounded to nearest.
220 1 ==> like 0, but with Steele & White stopping rule;
221 e.g. with IEEE P754 arithmetic , mode 0 gives
222 1e23 whereas mode 1 gives 9.999999999999999e22.
223 2 ==> max(1,ndigits) significant digits. This gives a
224 return value similar to that of ecvt, except
225 that trailing zeros are suppressed.
226 3 ==> through ndigits past the decimal point. This
227 gives a return value similar to that from fcvt,
228 except that trailing zeros are suppressed, and
229 ndigits can be negative.
230 4-9 should give the same return values as 2-3, i.e.,
231 4 <= mode <= 9 ==> same return as mode
232 2 + (mode & 1). These modes are mainly for
233 debugging; often they run slower but sometimes
234 faster than modes 2-3.
235 4,5,8,9 ==> left-to-right digit generation.
236 6-9 ==> don't try fast floating-point estimate
239 > 16 ==> Floating-point arg is treated as single precision.
241 Values of mode other than 0-9 are treated as mode 0.
243 Sufficient space is allocated to the return value
244 to hold the suppressed trailing zeros.
247 int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1, j, j1, k, k0,
248 k_check, leftright, m2, m5, s2, s5, spec_case, try_quick;
249 union double_union d, d2, eps;
251 #ifndef Sudden_Underflow
255 _Jv_Bigint *b, *b1, *delta, *mlo, *mhi, *S;
263 ptr->_result->_k = ptr->_result_k;
264 ptr->_result->_maxwds = 1 << ptr->_result_k;
265 Bfree (ptr, ptr->_result);
269 if (word0 (d) & Sign_bit)
271 /* set sign for everything, including 0's and NaNs */
273 word0 (d) &= ~Sign_bit; /* clear sign bit */
278 #if defined(IEEE_Arith) + defined(VAX)
280 if ((word0 (d) & Exp_mask) == Exp_mask)
282 if (word0 (d) == 0x8000)
285 /* Infinity or NaN */
289 !word1 (d) && !(word0 (d) & 0xfffff) ? "Infinity" :
302 d.d += 0; /* normalize */
313 b = d2b (ptr, d.d, &be, &bbits);
314 #ifdef Sudden_Underflow
315 i = (int) (word0 (d) >> Exp_shift1 & (Exp_mask >> Exp_shift1));
317 if ((i = (int) (word0 (d) >> Exp_shift1 & (Exp_mask >> Exp_shift1))))
321 word0 (d2) &= Frac_mask1;
322 word0 (d2) |= Exp_11;
324 if (j = 11 - hi0bits (word0 (d2) & Frac_mask))
328 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
329 * log10(x) = log(x) / log(10)
330 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
331 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
333 * This suggests computing an approximation k to log10(d) by
335 * k = (i - Bias)*0.301029995663981
336 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
338 * We want k to be too large rather than too small.
339 * The error in the first-order Taylor series approximation
340 * is in our favor, so we just round up the constant enough
341 * to compensate for any error in the multiplication of
342 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
343 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
344 * adding 1e-13 to the constant term more than suffices.
345 * Hence we adjust the constant term to 0.1760912590558.
346 * (We could get a more accurate k by invoking log10,
347 * but this is probably not worthwhile.)
355 #ifndef Sudden_Underflow
360 /* d is denormalized */
362 i = bbits + be + (Bias + (P - 1) - 1);
363 x = i > 32 ? word0 (d) << (64 - i) | word1 (d) >> (i - 32)
364 : word1 (d) << (32 - i);
366 word0 (d2) -= 31 * Exp_msk1; /* adjust exponent */
367 i -= (Bias + (P - 1) - 1) + 1;
371 ds = (d2.d - 1.5) * 0.289529654602168 + 0.1760912590558 + i * 0.301029995663981;
373 if (ds < 0. && ds != k)
374 k--; /* want k = floor(ds) */
376 if (k >= 0 && k <= Ten_pmax)
405 if (mode < 0 || mode > 9)
428 ilim = ilim1 = i = ndigits;
440 j = sizeof (unsigned long);
441 for (ptr->_result_k = 0; (int) (sizeof (_Jv_Bigint) - sizeof (unsigned long)) + j <= i;
444 ptr->_result = Balloc (ptr, ptr->_result_k);
445 s = s0 = (char *) ptr->_result;
447 if (ilim >= 0 && ilim <= Quick_max && try_quick)
449 /* Try to get by with floating-point arithmetic. */
455 ieps = 2; /* conservative */
462 /* prevent overflows */
464 d.d /= bigtens[n_bigtens - 1];
467 for (; j; j >>= 1, i++)
477 d.d *= tens[j1 & 0xf];
478 for (j = j1 >> 4; j; j >>= 1, i++)
485 if (k_check && d.d < 1. && ilim > 0)
494 eps.d = ieps * d.d + 7.;
495 word0 (eps) -= (P - 1) * Exp_msk1;
509 /* Use Steele & White method of only
510 * generating digits needed.
512 eps.d = 0.5 / tens[ilim - 1] - eps.d;
517 *s++ = '0' + (int) L;
520 if (1. - d.d < eps.d)
531 /* Generate ilim digits, then fix them up. */
532 eps.d *= tens[ilim - 1];
533 for (i = 1;; i++, d.d *= 10.)
537 *s++ = '0' + (int) L;
540 if (d.d > 0.5 + eps.d)
542 else if (d.d < 0.5 - eps.d)
561 /* Do we have a "small" integer? */
563 if (be >= 0 && k <= Int_max)
567 if (ndigits < 0 && ilim <= 0)
570 if (ilim < 0 || d.d <= 5 * ds)
578 #ifdef Check_FLT_ROUNDS
579 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
586 *s++ = '0' + (int) L;
590 if (d.d > ds || (d.d == ds && L & 1))
618 #ifndef Sudden_Underflow
619 denorm ? be + (Bias + (P - 1) - 1 + 1) :
622 1 + 4 * P - 3 - bbits + ((bbits + be - 1) & 3);
648 if (m2 > 0 && s2 > 0)
650 i = m2 < s2 ? m2 : s2;
661 mhi = pow5mult (ptr, mhi, m5);
662 b1 = mult (ptr, mhi, b);
667 b = pow5mult (ptr, b, j);
670 b = pow5mult (ptr, b, b5);
674 S = pow5mult (ptr, S, s5);
676 /* Check for special case that d is a normalized power of 2. */
680 if (!word1 (d) && !(word0 (d) & Bndry_mask)
681 #ifndef Sudden_Underflow
682 && word0 (d) & Exp_mask
686 /* The special case */
695 /* Arrange for convenient computation of quotients:
696 * shift left if necessary so divisor has 4 leading 0 bits.
698 * Perhaps we should just compute leading 28 bits of S once
699 * and for all and pass them and a shift to quorem, so it
700 * can do shifts and ors to compute the numerator for q.
704 if ((i = ((s5 ? 32 - hi0bits (S->_x[S->_wds - 1]) : 1) + s2) & 0x1f))
707 if ((i = ((s5 ? 32 - hi0bits (S->_x[S->_wds - 1]) : 1) + s2) & 0xf))
725 b = lshift (ptr, b, b2);
727 S = lshift (ptr, S, s2);
733 b = multadd (ptr, b, 10, 0); /* we botched the k estimate */
735 mhi = multadd (ptr, mhi, 10, 0);
739 if (ilim <= 0 && mode > 2)
741 if (ilim < 0 || cmp (b, S = multadd (ptr, S, 5, 0)) <= 0)
743 /* no digits, fcvt style */
756 mhi = lshift (ptr, mhi, m2);
758 /* Single precision case, */
760 mhi = lshift (ptr, mhi, 29);
762 /* Compute mlo -- check for special case
763 * that d is a normalized power of 2.
769 mhi = Balloc (ptr, mhi->_k);
771 mhi = lshift (ptr, mhi, Log2P);
776 dig = quorem (b, S) + '0';
777 /* Do we yet have the shortest decimal string
778 * that will round to d?
781 delta = diff (ptr, S, mhi);
782 j1 = delta->_sign ? 1 : cmp (b, delta);
785 if (j1 == 0 && !mode && !(word1 (d) & 1))
795 if (j < 0 || (j == 0 && !mode
803 b = lshift (ptr, b, 1);
805 if ((j1 > 0 || (j1 == 0 && dig & 1))
815 { /* possible if i == 1 */
826 b = multadd (ptr, b, 10, 0);
828 mlo = mhi = multadd (ptr, mhi, 10, 0);
831 mlo = multadd (ptr, mlo, 10, 0);
832 mhi = multadd (ptr, mhi, 10, 0);
839 *s++ = dig = quorem (b, S) + '0';
842 b = multadd (ptr, b, 10, 0);
845 /* Round off last digit */
847 b = lshift (ptr, b, 1);
849 if (j > 0 || (j == 0 && dig & 1))
870 if (mlo && mlo != mhi)
886 (_d, mode, ndigits, decpt, sign, rve, buf, float_type),
896 struct _Jv_reent reent;
898 memset (&reent, 0, sizeof reent);
900 p = _dtoa_r (&reent, _d, mode, ndigits, decpt, sign, rve, float_type);