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
, 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
;
900 memset (&reent
, 0, sizeof reent
);
902 p
= _dtoa_r (&reent
, _d
, mode
, ndigits
, decpt
, sign
, rve
, float_type
);