tstrtod.c - plan9port - [fork] Plan 9 from user space
HTML git clone git://src.adamsgaard.dk/plan9port
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
tstrtod.c (8562B)
---
1 /* Copyright (c) 2002-2006 Lucent Technologies; see LICENSE */
2 #include <stdlib.h>
3 #include <math.h>
4 #include <ctype.h>
5 #include <stdlib.h>
6 #include <string.h>
7 #include <errno.h>
8 #include "plan9.h"
9 #include "fmt.h"
10 #include "fmtdef.h"
11
12 static ulong
13 umuldiv(ulong a, ulong b, ulong c)
14 {
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 */
36 enum
37 {
38 Nbits = 28, /* bits safely represented in a ulong */
39 Nmant = 53, /* bits of precision required */
40 Prec = (Nmant+Nbits+1)/Nbits, /* words of Nbits each to represent mantissa */
41 Sigbit = 1<<(Prec*Nbits-Nmant), /* first significant bit of Prec-th word */
42 Ndig = 1500,
43 One = (ulong)(1<<Nbits),
44 Half = (ulong)(One>>1),
45 Maxe = 310,
46
47 Fsign = 1<<0, /* found - */
48 Fesign = 1<<1, /* found e- */
49 Fdpoint = 1<<2, /* found . */
50
51 S0 = 0, /* _ _S0 +S1 #S2 .S3 */
52 S1, /* _+ #S2 .S3 */
53 S2, /* _+# #S2 .S4 eS5 */
54 S3, /* _+. #S4 */
55 S4, /* _+#.# #S4 eS5 */
56 S5, /* _+#.#e +S6 #S7 */
57 S6, /* _+#.#e+ #S7 */
58 S7 /* _+#.#e+# #S7 */
59 };
60
61 static int xcmp(char*, char*);
62 static int fpcmp(char*, ulong*);
63 static void frnorm(ulong*);
64 static void divascii(char*, int*, int*, int*);
65 static void mulascii(char*, int*, int*, int*);
66
67 typedef struct Tab Tab;
68 struct Tab
69 {
70 int bp;
71 int siz;
72 char* cmp;
73 };
74
75 double
76 fmtstrtod(const char *as, char **aas)
77 {
78 int na, ex, dp, bp, c, i, flag, state;
79 ulong low[Prec], hig[Prec], mid[Prec];
80 double d;
81 char *s, a[Ndig];
82
83 flag = 0; /* Fsign, Fesign, Fdpoint */
84 na = 0; /* number of digits of a[] */
85 dp = 0; /* na of decimal point */
86 ex = 0; /* exonent */
87
88 state = S0;
89 for(s=(char*)as;; s++) {
90 c = *s;
91 if(c >= '0' && c <= '9') {
92 switch(state) {
93 case S0:
94 case S1:
95 case S2:
96 state = S2;
97 break;
98 case S3:
99 case S4:
100 state = S4;
101 break;
102
103 case S5:
104 case S6:
105 case S7:
106 state = S7;
107 ex = ex*10 + (c-'0');
108 continue;
109 }
110 if(na == 0 && c == '0') {
111 dp--;
112 continue;
113 }
114 if(na < Ndig-50)
115 a[na++] = c;
116 continue;
117 }
118 switch(c) {
119 case '\t':
120 case '\n':
121 case '\v':
122 case '\f':
123 case '\r':
124 case ' ':
125 if(state == S0)
126 continue;
127 break;
128 case '-':
129 if(state == S0)
130 flag |= Fsign;
131 else
132 flag |= Fesign;
133 case '+':
134 if(state == S0)
135 state = S1;
136 else
137 if(state == S5)
138 state = S6;
139 else
140 break; /* syntax */
141 continue;
142 case '.':
143 flag |= Fdpoint;
144 dp = na;
145 if(state == S0 || state == S1) {
146 state = S3;
147 continue;
148 }
149 if(state == S2) {
150 state = S4;
151 continue;
152 }
153 break;
154 case 'e':
155 case 'E':
156 if(state == S2 || state == S4) {
157 state = S5;
158 continue;
159 }
160 break;
161 }
162 break;
163 }
164
165 /*
166 * clean up return char-pointer
167 */
168 switch(state) {
169 case S0:
170 if(xcmp(s, "nan") == 0) {
171 if(aas != nil)
172 *aas = s+3;
173 goto retnan;
174 }
175 case S1:
176 if(xcmp(s, "infinity") == 0) {
177 if(aas != nil)
178 *aas = s+8;
179 goto retinf;
180 }
181 if(xcmp(s, "inf") == 0) {
182 if(aas != nil)
183 *aas = s+3;
184 goto retinf;
185 }
186 case S3:
187 if(aas != nil)
188 *aas = (char*)as;
189 goto ret0; /* no digits found */
190 case S6:
191 s--; /* back over +- */
192 case S5:
193 s--; /* back over e */
194 break;
195 }
196 if(aas != nil)
197 *aas = s;
198
199 if(flag & Fdpoint)
200 while(na > 0 && a[na-1] == '0')
201 na--;
202 if(na == 0)
203 goto ret0; /* zero */
204 a[na] = 0;
205 if(!(flag & Fdpoint))
206 dp = na;
207 if(flag & Fesign)
208 ex = -ex;
209 dp += ex;
210 if(dp < -Maxe){
211 errno = ERANGE;
212 goto ret0; /* underflow by exp */
213 } else
214 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
291 ret0:
292 return 0;
293
294 retnan:
295 return __NaN();
296
297 retinf:
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
305 out:
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
318 static void
319 frnorm(ulong *f)
320 {
321 int i, c;
322
323 c = 0;
324 for(i=Prec-1; i>0; i--) {
325 f[i] += c;
326 c = f[i] >> Nbits;
327 f[i] &= One-1;
328 }
329 f[0] += c;
330 }
331
332 static int
333 fpcmp(char *a, ulong* f)
334 {
335 ulong tf[Prec];
336 int i, d, c;
337
338 for(i=0; i<Prec; i++)
339 tf[i] = f[i];
340
341 for(;;) {
342 /* tf *= 10 */
343 for(i=0; i<Prec; i++)
344 tf[i] = tf[i]*10;
345 frnorm(tf);
346 d = (tf[0] >> Nbits) + '0';
347 tf[0] &= One-1;
348
349 /* compare next digit */
350 c = *a;
351 if(c == 0) {
352 if('0' < d)
353 return -1;
354 if(tf[0] != 0)
355 goto cont;
356 for(i=1; i<Prec; i++)
357 if(tf[i] != 0)
358 goto cont;
359 return 0;
360 }
361 if(c > d)
362 return +1;
363 if(c < d)
364 return -1;
365 a++;
366 cont:;
367 }
368 }
369
370 static void
371 divby(char *a, int *na, int b)
372 {
373 int n, c;
374 char *p;
375
376 p = a;
377 n = 0;
378 while(n>>b == 0) {
379 c = *a++;
380 if(c == 0) {
381 while(n) {
382 c = n*10;
383 if(c>>b)
384 break;
385 n = c;
386 }
387 goto xx;
388 }
389 n = n*10 + c-'0';
390 (*na)--;
391 }
392 for(;;) {
393 c = n>>b;
394 n -= c<<b;
395 *p++ = c + '0';
396 c = *a++;
397 if(c == 0)
398 break;
399 n = n*10 + c-'0';
400 }
401 (*na)++;
402 xx:
403 while(n) {
404 n = n*10;
405 c = n>>b;
406 n -= c<<b;
407 *p++ = c + '0';
408 (*na)++;
409 }
410 *p = 0;
411 }
412
413 static Tab tab1[] =
414 {
415 1, 0, "",
416 3, 1, "7",
417 6, 2, "63",
418 9, 3, "511",
419 13, 4, "8191",
420 16, 5, "65535",
421 19, 6, "524287",
422 23, 7, "8388607",
423 26, 8, "67108863",
424 27, 9, "134217727",
425 };
426
427 static void
428 divascii(char *a, int *na, int *dp, int *bp)
429 {
430 int b, d;
431 Tab *t;
432
433 d = *dp;
434 if(d >= (int)(nelem(tab1)))
435 d = (int)(nelem(tab1))-1;
436 t = tab1 + d;
437 b = t->bp;
438 if(memcmp(a, t->cmp, t->siz) > 0)
439 d--;
440 *dp -= d;
441 *bp += b;
442 divby(a, na, b);
443 }
444
445 static void
446 mulby(char *a, char *p, char *q, int b)
447 {
448 int n, c;
449
450 n = 0;
451 *p = 0;
452 for(;;) {
453 q--;
454 if(q < a)
455 break;
456 c = *q - '0';
457 c = (c<<b) + n;
458 n = c/10;
459 c -= n*10;
460 p--;
461 *p = c + '0';
462 }
463 while(n) {
464 c = n;
465 n = c/10;
466 c -= n*10;
467 p--;
468 *p = c + '0';
469 }
470 }
471
472 static Tab tab2[] =
473 {
474 1, 1, "", /* dp = 0-0 */
475 3, 3, "125",
476 6, 5, "15625",
477 9, 7, "1953125",
478 13, 10, "1220703125",
479 16, 12, "152587890625",
480 19, 14, "19073486328125",
481 23, 17, "11920928955078125",
482 26, 19, "1490116119384765625",
483 27, 19, "7450580596923828125", /* dp 8-9 */
484 };
485
486 static void
487 mulascii(char *a, int *na, int *dp, int *bp)
488 {
489 char *p;
490 int d, b;
491 Tab *t;
492
493 d = -*dp;
494 if(d >= (int)(nelem(tab2)))
495 d = (int)(nelem(tab2))-1;
496 t = tab2 + d;
497 b = t->bp;
498 if(memcmp(a, t->cmp, t->siz) < 0)
499 d--;
500 p = a + *na;
501 *bp -= b;
502 *dp += d;
503 *na += d;
504 mulby(a, p+d, p, b);
505 }
506
507 static int
508 xcmp(char *a, char *b)
509 {
510 int c1, c2;
511
512 while((c1 = *b++) != '\0') {
513 c2 = *a++;
514 if(isupper(c2))
515 c2 = tolower(c2);
516 if(c1 != c2)
517 return 1;
518 }
519 return 0;
520 }