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
);