![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to strtod.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10no / libc / stdio|
Return to strtod.c CVS log | Up to [Research Unix] / researchv10no / libc / stdio |
1.1 root 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
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 *
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
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
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
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: if (nd <= DBL_DIG
347: #ifndef RND_PRODQUOT
348: && FLT_ROUNDS == 1
349: #endif
350: ) {
351: if (!e)
352: goto ret;
353: if (e > 0) {
354: if (e <= Ten_pmax) {
355: #ifdef VAX
356: goto vax_ovfl_check;
357: #else
358: /* rv = */ rounded_product(rv.d, tens[e]);
359: goto ret;
360: #endif
361: }
362: i = DBL_DIG - nd;
363: if (e <= Ten_pmax + i) {
364: /* A fancier test would sometimes let us do
365: * this for larger i values.
366: */
367: e -= i;
368: rv.d *= tens[i];
369: #ifdef VAX
370: /* VAX exponent range is so narrow we must
371: * worry about overflow here...
372: */
373: vax_ovfl_check:
374: word0(rv) -= P*Exp_msk1;
375: /* rv = */ rounded_product(rv.d, tens[e]);
376: if ((word0(rv) & Exp_mask)
377: > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
378: goto ovfl;
379: word0(rv) += P*Exp_msk1;
380: #else
381: /* rv = */ rounded_product(rv.d, tens[e]);
382: #endif
383: goto ret;
384: }
385: }
386: else if (e >= -Ten_pmax) {
387: /* rv = */ rounded_quotient(rv.d, tens[-e]);
388: goto ret;
389: }
390: }
391: e1 += nd - k;
392:
393: /* Get starting approximation = rv * 10**e1 */
394:
395: if (e1 > 0) {
396: if (i = e1 & 15)
397: rv.d *= tens[i];
398: if (e1 &= ~15) {
399: if (e1 > DBL_MAX_10_EXP) {
400: ovfl:
401: errno = ERANGE;
402: rv.d = HUGE_VAL;
403: goto ret;
404: }
405: if (e1 >>= 4) {
406: for(j = 0; e1 > 1; j++, e1 >>= 1)
407: if (e1 & 1)
408: rv.d *= bigtens[j];
409: /* The last multiplication could overflow. */
410: word0(rv) -= P*Exp_msk1;
411: rv.d *= bigtens[j];
412: if ((z = word0(rv) & Exp_mask)
413: > Exp_msk1*(DBL_MAX_EXP+Bias-P))
414: goto ovfl;
415: if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
416: /* set to largest number */
417: /* (Can't trust DBL_MAX) */
418: word0(rv) = Big0;
419: word1(rv) = Big1;
420: }
421: else
422: word0(rv) += P*Exp_msk1;
423: }
424:
425: }
426: }
427: else if (e1 < 0) {
428: e1 = -e1;
429: if (i = e1 & 15)
430: rv.d /= tens[i];
431: if (e1 &= ~15) {
432: e1 >>= 4;
433: for(j = 0; e1 > 1; j++, e1 >>= 1)
434: if (e1 & 1)
435: rv.d *= tinytens[j];
436: /* The last multiplication could underflow. */
437: rv0.d = rv.d;
438: rv.d *= tinytens[j];
439: if (rv.d == 0) {
440: rv.d = 2.*rv0.d;
441: rv.d *= tinytens[j];
442: if (rv.d == 0) {
443: undfl:
444: rv.d = 0.;
445: errno = ERANGE;
446: goto ret;
447: }
448: word0(rv) = Tiny0;
449: word1(rv) = Tiny1;
450: /* The refinement below will clean
451: * this approximation up.
452: */
453: }
454: }
455: }
456:
457: /* Now the hard part -- adjusting rv to the correct value.*/
458:
459: /* Put digits into bd: true value = bd * 10^e */
460:
461: bd0 = s2b(s0, nd0, nd, y);
462:
463: for(;;) {
464: bd = Balloc(bd0->k);
465: Bcopy(bd, bd0);
466: bb = d2b(rv.d, &bbe, &bbbits); /* rv = bb * 2^bbe */
467: bs = i2b(1);
468:
469: if (e >= 0) {
470: bb2 = bb5 = 0;
471: bd2 = bd5 = e;
472: }
473: else {
474: bb2 = bb5 = -e;
475: bd2 = bd5 = 0;
476: }
477: if (bbe >= 0)
478: bb2 += bbe;
479: else
480: bd2 -= bbe;
481: bs2 = bb2;
482: #ifdef Sudden_Underflow
483: #ifdef IBM
484: j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
485: #else
486: j = P + 1 - bbbits;
487: #endif
488: #else
489: i = bbe + bbbits - 1; /* logb(rv) */
490: if (i < Emin) /* denormal */
491: j = bbe + (P-Emin);
492: else
493: j = P + 1 - bbbits;
494: #endif
495: bb2 += j;
496: bd2 += j;
497: i = bb2 < bd2 ? bb2 : bd2;
498: if (i > bs2)
499: i = bs2;
500: if (i > 0) {
501: bb2 -= i;
502: bd2 -= i;
503: bs2 -= i;
504: }
505: if (bb5 > 0) {
506: bs = pow5mult(bs, bb5);
507: bb1 = mult(bs, bb);
508: Bfree(bb);
509: bb = bb1;
510: }
511: if (bb2 > 0)
512: bb = lshift(bb, bb2);
513: if (bd5 > 0)
514: bd = pow5mult(bd, bd5);
515: if (bd2 > 0)
516: bd = lshift(bd, bd2);
517: if (bs2 > 0)
518: bs = lshift(bs, bs2);
519: delta = diff(bb, bd);
520: dsign = delta->sign;
521: delta->sign = 0;
522: i = cmp(delta, bs);
523: if (i < 0) {
524: /* Error is less than half an ulp -- check for
525: * special case of mantissa a power of two.
526: */
527: if (dsign || word1(rv) || word0(rv) & Bndry_mask)
528: break;
529: delta = lshift(delta,Log2P);
530: if (cmp(delta, bs) > 0)
531: goto drop_down;
532: break;
533: }
534: if (i == 0) {
535: /* exactly half-way between */
536: if (dsign) {
537: if ((word0(rv) & Bndry_mask1) == Bndry_mask1
538: && word1(rv) == 0xffffffff) {
539: /*boundary case -- increment exponent*/
540: word0(rv) = (word0(rv) & Exp_mask)
541: + Exp_msk1
542: #ifdef IBM
543: | Exp_msk1 >> 4
544: #endif
545: ;
546: word1(rv) = 0;
547: break;
548: }
549: }
550: else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
551: drop_down:
552: /* boundary case -- decrement exponent */
553: #ifdef Sudden_Underflow
554: L = word0(rv) & Exp_mask;
555: #ifdef IBM
556: if (L < Exp_msk1)
557: #else
558: if (L <= Exp_msk1)
559: #endif
560: goto undfl;
561: L -= Exp_msk1;
562: #else
563: L = (word0(rv) & Exp_mask) - Exp_msk1;
564: #endif
565: word0(rv) = L | Bndry_mask1;
566: word1(rv) = 0xffffffff;
567: #ifdef IBM
568: goto cont;
569: #else
570: break;
571: #endif
572: }
573: #ifndef ROUND_BIASED
574: if (!(word1(rv) & LSB))
575: break;
576: #endif
577: if (dsign)
578: rv.d += ulp(rv.d);
579: #ifndef ROUND_BIASED
580: else {
581: rv.d -= ulp(rv.d);
582: #ifndef Sudden_Underflow
583: if (rv.d == 0)
584: goto undfl;
585: #endif
586: }
587: #endif
588: break;
589: }
590: if ((aadj = ratio(delta, bs)) <= 2.) {
591: if (dsign)
592: aadj = aadj1 = 1.;
593: else if (word1(rv) || word0(rv) & Bndry_mask) {
594: #ifndef Sudden_Underflow
595: if (word1(rv) == Tiny1 && !word0(rv))
596: goto undfl;
597: #endif
598: aadj = 1.;
599: aadj1 = -1.;
600: }
601: else {
602: /* special case -- power of FLT_RADIX to be */
603: /* rounded down... */
604:
605: if (aadj < 2./FLT_RADIX)
606: aadj = 1./FLT_RADIX;
607: else
608: aadj *= 0.5;
609: aadj1 = -aadj;
610: }
611: }
612: else {
613: aadj *= 0.5;
614: aadj1 = dsign ? aadj : -aadj;
615: #ifdef Check_FLT_ROUNDS
616: switch(FLT_ROUNDS) {
617: case 2: /* towards +infinity */
618: aadj1 -= 0.5;
619: break;
620: case 0: /* towards 0 */
621: case 3: /* towards -infinity */
622: aadj1 += 0.5;
623: }
624: #else
625: if (FLT_ROUNDS == 0)
626: aadj1 += 0.5;
627: #endif
628: }
629: y = word0(rv) & Exp_mask;
630:
631: /* Check for overflow */
632:
633: if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
634: rv0.d = rv.d;
635: word0(rv) -= P*Exp_msk1;
636: adj = aadj1 * ulp(rv.d);
637: rv.d += adj;
638: if ((word0(rv) & Exp_mask) >=
639: Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
640: if (word0(rv0) == Big0 && word1(rv0) == Big1)
641: goto ovfl;
642: word0(rv) = Big0;
643: word1(rv) = Big1;
644: goto cont;
645: }
646: else
647: word0(rv) += P*Exp_msk1;
648: }
649: else {
650: #ifdef Sudden_Underflow
651: if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
652: rv0.d = rv.d;
653: word0(rv) += P*Exp_msk1;
654: adj = aadj1 * ulp(rv.d);
655: rv.d += adj;
656: #ifdef IBM
657: if ((word0(rv) & Exp_mask) < P*Exp_msk1)
658: #else
659: if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
660: #endif
661: {
662: if (word0(rv0) == Tiny0
663: && word1(rv0) == Tiny1)
664: goto undfl;
665: word0(rv) = Tiny0;
666: word1(rv) = Tiny1;
667: goto cont;
668: }
669: else
670: word0(rv) -= P*Exp_msk1;
671: }
672: else {
673: adj = aadj1 * ulp(rv.d);
674: rv.d += adj;
675: }
676: #else
677: /* Compute adj so that the IEEE rounding rules will
678: * correctly round rv + adj in some half-way cases.
679: * If rv * ulp(rv) is denormalized (i.e.,
680: * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
681: * trouble from bits lost to denormalization;
682: * example: 1.2e-307 .
683: */
684: if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
685: aadj1 = (double)(int)(aadj + 0.5);
686: if (!dsign)
687: aadj1 = -aadj1;
688: }
689: adj = aadj1 * ulp(rv.d);
690: rv.d += adj;
691: #endif
692: }
693: z = word0(rv) & Exp_mask;
694: if (y == z) {
695: /* Can we stop now? */
696: L = aadj;
697: aadj -= L;
698: /* The tolerances below are conservative. */
699: if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
700: if (aadj < .4999999 || aadj > .5000001)
701: break;
702: }
703: else if (aadj < .4999999/FLT_RADIX)
704: break;
705: }
706: cont:
707: Bfree(bb);
708: Bfree(bd);
709: Bfree(bs);
710: Bfree(delta);
711: }
712: Bfree(bb);
713: Bfree(bd);
714: Bfree(bs);
715: Bfree(bd0);
716: Bfree(delta);
717: ret:
718: if (se)
719: *se = (char *)s;
720: return sign ? -rv.d : rv.d;
721: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.