strtod.c (9133B) download
1/* Copyright (c) 2002-2006 Lucent Technologies; see LICENSE */
2#include "common.h"
3#include "fmt.h"
4#include "fmtdef.h"
5#include "plan9.h"
6
7#include <ctype.h>
8#include <errno.h>
9#include <math.h>
10#include <stdlib.h>
11#include <string.h>
12
13static ulong
14umuldiv(ulong a, ulong b, ulong c) {
15 double d;
16
17 d = ((double) a * (double) b) / (double) c;
18 if (d >= 4294967295.)
19 d = 4294967295.;
20 return (ulong) d;
21}
22
23/*
24 * This routine will convert to arbitrary precision
25 * floating point entirely in multi-precision fixed.
26 * The answer is the closest floating point number to
27 * the given decimal number. Exactly half way are
28 * rounded ala ieee rules.
29 * Method is to scale input decimal between .500 and .999...
30 * with external power of 2, then binary search for the
31 * closest mantissa to this decimal number.
32 * Nmant is is the required precision. (53 for ieee dp)
33 * Nbits is the max number of bits/word. (must be <= 28)
34 * Prec is calculated - the number of words of fixed mantissa.
35 */
36enum {
37 Nbits = 28, /* bits safely represented in a ulong */
38 Nmant = 53, /* bits of precision required */
39 Prec = (Nmant + Nbits + 1) / Nbits, /* words of Nbits each to represent mantissa */
40 Sigbit = 1 << (Prec * Nbits - Nmant), /* first significant bit of Prec-th word */
41 Ndig = 1500,
42 One = (ulong) (1 << Nbits),
43 Half = (ulong) (One >> 1),
44 Maxe = 310,
45
46 Fsign = 1 << 0, /* found - */
47 Fesign = 1 << 1, /* found e- */
48 Fdpoint = 1 << 2, /* found . */
49
50 S0 = 0, /* _ _S0 +S1 #S2 .S3 */
51 S1, /* _+ #S2 .S3 */
52 S2, /* _+# #S2 .S4 eS5 */
53 S3, /* _+. #S4 */
54 S4, /* _+#.# #S4 eS5 */
55 S5, /* _+#.#e +S6 #S7 */
56 S6, /* _+#.#e+ #S7 */
57 S7 /* _+#.#e+# #S7 */
58};
59
60static int xcmp(char*, char*);
61static int fpcmp(char*, ulong*);
62static void frnorm(ulong*);
63static void divascii(char*, int*, int*, int*);
64static void mulascii(char*, int*, int*, int*);
65
66typedef struct Tab Tab;
67struct Tab {
68 int bp;
69 int siz;
70 char* cmp;
71};
72
73double
74fmtstrtod(const char* as, char** aas) {
75 int na, ex, dp, bp, c, i, flag, state;
76 ulong low[Prec], hig[Prec], mid[Prec];
77 double d;
78 char * s, a[Ndig];
79
80 flag = 0; /* Fsign, Fesign, Fdpoint */
81 na = 0; /* number of digits of a[] */
82 dp = 0; /* na of decimal point */
83 ex = 0; /* exonent */
84
85 state = S0;
86 for (s = (char*) as;; s++) {
87 c = *s;
88 if (c >= '0' && c <= '9') {
89 switch (state) {
90 case S0:
91 case S1:
92 case S2:
93 state = S2;
94 break;
95 case S3:
96 case S4:
97 state = S4;
98 break;
99
100 case S5:
101 case S6:
102 case S7:
103 state = S7;
104 ex = ex * 10 + (c - '0');
105 continue;
106 }
107 if (na == 0 && c == '0') {
108 dp--;
109 continue;
110 }
111 if (na < Ndig - 50)
112 a[na++] = c;
113 continue;
114 }
115 switch (c) {
116 case '\t':
117 case '\n':
118 case '\v':
119 case '\f':
120 case '\r':
121 case ' ':
122 if (state == S0)
123 continue;
124 break;
125 case '-':
126 if (state == S0)
127 flag |= Fsign;
128 else
129 flag |= Fesign;
130 FALLTHROUGH;
131 case '+':
132 if (state == S0)
133 state = S1;
134 else if (state == S5)
135 state = S6;
136 else
137 break; /* syntax */
138 continue;
139 case '.':
140 flag |= Fdpoint;
141 dp = na;
142 if (state == S0 || state == S1) {
143 state = S3;
144 continue;
145 }
146 if (state == S2) {
147 state = S4;
148 continue;
149 }
150 break;
151 case 'e':
152 case 'E':
153 if (state == S2 || state == S4) {
154 state = S5;
155 continue;
156 }
157 break;
158 }
159 break;
160 }
161
162 /*
163 * clean up return char-pointer
164 */
165 switch (state) {
166 case S0:
167 if (xcmp(s, "nan") == 0) {
168 if (aas != nil)
169 *aas = s + 3;
170 goto retnan;
171 }
172 FALLTHROUGH;
173 case S1:
174 if (xcmp(s, "infinity") == 0) {
175 if (aas != nil)
176 *aas = s + 8;
177 goto retinf;
178 }
179 if (xcmp(s, "inf") == 0) {
180 if (aas != nil)
181 *aas = s + 3;
182 goto retinf;
183 }
184 FALLTHROUGH;
185 case S3:
186 if (aas != nil)
187 *aas = (char*) as;
188 goto ret0; /* no digits found */
189 FALLTHROUGH;
190 case S6:
191 s--; /* back over +- */
192 FALLTHROUGH;
193 case S5:
194 s--; /* back over e */
195 break;
196 }
197 if (aas != nil)
198 *aas = s;
199
200 if (flag & Fdpoint)
201 while (na > 0 && a[na - 1] == '0')
202 na--;
203 if (na == 0)
204 goto ret0; /* zero */
205 a[na] = 0;
206 if (!(flag & Fdpoint))
207 dp = na;
208 if (flag & Fesign)
209 ex = -ex;
210 dp += ex;
211 if (dp < -Maxe) {
212 errno = ERANGE;
213 goto ret0; /* underflow by exp */
214 } else if (dp > +Maxe)
215 goto retinf; /* overflow by exp */
216
217 /*
218 * normalize the decimal ascii number
219 * to range .[5-9][0-9]* e0
220 */
221 bp = 0; /* binary exponent */
222 while (dp > 0)
223 divascii(a, &na, &dp, &bp);
224 while (dp < 0 || a[0] < '5')
225 mulascii(a, &na, &dp, &bp);
226
227 /* close approx by naive conversion */
228 mid[0] = 0;
229 mid[1] = 1;
230 for (i = 0; (c = a[i]) != '\0'; i++) {
231 mid[0] = mid[0] * 10 + (c - '0');
232 mid[1] = mid[1] * 10;
233 if (i >= 8)
234 break;
235 }
236 low[0] = umuldiv(mid[0], One, mid[1]);
237 hig[0] = umuldiv(mid[0] + 1, One, mid[1]);
238 for (i = 1; i < Prec; i++) {
239 low[i] = 0;
240 hig[i] = One - 1;
241 }
242
243 /* binary search for closest mantissa */
244 for (;;) {
245 /* mid = (hig + low) / 2 */
246 c = 0;
247 for (i = 0; i < Prec; i++) {
248 mid[i] = hig[i] + low[i];
249 if (c)
250 mid[i] += One;
251 c = mid[i] & 1;
252 mid[i] >>= 1;
253 }
254 frnorm(mid);
255
256 /* compare */
257 c = fpcmp(a, mid);
258 if (c > 0) {
259 c = 1;
260 for (i = 0; i < Prec; i++)
261 if (low[i] != mid[i]) {
262 c = 0;
263 low[i] = mid[i];
264 }
265 if (c)
266 break; /* between mid and hig */
267 continue;
268 }
269 if (c < 0) {
270 for (i = 0; i < Prec; i++)
271 hig[i] = mid[i];
272 continue;
273 }
274
275 /* only hard part is if even/odd roundings wants to go up */
276 c = mid[Prec - 1] & (Sigbit - 1);
277 if (c == Sigbit / 2 && (mid[Prec - 1] & Sigbit) == 0)
278 mid[Prec - 1] -= c;
279 break; /* exactly mid */
280 }
281
282 /* normal rounding applies */
283 c = mid[Prec - 1] & (Sigbit - 1);
284 mid[Prec - 1] -= c;
285 if (c >= Sigbit / 2) {
286 mid[Prec - 1] += Sigbit;
287 frnorm(mid);
288 }
289 goto out;
290
291ret0:
292 return 0;
293
294retnan:
295 return __NaN();
296
297retinf:
298 /*
299 * Unix strtod requires these. Plan 9 would return Inf(0) or Inf(-1). */
300 errno = ERANGE;
301 if (flag & Fsign)
302 return -HUGE_VAL;
303 return HUGE_VAL;
304
305out:
306 d = 0;
307 for (i = 0; i < Prec; i++)
308 d = d * One + mid[i];
309 if (flag & Fsign)
310 d = -d;
311 d = ldexp(d, bp - Prec * Nbits);
312 if (d == 0) { /* underflow */
313 errno = ERANGE;
314 }
315 return d;
316}
317
318static void
319frnorm(ulong* f) {
320 int i, c;
321
322 c = 0;
323 for (i = Prec - 1; i > 0; i--) {
324 f[i] += c;
325 c = f[i] >> Nbits;
326 f[i] &= One - 1;
327 }
328 f[0] += c;
329}
330
331static int
332fpcmp(char* a, ulong* f) {
333 ulong tf[Prec];
334 int i, d, c;
335
336 for (i = 0; i < Prec; i++)
337 tf[i] = f[i];
338
339 for (;;) {
340 /* tf *= 10 */
341 for (i = 0; i < Prec; i++)
342 tf[i] = tf[i] * 10;
343 frnorm(tf);
344 d = (tf[0] >> Nbits) + '0';
345 tf[0] &= One - 1;
346
347 /* compare next digit */
348 c = *a;
349 if (c == 0) {
350 if ('0' < d)
351 return -1;
352 if (tf[0] != 0)
353 goto cont;
354 for (i = 1; i < Prec; i++)
355 if (tf[i] != 0)
356 goto cont;
357 return 0;
358 }
359 if (c > d)
360 return +1;
361 if (c < d)
362 return -1;
363 a++;
364 cont:;
365 }
366}
367
368static void
369divby(char* a, int* na, int b) {
370 int n, c;
371 char* p;
372
373 p = a;
374 n = 0;
375 while (n >> b == 0) {
376 c = *a++;
377 if (c == 0) {
378 while (n) {
379 c = n * 10;
380 if (c >> b)
381 break;
382 n = c;
383 }
384 goto xx;
385 }
386 n = n * 10 + c - '0';
387 (*na)--;
388 }
389 for (;;) {
390 c = n >> b;
391 n -= c << b;
392 *p++ = c + '0';
393 c = *a++;
394 if (c == 0)
395 break;
396 n = n * 10 + c - '0';
397 }
398 (*na)++;
399xx:
400 while (n) {
401 n = n * 10;
402 c = n >> b;
403 n -= c << b;
404 *p++ = c + '0';
405 (*na)++;
406 }
407 *p = 0;
408}
409
410static Tab tab1[] = {
411 { 1, 0, "" },
412 { 3, 1, "7" },
413 { 6, 2, "63" },
414 { 9, 3, "511" },
415 { 13, 4, "8191" },
416 { 16, 5, "65535" },
417 { 19, 6, "524287" },
418 { 23, 7, "8388607" },
419 { 26, 8, "67108863" },
420 { 27, 9, "134217727" }
421};
422
423static void
424divascii(char* a, int* na, int* dp, int* bp) {
425 int b, d;
426 Tab* t;
427
428 d = *dp;
429 if (d >= (int) (nelem(tab1)))
430 d = (int) (nelem(tab1)) - 1;
431 t = tab1 + d;
432 b = t->bp;
433 if (memcmp(a, t->cmp, t->siz) > 0)
434 d--;
435 *dp -= d;
436 *bp += b;
437 divby(a, na, b);
438}
439
440static void
441mulby(char* a, char* p, char* q, int b) {
442 int n, c;
443
444 n = 0;
445 *p = 0;
446 for (;;) {
447 q--;
448 if (q < a)
449 break;
450 c = *q - '0';
451 c = (c << b) + n;
452 n = c / 10;
453 c -= n * 10;
454 p--;
455 *p = c + '0';
456 }
457 while (n) {
458 c = n;
459 n = c / 10;
460 c -= n * 10;
461 p--;
462 *p = c + '0';
463 }
464}
465
466static Tab tab2[] = {
467 { 1, 1, "" }, /* dp = 0-0 */
468 { 3, 3, "125" },
469 { 6, 5, "15625" },
470 { 9, 7, "1953125" },
471 { 13, 10, "1220703125" },
472 { 16, 12, "152587890625" },
473 { 19, 14, "19073486328125" },
474 { 23, 17, "11920928955078125" },
475 { 26, 19, "1490116119384765625" },
476 { 27, 19, "7450580596923828125" }, /* dp 8-9 */
477};
478
479static void
480mulascii(char* a, int* na, int* dp, int* bp) {
481 char* p;
482 int d, b;
483 Tab* t;
484
485 d = -*dp;
486 if (d >= (int) (nelem(tab2)))
487 d = (int) (nelem(tab2)) - 1;
488 t = tab2 + d;
489 b = t->bp;
490 if (memcmp(a, t->cmp, t->siz) < 0)
491 d--;
492 p = a + *na;
493 *bp -= b;
494 *dp += d;
495 *na += d;
496 mulby(a, p + d, p, b);
497}
498
499static int
500xcmp(char* a, char* b) {
501 int c1, c2;
502
503 while ((c1 = *b++) != '\0') {
504 c2 = *a++;
505 if (isupper(c2))
506 c2 = tolower(c2);
507 if (c1 != c2)
508 return 1;
509 }
510 return 0;
511}