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 sign = nz0 = nz = 0;
223 rv.d = 0.;
224 for(s = s00;;s++) switch(*s) {
225 case '-':
226 sign = 1;
227 /* no break */
228 case '+':
229 if (*++s)
230 goto break2;
231 /* no break */
232 case 0:
233 s = s00;
234 goto ret;
235 case '\t':
236 case '\n':
237 case '\v':
238 case '\f':
239 case '\r':
240 case ' ':
241 continue;
242 default:
243 goto break2;
244 }
245 break2:
246 if (*s == '0') {
247 nz0 = 1;
248 while(*++s == '0') ;
249 if (!*s)
250 goto ret;
251 }
252 s0 = s;
253 y = z = 0;
254 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
255 if (nd < 9)
256 y = 10*y + c - '0';
257 else if (nd < 16)
258 z = 10*z + c - '0';
259 nd0 = nd;
260 if (c == '.') {
261 c = *++s;
262 if (!nd) {
263 for(; c == '0'; c = *++s)
264 nz++;
265 if (c > '0' && c <= '9') {
266 s0 = s;
267 nf += nz;
268 nz = 0;
269 goto have_dig;
270 }
271 goto dig_done;
272 }
273 for(; c >= '0' && c <= '9'; c = *++s) {
274 have_dig:
275 nz++;
276 if (c -= '0') {
277 nf += nz;
278 for(i = 1; i < nz; i++)
279 if (nd++ < 9)
280 y *= 10;
281 else if (nd <= DBL_DIG + 1)
282 z *= 10;
283 if (nd++ < 9)
284 y = 10*y + c;
285 else if (nd <= DBL_DIG + 1)
286 z = 10*z + c;
287 nz = 0;
288 }
289 }
290 }
291 dig_done:
292 e = 0;
293 if (c == 'e' || c == 'E') {
294 if (!nd && !nz && !nz0) {
295 s = s00;
296 goto ret;
297 }
298 s00 = s;
299 esign = 0;
300 switch(c = *++s) {
301 case '-':
302 esign = 1;
303 case '+':
304 c = *++s;
305 }
306 if (c >= '0' && c <= '9') {
307 while(c == '0')
308 c = *++s;
309 if (c > '0' && c <= '9') {
310 e = c - '0';
311 s1 = s;
312 while((c = *++s) >= '0' && c <= '9')
313 e = 10*e + c - '0';
314 if (s - s1 > 8)
315 /* Avoid confusion from exponents
316 * so large that e might overflow.
317 */
318 e = 9999999;
319 if (esign)
320 e = -e;
321 }
322 else
323 e = 0;
324 }
325 else
326 s = s00;
327 }
328 if (!nd) {
329 if (!nz && !nz0)
330 s = s00;
331 goto ret;
332 }
333 e1 = e -= nf;
334
335 /* Now we have nd0 digits, starting at s0, followed by a
336 * decimal point, followed by nd-nd0 digits. The number we're
337 * after is the integer represented by those digits times
338 * 10**e */
339
340 if (!nd0)
341 nd0 = nd;
342 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
343 rv.d = y;
344 if (k > 9)
345 rv.d = tens[k - 9] * rv.d + z;
346 bd0 = 0;
347 if (nd <= DBL_DIG
348 #ifndef RND_PRODQUOT
349 && FLT_ROUNDS == 1
350 #endif
351 ) {
352 if (!e)
353 goto ret;
354 if (e > 0) {
355 if (e <= Ten_pmax) {
356 #ifdef VAX
357 goto vax_ovfl_check;
358 #else
359 /* rv = */ rounded_product(rv.d, tens[e]);
360 goto ret;
361 #endif
362 }
363 i = DBL_DIG - nd;
364 if (e <= Ten_pmax + i) {
365 /* A fancier test would sometimes let us do
366 * this for larger i values.
367 */
368 e -= i;
369 rv.d *= tens[i];
370 #ifdef VAX
371 /* VAX exponent range is so narrow we must
372 * worry about overflow here...
373 */
374 vax_ovfl_check:
375 word0(rv) -= P*Exp_msk1;
376 /* rv = */ rounded_product(rv.d, tens[e]);
377 if ((word0(rv) & Exp_mask)
378 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
379 goto ovfl;
380 word0(rv) += P*Exp_msk1;
381 #else
382 /* rv = */ rounded_product(rv.d, tens[e]);
383 #endif
384 goto ret;
385 }
386 }
387 else if (e >= -Ten_pmax) {
388 /* rv = */ rounded_quotient(rv.d, tens[-e]);
389 goto ret;
390 }
391 }
392 e1 += nd - k;
393
394 /* Get starting approximation = rv * 10**e1 */
395
396 if (e1 > 0) {
397 if (nd0 + e1 - 1 > DBL_MAX_10_EXP)
398 goto ovfl;
399 if (i = e1 & 15)
400 rv.d *= tens[i];
401 if (e1 &= ~15) {
402 if (e1 > DBL_MAX_10_EXP) {
403 ovfl:
404 errno = ERANGE;
405 rv.d = HUGE_VAL;
406 if (bd0)
407 goto retfree;
408 goto ret;
409 }
410 if (e1 >>= 4) {
411 for(j = 0; e1 > 1; j++, e1 >>= 1)
412 if (e1 & 1)
413 rv.d *= bigtens[j];
414 /* The last multiplication could overflow. */
415 word0(rv) -= P*Exp_msk1;
416 rv.d *= bigtens[j];
417 if ((z = word0(rv) & Exp_mask)
418 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
419 goto ovfl;
420 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
421 /* set to largest number */
422 /* (Can't trust DBL_MAX) */
423 word0(rv) = Big0;
424 word1(rv) = Big1;
425 }
426 else
427 word0(rv) += P*Exp_msk1;
428 }
429
430 }
431 }
432 else if (e1 < 0) {
433 e1 = -e1;
434 if (i = e1 & 15)
435 rv.d /= tens[i];
436 if (e1 &= ~15) {
437 e1 >>= 4;
438 if (e1 >= 1 << n_bigtens)
439 goto undfl;
440 for(j = 0; e1 > 1; j++, e1 >>= 1)
441 if (e1 & 1)
442 rv.d *= tinytens[j];
443 /* The last multiplication could underflow. */
444 rv0.d = rv.d;
445 rv.d *= tinytens[j];
446 if (rv.d == 0) {
447 rv.d = 2.*rv0.d;
448 rv.d *= tinytens[j];
449 if (rv.d == 0) {
450 undfl:
451 rv.d = 0.;
452 errno = ERANGE;
453 if (bd0)
454 goto retfree;
455 goto ret;
456 }
457 word0(rv) = Tiny0;
458 word1(rv) = Tiny1;
459 /* The refinement below will clean
460 * this approximation up.
461 */
462 }
463 }
464 }
465
466 /* Now the hard part -- adjusting rv to the correct value.*/
467
468 /* Put digits into bd: true value = bd * 10^e */
469
470 bd0 = s2b(s0, nd0, nd, y);
471
472 for(;;) {
473 bd = Balloc(bd0->k);
474 Bcopy(bd, bd0);
475 bb = d2b(rv.d, &bbe, &bbbits); /* rv = bb * 2^bbe */
476 bs = i2b(1);
477
478 if (e >= 0) {
479 bb2 = bb5 = 0;
480 bd2 = bd5 = e;
481 }
482 else {
483 bb2 = bb5 = -e;
484 bd2 = bd5 = 0;
485 }
486 if (bbe >= 0)
487 bb2 += bbe;
488 else
489 bd2 -= bbe;
490 bs2 = bb2;
491 #ifdef Sudden_Underflow
492 #ifdef IBM
493 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
494 #else
495 j = P + 1 - bbbits;
496 #endif
497 #else
498 i = bbe + bbbits - 1; /* logb(rv) */
499 if (i < Emin) /* denormal */
500 j = bbe + (P-Emin);
501 else
502 j = P + 1 - bbbits;
503 #endif
504 bb2 += j;
505 bd2 += j;
506 i = bb2 < bd2 ? bb2 : bd2;
507 if (i > bs2)
508 i = bs2;
509 if (i > 0) {
510 bb2 -= i;
511 bd2 -= i;
512 bs2 -= i;
513 }
514 if (bb5 > 0) {
515 bs = pow5mult(bs, bb5);
516 bb1 = mult(bs, bb);
517 Bfree(bb);
518 bb = bb1;
519 }
520 if (bb2 > 0)
521 bb = lshift(bb, bb2);
522 if (bd5 > 0)
523 bd = pow5mult(bd, bd5);
524 if (bd2 > 0)
525 bd = lshift(bd, bd2);
526 if (bs2 > 0)
527 bs = lshift(bs, bs2);
528 delta = diff(bb, bd);
529 dsign = delta->sign;
530 delta->sign = 0;
531 i = cmp(delta, bs);
532 if (i < 0) {
533 /* Error is less than half an ulp -- check for
534 * special case of mantissa a power of two.
535 */
536 if (dsign || word1(rv) || word0(rv) & Bndry_mask)
537 break;
538 delta = lshift(delta,Log2P);
539 if (cmp(delta, bs) > 0)
540 goto drop_down;
541 break;
542 }
543 if (i == 0) {
544 /* exactly half-way between */
545 if (dsign) {
546 if ((word0(rv) & Bndry_mask1) == Bndry_mask1
547 && word1(rv) == 0xffffffff) {
548 /*boundary case -- increment exponent*/
549 word0(rv) = (word0(rv) & Exp_mask)
550 + Exp_msk1
551 #ifdef IBM
552 | Exp_msk1 >> 4
553 #endif
554 ;
555 word1(rv) = 0;
556 break;
557 }
558 }
559 else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
560 drop_down:
561 /* boundary case -- decrement exponent */
562 #ifdef Sudden_Underflow
563 L = word0(rv) & Exp_mask;
564 #ifdef IBM
565 if (L < Exp_msk1)
566 #else
567 if (L <= Exp_msk1)
568 #endif
569 goto undfl;
570 L -= Exp_msk1;
571 #else
572 L = (word0(rv) & Exp_mask) - Exp_msk1;
573 #endif
574 word0(rv) = L | Bndry_mask1;
575 word1(rv) = 0xffffffff;
576 #ifdef IBM
577 goto cont;
578 #else
579 break;
580 #endif
581 }
582 #ifndef ROUND_BIASED
583 if (!(word1(rv) & LSB))
584 break;
585 #endif
586 if (dsign)
587 rv.d += ulp(rv.d);
588 #ifndef ROUND_BIASED
589 else {
590 rv.d -= ulp(rv.d);
591 #ifndef Sudden_Underflow
592 if (rv.d == 0)
593 goto undfl;
594 #endif
595 }
596 #endif
597 break;
598 }
599 if ((aadj = ratio(delta, bs)) <= 2.) {
600 if (dsign)
601 aadj = aadj1 = 1.;
602 else if (word1(rv) || word0(rv) & Bndry_mask) {
603 #ifndef Sudden_Underflow
604 if (word1(rv) == Tiny1 && !word0(rv))
605 goto undfl;
606 #endif
607 aadj = 1.;
608 aadj1 = -1.;
609 }
610 else {
611 /* special case -- power of FLT_RADIX to be */
612 /* rounded down... */
613
614 if (aadj < 2./FLT_RADIX)
615 aadj = 1./FLT_RADIX;
616 else
617 aadj *= 0.5;
618 aadj1 = -aadj;
619 }
620 }
621 else {
622 aadj *= 0.5;
623 aadj1 = dsign ? aadj : -aadj;
624 #ifdef Check_FLT_ROUNDS
625 switch(FLT_ROUNDS) {
626 case 2: /* towards +infinity */
627 aadj1 -= 0.5;
628 break;
629 case 0: /* towards 0 */
630 case 3: /* towards -infinity */
631 aadj1 += 0.5;
632 }
633 #else
634 if (FLT_ROUNDS == 0)
635 aadj1 += 0.5;
636 #endif
637 }
638 y = word0(rv) & Exp_mask;
639
640 /* Check for overflow */
641
642 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
643 rv0.d = rv.d;
644 word0(rv) -= P*Exp_msk1;
645 adj = aadj1 * ulp(rv.d);
646 rv.d += adj;
647 if ((word0(rv) & Exp_mask) >=
648 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
649 if (word0(rv0) == Big0 && word1(rv0) == Big1)
650 goto ovfl;
651 word0(rv) = Big0;
652 word1(rv) = Big1;
653 goto cont;
654 }
655 else
656 word0(rv) += P*Exp_msk1;
657 }
658 else {
659 #ifdef Sudden_Underflow
660 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
661 rv0.d = rv.d;
662 word0(rv) += P*Exp_msk1;
663 adj = aadj1 * ulp(rv.d);
664 rv.d += adj;
665 #ifdef IBM
666 if ((word0(rv) & Exp_mask) < P*Exp_msk1)
667 #else
668 if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
669 #endif
670 {
671 if (word0(rv0) == Tiny0
672 && word1(rv0) == Tiny1)
673 goto undfl;
674 word0(rv) = Tiny0;
675 word1(rv) = Tiny1;
676 goto cont;
677 }
678 else
679 word0(rv) -= P*Exp_msk1;
680 }
681 else {
682 adj = aadj1 * ulp(rv.d);
683 rv.d += adj;
684 }
685 #else
686 /* Compute adj so that the IEEE rounding rules will
687 * correctly round rv + adj in some half-way cases.
688 * If rv * ulp(rv) is denormalized (i.e.,
689 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
690 * trouble from bits lost to denormalization;
691 * example: 1.2e-307 .
692 */
693 if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
694 aadj1 = (double)(int)(aadj + 0.5);
695 if (!dsign)
696 aadj1 = -aadj1;
697 }
698 adj = aadj1 * ulp(rv.d);
699 rv.d += adj;
700 #endif
701 }
702 z = word0(rv) & Exp_mask;
703 if (y == z) {
704 /* Can we stop now? */
705 L = aadj;
706 aadj -= L;
707 /* The tolerances below are conservative. */
708 if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
709 if (aadj < .4999999 || aadj > .5000001)
710 break;
711 }
712 else if (aadj < .4999999/FLT_RADIX)
713 break;
714 }
715 cont:
716 Bfree(bb);
717 Bfree(bd);
718 Bfree(bs);
719 Bfree(delta);
720 }
721 retfree:
722 Bfree(bb);
723 Bfree(bd);
724 Bfree(bs);
725 Bfree(bd0);
726 Bfree(delta);
727 ret:
728 if (se)
729 *se = (char *)s;
730 return sign ? -rv.d : rv.d;
731 }
732