1 #include "fconv.h"
2
3 /* strtod for IEEE-, VAX-, and IBM-arithmetic machines (dmg).
4 *
5 * This strtod returns a nearest machine number to the input decimal
6 * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
7 * broken by the IEEE round-even rule. Otherwise ties are broken by
8 * biased rounding (add half and chop).
9 *
10 * Inspired loosely by William D. Clinger's paper "How to Read Floating
11 * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
12 *
13 * Modifications:
14 *
15 * 1. We only require IEEE, IBM, or VAX double-precision
16 * arithmetic (not IEEE double-extended).
17 * 2. We get by with floating-point arithmetic in a case that
18 * Clinger missed -- when we're computing d * 10^n
19 * for a small integer d and the integer n is not too
20 * much larger than 22 (the maximum integer k for which
21 * we can represent 10^k exactly), we may be able to
22 * compute (d*10^k) * 10^(e-k) with just one roundoff.
23 * 3. Rather than a bit-at-a-time adjustment of the binary
24 * result in the hard case, we use floating-point
25 * arithmetic to determine the adjustment to within
26 * one bit; only in really hard cases do we need to
27 * compute a second residual.
28 * 4. Because of 3., we don't need a large table of powers of 10
29 * for ten-to-e (just some small tables, e.g. of 10^k
30 * for 0 <= k <= 22).
31 */
32
33 #ifdef RND_PRODQUOT
34 #define rounded_product(a,b) a = rnd_prod(a, b)
35 #define rounded_quotient(a,b) a = rnd_quot(a, b)
36 extern double rnd_prod(double, double), rnd_quot(double, double);
37 #else
38 #define rounded_product(a,b) a *= b
39 #define rounded_quotient(a,b) a /= b
40 #endif
41
42 static double
ulp(double xarg)43 ulp(double xarg)
44 {
45 register long L;
46 Dul a;
47 Dul x;
48
49 x.d = xarg;
50 L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
51 #ifndef Sudden_Underflow
52 if (L > 0) {
53 #endif
54 #ifdef IBM
55 L |= Exp_msk1 >> 4;
56 #endif
57 word0(a) = L;
58 word1(a) = 0;
59 #ifndef Sudden_Underflow
60 }
61 else {
62 L = -L >> Exp_shift;
63 if (L < Exp_shift) {
64 word0(a) = 0x80000 >> L;
65 word1(a) = 0;
66 }
67 else {
68 word0(a) = 0;
69 L -= Exp_shift;
70 word1(a) = L >= 31 ? 1 : 1 << 31 - L;
71 }
72 }
73 #endif
74 return a.d;
75 }
76
77 static Bigint *
s2b(CONST char * s,int nd0,int nd,unsigned long y9)78 s2b(CONST char *s, int nd0, int nd, unsigned long y9)
79 {
80 Bigint *b;
81 int i, k;
82 long x, y;
83
84 x = (nd + 8) / 9;
85 for(k = 0, y = 1; x > y; y <<= 1, k++) ;
86 #ifdef Pack_32
87 b = Balloc(k);
88 b->x[0] = y9;
89 b->wds = 1;
90 #else
91 b = Balloc(k+1);
92 b->x[0] = y9 & 0xffff;
93 b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
94 #endif
95
96 i = 9;
97 if (9 < nd0) {
98 s += 9;
99 do b = multadd(b, 10, *s++ - '0');
100 while(++i < nd0);
101 s++;
102 }
103 else
104 s += 10;
105 for(; i < nd; i++)
106 b = multadd(b, 10, *s++ - '0');
107 return b;
108 }
109
110 static double
b2d(Bigint * a,int * e)111 b2d(Bigint *a, int *e)
112 {
113 unsigned long *xa, *xa0, w, y, z;
114 int k;
115 Dul d;
116 #ifdef VAX
117 unsigned long d0, d1;
118 #else
119 #define d0 word0(d)
120 #define d1 word1(d)
121 #endif
122
123 xa0 = a->x;
124 xa = xa0 + a->wds;
125 y = *--xa;
126 #ifdef DEBUG
127 if (!y) Bug("zero y in b2d");
128 #endif
129 k = hi0bits(y);
130 *e = 32 - k;
131 #ifdef Pack_32
132 if (k < Ebits) {
133 d0 = Exp_1 | y >> Ebits - k;
134 w = xa > xa0 ? *--xa : 0;
135 d1 = y << (32-Ebits) + k | w >> Ebits - k;
136 goto ret_d;
137 }
138 z = xa > xa0 ? *--xa : 0;
139 if (k -= Ebits) {
140 d0 = Exp_1 | y << k | z >> 32 - k;
141 y = xa > xa0 ? *--xa : 0;
142 d1 = z << k | y >> 32 - k;
143 }
144 else {
145 d0 = Exp_1 | y;
146 d1 = z;
147 }
148 #else
149 if (k < Ebits + 16) {
150 z = xa > xa0 ? *--xa : 0;
151 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
152 w = xa > xa0 ? *--xa : 0;
153 y = xa > xa0 ? *--xa : 0;
154 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
155 goto ret_d;
156 }
157 z = xa > xa0 ? *--xa : 0;
158 w = xa > xa0 ? *--xa : 0;
159 k -= Ebits + 16;
160 d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
161 y = xa > xa0 ? *--xa : 0;
162 d1 = w << k + 16 | y << k;
163 #endif
164 ret_d:
165 #ifdef VAX
166 word0(d) = d0 >> 16 | d0 << 16;
167 word1(d) = d1 >> 16 | d1 << 16;
168 #else
169 #undef d0
170 #undef d1
171 #endif
172 return d.d;
173 }
174
175 static double
ratio(Bigint * a,Bigint * b)176 ratio(Bigint *a, Bigint *b)
177 {
178 Dul da, db;
179 int k, ka, kb;
180
181 da.d = b2d(a, &ka);
182 db.d = b2d(b, &kb);
183 #ifdef Pack_32
184 k = ka - kb + 32*(a->wds - b->wds);
185 #else
186 k = ka - kb + 16*(a->wds - b->wds);
187 #endif
188 #ifdef IBM
189 if (k > 0) {
190 word0(da) += (k >> 2)*Exp_msk1;
191 if (k &= 3)
192 da *= 1 << k;
193 }
194 else {
195 k = -k;
196 word0(db) += (k >> 2)*Exp_msk1;
197 if (k &= 3)
198 db *= 1 << k;
199 }
200 #else
201 if (k > 0)
202 word0(da) += k*Exp_msk1;
203 else {
204 k = -k;
205 word0(db) += k*Exp_msk1;
206 }
207 #endif
208 return da.d / db.d;
209 }
210
211 double
strtod(CONST char * s00,char ** se)212 strtod(CONST char *s00, char **se)
213 {
214 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
215 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
216 CONST char *s, *s0, *s1;
217 double aadj, aadj1, adj;
218 Dul rv, rv0;
219 long L;
220 unsigned long y, z;
221 Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
222
223 bb = bd = bs = delta = 0;
224 sign = nz0 = nz = 0;
225 rv.d = 0.;
226 for(s = s00;;s++) switch(*s) {
227 case '-':
228 sign = 1;
229 /* no break */
230 case '+':
231 if (*++s)
232 goto break2;
233 /* no break */
234 case 0:
235 s = s00;
236 goto ret;
237 case '\t':
238 case '\n':
239 case '\v':
240 case '\f':
241 case '\r':
242 case ' ':
243 continue;
244 default:
245 goto break2;
246 }
247 break2:
248 if (*s == '0') {
249 nz0 = 1;
250 while(*++s == '0') ;
251 if (!*s)
252 goto ret;
253 }
254 s0 = s;
255 y = z = 0;
256 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
257 if (nd < 9)
258 y = 10*y + c - '0';
259 else if (nd < 16)
260 z = 10*z + c - '0';
261 nd0 = nd;
262 if (c == '.') {
263 c = *++s;
264 if (!nd) {
265 for(; c == '0'; c = *++s)
266 nz++;
267 if (c > '0' && c <= '9') {
268 s0 = s;
269 nf += nz;
270 nz = 0;
271 goto have_dig;
272 }
273 goto dig_done;
274 }
275 for(; c >= '0' && c <= '9'; c = *++s) {
276 have_dig:
277 nz++;
278 if (c -= '0') {
279 nf += nz;
280 for(i = 1; i < nz; i++)
281 if (nd++ < 9)
282 y *= 10;
283 else if (nd <= DBL_DIG + 1)
284 z *= 10;
285 if (nd++ < 9)
286 y = 10*y + c;
287 else if (nd <= DBL_DIG + 1)
288 z = 10*z + c;
289 nz = 0;
290 }
291 }
292 }
293 dig_done:
294 e = 0;
295 if (c == 'e' || c == 'E') {
296 if (!nd && !nz && !nz0) {
297 s = s00;
298 goto ret;
299 }
300 s00 = s;
301 esign = 0;
302 switch(c = *++s) {
303 case '-':
304 esign = 1;
305 case '+':
306 c = *++s;
307 }
308 if (c >= '0' && c <= '9') {
309 while(c == '0')
310 c = *++s;
311 if (c > '0' && c <= '9') {
312 e = c - '0';
313 s1 = s;
314 while((c = *++s) >= '0' && c <= '9')
315 e = 10*e + c - '0';
316 if (s - s1 > 8)
317 /* Avoid confusion from exponents
318 * so large that e might overflow.
319 */
320 e = 9999999;
321 if (esign)
322 e = -e;
323 }
324 else
325 e = 0;
326 }
327 else
328 s = s00;
329 }
330 if (!nd) {
331 if (!nz && !nz0)
332 s = s00;
333 goto ret;
334 }
335 e1 = e -= nf;
336
337 /* Now we have nd0 digits, starting at s0, followed by a
338 * decimal point, followed by nd-nd0 digits. The number we're
339 * after is the integer represented by those digits times
340 * 10**e */
341
342 if (!nd0)
343 nd0 = nd;
344 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
345 rv.d = y;
346 if (k > 9)
347 rv.d = tens[k - 9] * rv.d + z;
348 bd0 = 0;
349 if (nd <= DBL_DIG
350 #ifndef RND_PRODQUOT
351 && FLT_ROUNDS == 1
352 #endif
353 ) {
354 if (!e)
355 goto ret;
356 if (e > 0) {
357 if (e <= Ten_pmax) {
358 #ifdef VAX
359 goto vax_ovfl_check;
360 #else
361 /* rv = */ rounded_product(rv.d, tens[e]);
362 goto ret;
363 #endif
364 }
365 i = DBL_DIG - nd;
366 if (e <= Ten_pmax + i) {
367 /* A fancier test would sometimes let us do
368 * this for larger i values.
369 */
370 e -= i;
371 rv.d *= tens[i];
372 #ifdef VAX
373 /* VAX exponent range is so narrow we must
374 * worry about overflow here...
375 */
376 vax_ovfl_check:
377 word0(rv) -= P*Exp_msk1;
378 /* rv = */ rounded_product(rv.d, tens[e]);
379 if ((word0(rv) & Exp_mask)
380 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
381 goto ovfl;
382 word0(rv) += P*Exp_msk1;
383 #else
384 /* rv = */ rounded_product(rv.d, tens[e]);
385 #endif
386 goto ret;
387 }
388 }
389 else if (e >= -Ten_pmax) {
390 /* rv = */ rounded_quotient(rv.d, tens[-e]);
391 goto ret;
392 }
393 }
394 e1 += nd - k;
395
396 /* Get starting approximation = rv * 10**e1 */
397
398 if (e1 > 0) {
399 if (nd0 + e1 - 1 > DBL_MAX_10_EXP)
400 goto ovfl;
401 if (i = e1 & 15)
402 rv.d *= tens[i];
403 if (e1 &= ~15) {
404 if (e1 > DBL_MAX_10_EXP) {
405 ovfl:
406 errno = ERANGE;
407 rv.d = HUGE_VAL;
408 if (bd0)
409 goto retfree;
410 goto ret;
411 }
412 if (e1 >>= 4) {
413 for(j = 0; e1 > 1; j++, e1 >>= 1)
414 if (e1 & 1)
415 rv.d *= bigtens[j];
416 /* The last multiplication could overflow. */
417 word0(rv) -= P*Exp_msk1;
418 rv.d *= bigtens[j];
419 if ((z = word0(rv) & Exp_mask)
420 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
421 goto ovfl;
422 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
423 /* set to largest number */
424 /* (Can't trust DBL_MAX) */
425 word0(rv) = Big0;
426 word1(rv) = Big1;
427 }
428 else
429 word0(rv) += P*Exp_msk1;
430 }
431
432 }
433 }
434 else if (e1 < 0) {
435 e1 = -e1;
436 if (i = e1 & 15)
437 rv.d /= tens[i];
438 if (e1 &= ~15) {
439 e1 >>= 4;
440 if (e1 >= 1 << n_bigtens)
441 goto undfl;
442 for(j = 0; e1 > 1; j++, e1 >>= 1)
443 if (e1 & 1)
444 rv.d *= tinytens[j];
445 /* The last multiplication could underflow. */
446 rv0.d = rv.d;
447 rv.d *= tinytens[j];
448 if (rv.d == 0) {
449 rv.d = 2.*rv0.d;
450 rv.d *= tinytens[j];
451 if (rv.d == 0) {
452 undfl:
453 rv.d = 0.;
454 errno = ERANGE;
455 if (bd0)
456 goto retfree;
457 goto ret;
458 }
459 word0(rv) = Tiny0;
460 word1(rv) = Tiny1;
461 /* The refinement below will clean
462 * this approximation up.
463 */
464 }
465 }
466 }
467
468 /* Now the hard part -- adjusting rv to the correct value.*/
469
470 /* Put digits into bd: true value = bd * 10^e */
471
472 bd0 = s2b(s0, nd0, nd, y);
473
474 for(;;) {
475 bd = Balloc(bd0->k);
476 Bcopy(bd, bd0);
477 bb = d2b(rv.d, &bbe, &bbbits); /* rv = bb * 2^bbe */
478 bs = i2b(1);
479
480 if (e >= 0) {
481 bb2 = bb5 = 0;
482 bd2 = bd5 = e;
483 }
484 else {
485 bb2 = bb5 = -e;
486 bd2 = bd5 = 0;
487 }
488 if (bbe >= 0)
489 bb2 += bbe;
490 else
491 bd2 -= bbe;
492 bs2 = bb2;
493 #ifdef Sudden_Underflow
494 #ifdef IBM
495 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
496 #else
497 j = P + 1 - bbbits;
498 #endif
499 #else
500 i = bbe + bbbits - 1; /* logb(rv) */
501 if (i < Emin) /* denormal */
502 j = bbe + (P-Emin);
503 else
504 j = P + 1 - bbbits;
505 #endif
506 bb2 += j;
507 bd2 += j;
508 i = bb2 < bd2 ? bb2 : bd2;
509 if (i > bs2)
510 i = bs2;
511 if (i > 0) {
512 bb2 -= i;
513 bd2 -= i;
514 bs2 -= i;
515 }
516 if (bb5 > 0) {
517 bs = pow5mult(bs, bb5);
518 bb1 = mult(bs, bb);
519 Bfree(bb);
520 bb = bb1;
521 }
522 if (bb2 > 0)
523 bb = lshift(bb, bb2);
524 if (bd5 > 0)
525 bd = pow5mult(bd, bd5);
526 if (bd2 > 0)
527 bd = lshift(bd, bd2);
528 if (bs2 > 0)
529 bs = lshift(bs, bs2);
530 delta = diff(bb, bd);
531 dsign = delta->sign;
532 delta->sign = 0;
533 i = cmp(delta, bs);
534 if (i < 0) {
535 /* Error is less than half an ulp -- check for
536 * special case of mantissa a power of two.
537 */
538 if (dsign || word1(rv) || word0(rv) & Bndry_mask)
539 break;
540 delta = lshift(delta,Log2P);
541 if (cmp(delta, bs) > 0)
542 goto drop_down;
543 break;
544 }
545 if (i == 0) {
546 /* exactly half-way between */
547 if (dsign) {
548 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
549 && word1(rv) == 0xffffffff) {
550 /*boundary case -- increment exponent*/
551 word0(rv) = (word0(rv) & Exp_mask)
552 + Exp_msk1
553 #ifdef IBM
554 | Exp_msk1 >> 4
555 #endif
556 ;
557 word1(rv) = 0;
558 break;
559 }
560 }
561 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
562 drop_down:
563 /* boundary case -- decrement exponent */
564 #ifdef Sudden_Underflow
565 L = word0(rv) & Exp_mask;
566 #ifdef IBM
567 if (L < Exp_msk1)
568 #else
569 if (L <= Exp_msk1)
570 #endif
571 goto undfl;
572 L -= Exp_msk1;
573 #else
574 L = (word0(rv) & Exp_mask) - Exp_msk1;
575 #endif
576 word0(rv) = L | Bndry_mask1;
577 word1(rv) = 0xffffffff;
578 #ifdef IBM
579 goto cont;
580 #else
581 break;
582 #endif
583 }
584 #ifndef ROUND_BIASED
585 if (!(word1(rv) & LSB))
586 break;
587 #endif
588 if (dsign)
589 rv.d += ulp(rv.d);
590 #ifndef ROUND_BIASED
591 else {
592 rv.d -= ulp(rv.d);
593 #ifndef Sudden_Underflow
594 if (rv.d == 0)
595 goto undfl;
596 #endif
597 }
598 #endif
599 break;
600 }
601 if ((aadj = ratio(delta, bs)) <= 2.) {
602 if (dsign)
603 aadj = aadj1 = 1.;
604 else if (word1(rv) || word0(rv) & Bndry_mask) {
605 #ifndef Sudden_Underflow
606 if (word1(rv) == Tiny1 && !word0(rv))
607 goto undfl;
608 #endif
609 aadj = 1.;
610 aadj1 = -1.;
611 }
612 else {
613 /* special case -- power of FLT_RADIX to be */
614 /* rounded down... */
615
616 if (aadj < 2./FLT_RADIX)
617 aadj = 1./FLT_RADIX;
618 else
619 aadj *= 0.5;
620 aadj1 = -aadj;
621 }
622 }
623 else {
624 aadj *= 0.5;
625 aadj1 = dsign ? aadj : -aadj;
626 #ifdef Check_FLT_ROUNDS
627 switch(FLT_ROUNDS) {
628 case 2: /* towards +infinity */
629 aadj1 -= 0.5;
630 break;
631 case 0: /* towards 0 */
632 case 3: /* towards -infinity */
633 aadj1 += 0.5;
634 }
635 #else
636 if (FLT_ROUNDS == 0)
637 aadj1 += 0.5;
638 #endif
639 }
640 y = word0(rv) & Exp_mask;
641
642 /* Check for overflow */
643
644 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
645 rv0.d = rv.d;
646 word0(rv) -= P*Exp_msk1;
647 adj = aadj1 * ulp(rv.d);
648 rv.d += adj;
649 if ((word0(rv) & Exp_mask) >=
650 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
651 if (word0(rv0) == Big0 && word1(rv0) == Big1)
652 goto ovfl;
653 word0(rv) = Big0;
654 word1(rv) = Big1;
655 goto cont;
656 }
657 else
658 word0(rv) += P*Exp_msk1;
659 }
660 else {
661 #ifdef Sudden_Underflow
662 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
663 rv0.d = rv.d;
664 word0(rv) += P*Exp_msk1;
665 adj = aadj1 * ulp(rv.d);
666 rv.d += adj;
667 #ifdef IBM
668 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
669 #else
670 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
671 #endif
672 {
673 if (word0(rv0) == Tiny0
674 && word1(rv0) == Tiny1)
675 goto undfl;
676 word0(rv) = Tiny0;
677 word1(rv) = Tiny1;
678 goto cont;
679 }
680 else
681 word0(rv) -= P*Exp_msk1;
682 }
683 else {
684 adj = aadj1 * ulp(rv.d);
685 rv.d += adj;
686 }
687 #else
688 /* Compute adj so that the IEEE rounding rules will
689 * correctly round rv + adj in some half-way cases.
690 * If rv * ulp(rv) is denormalized (i.e.,
691 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
692 * trouble from bits lost to denormalization;
693 * example: 1.2e-307 .
694 */
695 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
696 aadj1 = (double)(int)(aadj + 0.5);
697 if (!dsign)
698 aadj1 = -aadj1;
699 }
700 adj = aadj1 * ulp(rv.d);
701 rv.d += adj;
702 #endif
703 }
704 z = word0(rv) & Exp_mask;
705 if (y == z) {
706 /* Can we stop now? */
707 L = aadj;
708 aadj -= L;
709 /* The tolerances below are conservative. */
710 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
711 if (aadj < .4999999 || aadj > .5000001)
712 break;
713 }
714 else if (aadj < .4999999/FLT_RADIX)
715 break;
716 }
717 cont:
718 Bfree(bb);
719 Bfree(bd);
720 Bfree(bs);
721 Bfree(delta);
722 }
723 retfree:
724 Bfree(bb);
725 Bfree(bd);
726 Bfree(bs);
727 Bfree(bd0);
728 Bfree(delta);
729 ret:
730 if (se)
731 *se = (char *)s;
732 return sign ? -rv.d : rv.d;
733 }
734