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
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
, ilim
, ilim0
, ilim1
, j
, j1
, k
, k0
,
249 k_check
, leftright
, m2
, m5
, s2
, s5
, spec_case
, try_quick
;
250 union double_union d
, d2
, eps
;
252 #ifndef Sudden_Underflow
256 _Jv_Bigint
*b
, *b1
, *delta
, *mlo
, *mhi
, *S
;
264 ptr
->_result
->_k
= ptr
->_result_k
;
265 ptr
->_result
->_maxwds
= 1 << ptr
->_result_k
;
266 Bfree (ptr
, ptr
->_result
);
270 if (word0 (d
) & Sign_bit
)
272 /* set sign for everything, including 0's and NaNs */
274 word0 (d
) &= ~Sign_bit
; /* clear sign bit */
279 #if defined(IEEE_Arith) + defined(VAX)
281 if ((word0 (d
) & Exp_mask
) == Exp_mask
)
283 if (word0 (d
) == 0x8000)
286 /* Infinity or NaN */
290 !word1 (d
) && !(word0 (d
) & 0xfffff) ? "Infinity" :
303 d
.d
+= 0; /* normalize */
314 b
= d2b (ptr
, d
.d
, &be
, &bbits
);
315 #ifdef Sudden_Underflow
316 i
= (int) (word0 (d
) >> Exp_shift1
& (Exp_mask
>> Exp_shift1
));
318 if ((i
= (int) (word0 (d
) >> Exp_shift1
& (Exp_mask
>> Exp_shift1
))))
322 word0 (d2
) &= Frac_mask1
;
323 word0 (d2
) |= Exp_11
;
325 if (j
= 11 - hi0bits (word0 (d2
) & Frac_mask
))
329 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
330 * log10(x) = log(x) / log(10)
331 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
332 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
334 * This suggests computing an approximation k to log10(d) by
336 * k = (i - Bias)*0.301029995663981
337 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
339 * We want k to be too large rather than too small.
340 * The error in the first-order Taylor series approximation
341 * is in our favor, so we just round up the constant enough
342 * to compensate for any error in the multiplication of
343 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
344 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
345 * adding 1e-13 to the constant term more than suffices.
346 * Hence we adjust the constant term to 0.1760912590558.
347 * (We could get a more accurate k by invoking log10,
348 * but this is probably not worthwhile.)
356 #ifndef Sudden_Underflow
361 /* d is denormalized */
363 i
= bbits
+ be
+ (Bias
+ (P
- 1) - 1);
364 x
= i
> 32 ? word0 (d
) << (64 - i
) | word1 (d
) >> (i
- 32)
365 : word1 (d
) << (32 - i
);
367 word0 (d2
) -= 31 * Exp_msk1
; /* adjust exponent */
368 i
-= (Bias
+ (P
- 1) - 1) + 1;
372 ds
= (d2
.d
- 1.5) * 0.289529654602168 + 0.1760912590558 + i
* 0.301029995663981;
374 if (ds
< 0. && ds
!= k
)
375 k
--; /* want k = floor(ds) */
377 if (k
>= 0 && k
<= Ten_pmax
)
406 if (mode
< 0 || mode
> 9)
429 ilim
= ilim1
= i
= ndigits
;
441 j
= sizeof (unsigned long);
442 for (ptr
->_result_k
= 0; (int) (sizeof (_Jv_Bigint
) - sizeof (unsigned long)) + j
<= i
;
445 ptr
->_result
= Balloc (ptr
, ptr
->_result_k
);
446 s
= s0
= (char *) ptr
->_result
;
448 if (ilim
>= 0 && ilim
<= Quick_max
&& try_quick
)
450 /* Try to get by with floating-point arithmetic. */
456 ieps
= 2; /* conservative */
463 /* prevent overflows */
465 d
.d
/= bigtens
[n_bigtens
- 1];
468 for (; j
; j
>>= 1, i
++)
478 d
.d
*= tens
[j1
& 0xf];
479 for (j
= j1
>> 4; j
; j
>>= 1, i
++)
486 if (k_check
&& d
.d
< 1. && ilim
> 0)
495 eps
.d
= ieps
* d
.d
+ 7.;
496 word0 (eps
) -= (P
- 1) * Exp_msk1
;
510 /* Use Steele & White method of only
511 * generating digits needed.
513 eps
.d
= 0.5 / tens
[ilim
- 1] - eps
.d
;
518 *s
++ = '0' + (int) L
;
521 if (1. - d
.d
< eps
.d
)
532 /* Generate ilim digits, then fix them up. */
533 eps
.d
*= tens
[ilim
- 1];
534 for (i
= 1;; i
++, d
.d
*= 10.)
538 *s
++ = '0' + (int) L
;
541 if (d
.d
> 0.5 + eps
.d
)
543 else if (d
.d
< 0.5 - eps
.d
)
562 /* Do we have a "small" integer? */
564 if (be
>= 0 && k
<= Int_max
)
568 if (ndigits
< 0 && ilim
<= 0)
571 if (ilim
< 0 || d
.d
<= 5 * ds
)
579 #ifdef Check_FLT_ROUNDS
580 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
587 *s
++ = '0' + (int) L
;
591 if (d
.d
> ds
|| (d
.d
== ds
&& L
& 1))
619 #ifndef Sudden_Underflow
620 denorm
? be
+ (Bias
+ (P
- 1) - 1 + 1) :
623 1 + 4 * P
- 3 - bbits
+ ((bbits
+ be
- 1) & 3);
649 if (m2
> 0 && s2
> 0)
651 i
= m2
< s2
? m2
: s2
;
662 mhi
= pow5mult (ptr
, mhi
, m5
);
663 b1
= mult (ptr
, mhi
, b
);
668 b
= pow5mult (ptr
, b
, j
);
671 b
= pow5mult (ptr
, b
, b5
);
675 S
= pow5mult (ptr
, S
, s5
);
677 /* Check for special case that d is a normalized power of 2. */
681 if (!word1 (d
) && !(word0 (d
) & Bndry_mask
)
682 #ifndef Sudden_Underflow
683 && word0(d
) & Exp_mask
687 /* The special case */
696 /* Arrange for convenient computation of quotients:
697 * shift left if necessary so divisor has 4 leading 0 bits.
699 * Perhaps we should just compute leading 28 bits of S once
700 * and for all and pass them and a shift to quorem, so it
701 * can do shifts and ors to compute the numerator for q.
705 if ((i
= ((s5
? 32 - hi0bits (S
->_x
[S
->_wds
- 1]) : 1) + s2
) & 0x1f))
708 if ((i
= ((s5
? 32 - hi0bits (S
->_x
[S
->_wds
- 1]) : 1) + s2
) & 0xf))
726 b
= lshift (ptr
, b
, b2
);
728 S
= lshift (ptr
, S
, s2
);
734 b
= multadd (ptr
, b
, 10, 0); /* we botched the k estimate */
736 mhi
= multadd (ptr
, mhi
, 10, 0);
740 if (ilim
<= 0 && mode
> 2)
742 if (ilim
< 0 || cmp (b
, S
= multadd (ptr
, S
, 5, 0)) <= 0)
744 /* no digits, fcvt style */
757 mhi
= lshift (ptr
, mhi
, m2
);
759 /* Single precision case, */
761 mhi
= lshift (ptr
, mhi
, 29);
763 /* Compute mlo -- check for special case
764 * that d is a normalized power of 2.
770 mhi
= Balloc (ptr
, mhi
->_k
);
772 mhi
= lshift (ptr
, mhi
, Log2P
);
777 dig
= quorem (b
, S
) + '0';
778 /* Do we yet have the shortest decimal string
779 * that will round to d?
782 delta
= diff (ptr
, S
, mhi
);
783 j1
= delta
->_sign
? 1 : cmp (b
, delta
);
786 if (j1
== 0 && !mode
&& !(word1 (d
) & 1))
796 if (j
< 0 || (j
== 0 && !mode
804 b
= lshift (ptr
, b
, 1);
806 if ((j1
> 0 || (j1
== 0 && dig
& 1))
816 { /* possible if i == 1 */
827 b
= multadd (ptr
, b
, 10, 0);
829 mlo
= mhi
= multadd (ptr
, mhi
, 10, 0);
832 mlo
= multadd (ptr
, mlo
, 10, 0);
833 mhi
= multadd (ptr
, mhi
, 10, 0);
840 *s
++ = dig
= quorem (b
, S
) + '0';
843 b
= multadd (ptr
, b
, 10, 0);
846 /* Round off last digit */
848 b
= lshift (ptr
, b
, 1);
850 if (j
> 0 || (j
== 0 && dig
& 1))
871 if (mlo
&& mlo
!= mhi
)
887 (_d
, mode
, ndigits
, decpt
, sign
, rve
, buf
, float_type
),
897 struct _Jv_reent reent
;
899 memset (&reent
, 0, sizeof reent
);
901 p
= _dtoa_r (&reent
, _d
, mode
, ndigits
, decpt
, sign
, rve
, float_type
);