]>
git.saurik.com Git - apple/libc.git/blob - gdtoa/FreeBSD/gdtoa-gdtoa.c
1 /****************************************************************
3 The author of this software is David M. Gay.
5 Copyright (C) 1998, 1999 by Lucent Technologies
8 Permission to use, copy, modify, and distribute this software and
9 its documentation for any purpose and without fee is hereby
10 granted, provided that the above copyright notice appear in all
11 copies and that both that the copyright notice and this
12 permission notice and warranty disclaimer appear in supporting
13 documentation, and that the name of Lucent or any of its entities
14 not be used in advertising or publicity pertaining to
15 distribution of the software without specific, written prior
18 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
19 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
20 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
21 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
22 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
23 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
24 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
27 ****************************************************************/
29 /* Please send bug reports to David M. Gay (dmg at acm dot org,
30 * with " at " changed at "@" and " dot " changed to "."). */
36 bitstob(bits
, nbits
, bbits
) ULong
*bits
; int nbits
; int *bbits
;
38 bitstob(ULong
*bits
, int nbits
, int *bbits
)
56 be
= bits
+ ((nbits
- 1) >> kshift
);
59 *x
++ = *bits
& ALL_ON
;
61 *x
++ = (*bits
>> 16) & ALL_ON
;
63 } while(++bits
<= be
);
72 *bbits
= i
*ULbits
+ 32 - hi0bits(b
->x
[i
]);
77 /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
79 * Inspired by "How to Print Floating-Point Numbers Accurately" by
80 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
83 * 1. Rather than iterating, we use a simple numeric overestimate
84 * to determine k = floor(log10(d)). We scale relevant
85 * quantities using O(log2(k)) rather than O(k) multiplications.
86 * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
87 * try to generate digits strictly left to right. Instead, we
88 * compute with fewer bits and propagate the carry if necessary
89 * when rounding the final digit up. This is often faster.
90 * 3. Under the assumption that input will be rounded nearest,
91 * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
92 * That is, we allow equality in stopping tests when the
93 * round-nearest rule will give the same floating-point value
94 * as would satisfaction of the stopping test with strict
96 * 4. We remove common factors of powers of 2 from relevant
98 * 5. When converting floating-point integers less than 1e16,
99 * we use floating-point arithmetic rather than resorting
100 * to multiple-precision integers.
101 * 6. When asked to produce fewer than 15 digits, we first try
102 * to get by with floating-point arithmetic; we resort to
103 * multiple-precision integer arithmetic only if we cannot
104 * guarantee that the floating-point calculation has given
105 * the correctly rounded result. For k requested digits and
106 * "uniformly" distributed input, the probability is
107 * something like 10^(k-15) that we must resort to the Long
114 (fpi
, be
, bits
, kindp
, mode
, ndigits
, decpt
, rve
)
115 FPI
*fpi
; int be
; ULong
*bits
;
116 int *kindp
, mode
, ndigits
, *decpt
; char **rve
;
118 (FPI
*fpi
, int be
, ULong
*bits
, int *kindp
, int mode
, int ndigits
, int *decpt
, char **rve
)
121 /* Arguments ndigits and decpt are similar to the second and third
122 arguments of ecvt and fcvt; trailing zeros are suppressed from
123 the returned string. If not null, *rve is set to point
124 to the end of the return value. If d is +-Infinity or NaN,
125 then *decpt is set to 9999.
128 0 ==> shortest string that yields d when read in
129 and rounded to nearest.
130 1 ==> like 0, but with Steele & White stopping rule;
131 e.g. with IEEE P754 arithmetic , mode 0 gives
132 1e23 whereas mode 1 gives 9.999999999999999e22.
133 2 ==> max(1,ndigits) significant digits. This gives a
134 return value similar to that of ecvt, except
135 that trailing zeros are suppressed.
136 3 ==> through ndigits past the decimal point. This
137 gives a return value similar to that from fcvt,
138 except that trailing zeros are suppressed, and
139 ndigits can be negative.
140 4-9 should give the same return values as 2-3, i.e.,
141 4 <= mode <= 9 ==> same return as mode
142 2 + (mode & 1). These modes are mainly for
143 debugging; often they run slower but sometimes
144 faster than modes 2-3.
145 4,5,8,9 ==> left-to-right digit generation.
146 6-9 ==> don't try fast floating-point estimate
149 Values of mode other than 0-9 are treated as mode 0.
151 Sufficient space is allocated to the return value
152 to hold the suppressed trailing zeros.
155 int bbits
, b2
, b5
, be0
, dig
, i
, ieps
, ilim
, ilim0
, ilim1
, inex
;
156 int j
, j1
, k
, k0
, k_check
, kind
, leftright
, m2
, m5
, nbits
;
157 int rdir
, s2
, s5
, spec_case
, try_quick
;
159 Bigint
*b
, *b1
, *delta
, *mlo
, *mhi
, *mhi1
, *S
;
160 double d
, d2
, ds
, eps
;
163 #ifndef MULTIPLE_THREADS
165 freedtoa(dtoa_result
);
170 kind
= *kindp
&= ~STRTOG_Inexact
;
171 switch(kind
& STRTOG_Retmask
) {
175 case STRTOG_Denormal
:
177 case STRTOG_Infinite
:
179 return nrv_alloc("Infinity", rve
, 8);
182 return nrv_alloc("NaN", rve
, 3);
186 b
= bitstob(bits
, nbits
= fpi
->nbits
, &bbits
);
188 if ( (i
= trailz(b
)) !=0) {
197 return nrv_alloc("0", rve
, 1);
200 dval(d
) = b2d(b
, &i
);
202 word0(d
) &= Frac_mask1
;
205 if ( (j
= 11 - hi0bits(word0(d
) & Frac_mask
)) !=0)
209 /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
210 * log10(x) = log(x) / log(10)
211 * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
212 * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
214 * This suggests computing an approximation k to log10(d) by
216 * k = (i - Bias)*0.301029995663981
217 * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
219 * We want k to be too large rather than too small.
220 * The error in the first-order Taylor series approximation
221 * is in our favor, so we just round up the constant enough
222 * to compensate for any error in the multiplication of
223 * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
224 * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
225 * adding 1e-13 to the constant term more than suffices.
226 * Hence we adjust the constant term to 0.1760912590558.
227 * (We could get a more accurate k by invoking log10,
228 * but this is probably not worthwhile.)
234 ds
= (dval(d
)-1.5)*0.289529654602168 + 0.1760912590558 + i
*0.301029995663981;
236 /* correct assumption about exponent range */
243 if (ds
< 0. && ds
!= k
)
244 k
--; /* want k = floor(ds) */
248 if ( (j1
= j
& 3) !=0)
250 word0(d
) += j
<< Exp_shift
- 2 & Exp_mask
;
252 word0(d
) += (be
+ bbits
- 1) << Exp_shift
;
254 if (k
>= 0 && k
<= Ten_pmax
) {
255 if (dval(d
) < tens
[k
])
278 if (mode
< 0 || mode
> 9)
290 i
= (int)(nbits
* .30103) + 3;
299 ilim
= ilim1
= i
= ndigits
;
311 s
= s0
= rv_alloc(i
);
313 if ( (rdir
= fpi
->rounding
- 1) !=0) {
316 if (kind
& STRTOG_Neg
)
320 /* Now rdir = 0 ==> round near, 1 ==> round up, 2 ==> round down. */
322 if (ilim
>= 0 && ilim
<= Quick_max
&& try_quick
&& !rdir
323 #ifndef IMPRECISE_INEXACT
328 /* Try to get by with floating-point arithmetic. */
333 if ( (j
= 11 - hi0bits(word0(d
) & Frac_mask
)) !=0)
338 ieps
= 2; /* conservative */
343 /* prevent overflows */
345 dval(d
) /= bigtens
[n_bigtens
-1];
348 for(; j
; j
>>= 1, i
++)
356 if ( (j1
= -k
) !=0) {
357 dval(d
) *= tens
[j1
& 0xf];
358 for(j
= j1
>> 4; j
; j
>>= 1, i
++)
361 dval(d
) *= bigtens
[i
];
365 if (k_check
&& dval(d
) < 1. && ilim
> 0) {
373 dval(eps
) = ieps
*dval(d
) + 7.;
374 word0(eps
) -= (P
-1)*Exp_msk1
;
378 if (dval(d
) > dval(eps
))
380 if (dval(d
) < -dval(eps
))
386 /* Use Steele & White method of only
387 * generating digits needed.
389 dval(eps
) = ds
*0.5/tens
[ilim
-1] - dval(eps
);
391 L
= (Long
)(dval(d
)/ds
);
394 if (dval(d
) < dval(eps
)) {
396 inex
= STRTOG_Inexlo
;
399 if (ds
- dval(d
) < dval(eps
))
409 /* Generate ilim digits, then fix them up. */
410 dval(eps
) *= tens
[ilim
-1];
411 for(i
= 1;; i
++, dval(d
) *= 10.) {
412 if ( (L
= (Long
)(dval(d
)/ds
)) !=0)
417 if (dval(d
) > ds
+ dval(eps
))
419 else if (dval(d
) < ds
- dval(eps
)) {
421 inex
= STRTOG_Inexlo
;
422 goto clear_trailing0
;
437 /* Do we have a "small" integer? */
439 if (be
>= 0 && k
<= Int_max
) {
442 if (ndigits
< 0 && ilim
<= 0) {
444 if (ilim
< 0 || dval(d
) <= 5*ds
)
448 for(i
= 1;; i
++, dval(d
) *= 10.) {
451 #ifdef Check_FLT_ROUNDS
452 /* If FLT_ROUNDS == 2, L will usually be high by 1 */
465 inex
= STRTOG_Inexlo
;
469 if (dval(d
) > ds
|| dval(d
) == ds
&& L
& 1) {
471 inex
= STRTOG_Inexhi
;
481 inex
= STRTOG_Inexlo
;
498 if (be
- i
++ < fpi
->emin
)
500 i
= be
- fpi
->emin
+ 1;
511 if ((i
= ilim
) < 0) {
520 if (m2
> 0 && s2
> 0) {
521 i
= m2
< s2
? m2
: s2
;
529 mhi
= pow5mult(mhi
, m5
);
534 if ( (j
= b5
- m5
) !=0)
544 /* Check for special case that d is a normalized power of 2. */
548 if (bbits
== 1 && be0
> fpi
->emin
+ 1) {
549 /* The special case */
556 /* Arrange for convenient computation of quotients:
557 * shift left if necessary so divisor has 4 leading 0 bits.
559 * Perhaps we should just compute leading 28 bits of S once
560 * and for all and pass them and a shift to quorem, so it
561 * can do shifts and ors to compute the numerator for q.
564 if ( (i
= ((s5
? 32 - hi0bits(S
->x
[S
->wds
-1]) : 1) + s2
) & 0x1f) !=0)
567 if ( (i
= ((s5
? 32 - hi0bits(S
->x
[S
->wds
-1]) : 1) + s2
) & 0xf) !=0)
589 b
= multadd(b
, 10, 0); /* we botched the k estimate */
591 mhi
= multadd(mhi
, 10, 0);
595 if (ilim
<= 0 && mode
> 2) {
596 if (ilim
< 0 || cmp(b
,S
= multadd(S
,5,0)) <= 0) {
597 /* no digits, fcvt style */
600 inex
= STRTOG_Inexlo
;
604 inex
= STRTOG_Inexhi
;
611 mhi
= lshift(mhi
, m2
);
613 /* Compute mlo -- check for special case
614 * that d is a normalized power of 2.
619 mhi
= Balloc(mhi
->k
);
621 mhi
= lshift(mhi
, 1);
625 dig
= quorem(b
,S
) + '0';
626 /* Do we yet have the shortest decimal string
627 * that will round to d?
630 delta
= diff(S
, mhi
);
631 j1
= delta
->sign
? 1 : cmp(b
, delta
);
634 if (j1
== 0 && !mode
&& !(bits
[0] & 1) && !rdir
) {
638 if (b
->wds
> 1 || b
->x
[0])
639 inex
= STRTOG_Inexlo
;
643 inex
= STRTOG_Inexhi
;
649 if (j
< 0 || j
== 0 && !mode
654 if (rdir
&& (b
->wds
> 1 || b
->x
[0])) {
656 inex
= STRTOG_Inexlo
;
659 while (cmp(S
,mhi
) > 0) {
661 mhi1
= multadd(mhi
, 10, 0);
665 b
= multadd(b
, 10, 0);
666 dig
= quorem(b
,S
) + '0';
670 inex
= STRTOG_Inexhi
;
676 if ((j1
> 0 || j1
== 0 && dig
& 1)
679 inex
= STRTOG_Inexhi
;
681 if (b
->wds
> 1 || b
->x
[0])
682 inex
= STRTOG_Inexlo
;
687 if (j1
> 0 && rdir
!= 2) {
688 if (dig
== '9') { /* possible if i == 1 */
691 inex
= STRTOG_Inexhi
;
694 inex
= STRTOG_Inexhi
;
701 b
= multadd(b
, 10, 0);
703 mlo
= mhi
= multadd(mhi
, 10, 0);
705 mlo
= multadd(mlo
, 10, 0);
706 mhi
= multadd(mhi
, 10, 0);
712 *s
++ = dig
= quorem(b
,S
) + '0';
715 b
= multadd(b
, 10, 0);
718 /* Round off last digit */
721 if (rdir
== 2 || b
->wds
<= 1 && !b
->x
[0])
727 if (j
> 0 || j
== 0 && dig
& 1) {
729 inex
= STRTOG_Inexhi
;
740 if (b
->wds
> 1 || b
->x
[0])
741 inex
= STRTOG_Inexlo
;
748 if (mlo
&& mlo
!= mhi
)