1 #include <u.h>
2 #include <libc.h>
3 #include "fmtdef.h"
4
5 static ulong
umuldiv(ulong a,ulong b,ulong c)6 umuldiv(ulong a, ulong b, ulong c)
7 {
8 double d;
9
10 d = ((double)a * (double)b) / (double)c;
11 if(d >= 4294967295.)
12 d = 4294967295.;
13 return (ulong)d;
14 }
15
16 /*
17 * This routine will convert to arbitrary precision
18 * floating point entirely in multi-precision fixed.
19 * The answer is the closest floating point number to
20 * the given decimal number. Exactly half way are
21 * rounded ala ieee rules.
22 * Method is to scale input decimal between .500 and .999...
23 * with external power of 2, then binary search for the
24 * closest mantissa to this decimal number.
25 * Nmant is is the required precision. (53 for ieee dp)
26 * Nbits is the max number of bits/word. (must be <= 28)
27 * Prec is calculated - the number of words of fixed mantissa.
28 */
29 enum
30 {
31 Nbits = 28, /* bits safely represented in a ulong */
32 Nmant = 53, /* bits of precision required */
33 Prec = (Nmant+Nbits+1)/Nbits, /* words of Nbits each to represent mantissa */
34 Sigbit = 1<<(Prec*Nbits-Nmant), /* first significant bit of Prec-th word */
35 Ndig = 1500,
36 One = (ulong)(1<<Nbits),
37 Half = (ulong)(One>>1),
38 Maxe = 310,
39
40 Fsign = 1<<0, /* found - */
41 Fesign = 1<<1, /* found e- */
42 Fdpoint = 1<<2, /* found . */
43
44 S0 = 0, /* _ _S0 +S1 #S2 .S3 */
45 S1, /* _+ #S2 .S3 */
46 S2, /* _+# #S2 .S4 eS5 */
47 S3, /* _+. #S4 */
48 S4, /* _+#.# #S4 eS5 */
49 S5, /* _+#.#e +S6 #S7 */
50 S6, /* _+#.#e+ #S7 */
51 S7, /* _+#.#e+# #S7 */
52 };
53
54 static int xcmp(char*, char*);
55 static int fpcmp(char*, ulong*);
56 static void frnorm(ulong*);
57 static void divascii(char*, int*, int*, int*);
58 static void mulascii(char*, int*, int*, int*);
59
60 typedef struct Tab Tab;
61 struct Tab
62 {
63 int bp;
64 int siz;
65 char* cmp;
66 };
67
68 #ifndef ERANGE
69 #define ERANGE 12345
70 #endif
71
72 double
fmtstrtod(const char * as,char ** aas)73 fmtstrtod(const char *as, char **aas)
74 {
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 case '+':
131 if(state == S0)
132 state = S1;
133 else
134 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 case S1:
173 if(xcmp(s, "infinity") == 0) {
174 if(aas != nil)
175 *aas = s+8;
176 goto retinf;
177 }
178 if(xcmp(s, "inf") == 0) {
179 if(aas != nil)
180 *aas = s+3;
181 goto retinf;
182 }
183 case S3:
184 if(aas != nil)
185 *aas = (char*)as;
186 goto ret0; /* no digits found */
187 case S6:
188 s--; /* back over +- */
189 case S5:
190 s--; /* back over e */
191 break;
192 }
193 if(aas != nil)
194 *aas = s;
195
196 if(flag & Fdpoint)
197 while(na > 0 && a[na-1] == '0')
198 na--;
199 if(na == 0)
200 goto ret0; /* zero */
201 a[na] = 0;
202 if(!(flag & Fdpoint))
203 dp = na;
204 if(flag & Fesign)
205 ex = -ex;
206 dp += ex;
207 if(dp < -Maxe){
208 errno = ERANGE;
209 goto ret0; /* underflow by exp */
210 } else
211 if(dp > +Maxe)
212 goto retinf; /* overflow by exp */
213
214 /*
215 * normalize the decimal ascii number
216 * to range .[5-9][0-9]* e0
217 */
218 bp = 0; /* binary exponent */
219 while(dp > 0)
220 divascii(a, &na, &dp, &bp);
221 while(dp < 0 || a[0] < '5')
222 mulascii(a, &na, &dp, &bp);
223
224 /* close approx by naive conversion */
225 mid[0] = 0;
226 mid[1] = 1;
227 for(i=0; (c=a[i]); i++) {
228 mid[0] = mid[0]*10 + (c-'0');
229 mid[1] = mid[1]*10;
230 if(i >= 8)
231 break;
232 }
233 low[0] = umuldiv(mid[0], One, mid[1]);
234 hig[0] = umuldiv(mid[0]+1, One, mid[1]);
235 for(i=1; i<Prec; i++) {
236 low[i] = 0;
237 hig[i] = One-1;
238 }
239
240 /* binary search for closest mantissa */
241 for(;;) {
242 /* mid = (hig + low) / 2 */
243 c = 0;
244 for(i=0; i<Prec; i++) {
245 mid[i] = hig[i] + low[i];
246 if(c)
247 mid[i] += One;
248 c = mid[i] & 1;
249 mid[i] >>= 1;
250 }
251 frnorm(mid);
252
253 /* compare */
254 c = fpcmp(a, mid);
255 if(c > 0) {
256 c = 1;
257 for(i=0; i<Prec; i++)
258 if(low[i] != mid[i]) {
259 c = 0;
260 low[i] = mid[i];
261 }
262 if(c)
263 break; /* between mid and hig */
264 continue;
265 }
266 if(c < 0) {
267 for(i=0; i<Prec; i++)
268 hig[i] = mid[i];
269 continue;
270 }
271
272 /* only hard part is if even/odd roundings wants to go up */
273 c = mid[Prec-1] & (Sigbit-1);
274 if(c == Sigbit/2 && (mid[Prec-1]&Sigbit) == 0)
275 mid[Prec-1] -= c;
276 break; /* exactly mid */
277 }
278
279 /* normal rounding applies */
280 c = mid[Prec-1] & (Sigbit-1);
281 mid[Prec-1] -= c;
282 if(c >= Sigbit/2) {
283 mid[Prec-1] += Sigbit;
284 frnorm(mid);
285 }
286 goto out;
287
288 ret0:
289 return 0;
290
291 retnan:
292 return __NaN();
293
294 retinf:
295 /*
296 * Unix strtod requires these. Plan 9 would return Inf(0) or Inf(-1). */
297 errno = ERANGE;
298 if(flag & Fsign)
299 return -HUGE_VAL;
300 return HUGE_VAL;
301
302 out:
303 d = 0;
304 for(i=0; i<Prec; i++)
305 d = d*One + mid[i];
306 if(flag & Fsign)
307 d = -d;
308 d = ldexp(d, bp - Prec*Nbits);
309 if(d == 0){ /* underflow */
310 errno = ERANGE;
311 }
312 return d;
313 }
314
315 static void
frnorm(ulong * f)316 frnorm(ulong *f)
317 {
318 int i, c;
319
320 c = 0;
321 for(i=Prec-1; i>0; i--) {
322 f[i] += c;
323 c = f[i] >> Nbits;
324 f[i] &= One-1;
325 }
326 f[0] += c;
327 }
328
329 static int
fpcmp(char * a,ulong * f)330 fpcmp(char *a, ulong* f)
331 {
332 ulong tf[Prec];
333 int i, d, c;
334
335 for(i=0; i<Prec; i++)
336 tf[i] = f[i];
337
338 for(;;) {
339 /* tf *= 10 */
340 for(i=0; i<Prec; i++)
341 tf[i] = tf[i]*10;
342 frnorm(tf);
343 d = (tf[0] >> Nbits) + '0';
344 tf[0] &= One-1;
345
346 /* compare next digit */
347 c = *a;
348 if(c == 0) {
349 if('0' < d)
350 return -1;
351 if(tf[0] != 0)
352 goto cont;
353 for(i=1; i<Prec; i++)
354 if(tf[i] != 0)
355 goto cont;
356 return 0;
357 }
358 if(c > d)
359 return +1;
360 if(c < d)
361 return -1;
362 a++;
363 cont:;
364 }
365 }
366
367 static void
divby(char * a,int * na,int b)368 divby(char *a, int *na, int b)
369 {
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)++;
399 xx:
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
410 static Tab tab1[] =
411 {
412 1, 0, "",
413 3, 1, "7",
414 6, 2, "63",
415 9, 3, "511",
416 13, 4, "8191",
417 16, 5, "65535",
418 19, 6, "524287",
419 23, 7, "8388607",
420 26, 8, "67108863",
421 27, 9, "134217727",
422 };
423
424 static void
divascii(char * a,int * na,int * dp,int * bp)425 divascii(char *a, int *na, int *dp, int *bp)
426 {
427 int b, d;
428 Tab *t;
429
430 d = *dp;
431 if(d >= (int)(nelem(tab1)))
432 d = (int)(nelem(tab1))-1;
433 t = tab1 + d;
434 b = t->bp;
435 if(memcmp(a, t->cmp, t->siz) > 0)
436 d--;
437 *dp -= d;
438 *bp += b;
439 divby(a, na, b);
440 }
441
442 static void
mulby(char * a,char * p,char * q,int b)443 mulby(char *a, char *p, char *q, int b)
444 {
445 int n, c;
446
447 n = 0;
448 *p = 0;
449 for(;;) {
450 q--;
451 if(q < a)
452 break;
453 c = *q - '0';
454 c = (c<<b) + n;
455 n = c/10;
456 c -= n*10;
457 p--;
458 *p = c + '0';
459 }
460 while(n) {
461 c = n;
462 n = c/10;
463 c -= n*10;
464 p--;
465 *p = c + '0';
466 }
467 }
468
469 static Tab tab2[] =
470 {
471 1, 1, "", /* dp = 0-0 */
472 3, 3, "125",
473 6, 5, "15625",
474 9, 7, "1953125",
475 13, 10, "1220703125",
476 16, 12, "152587890625",
477 19, 14, "19073486328125",
478 23, 17, "11920928955078125",
479 26, 19, "1490116119384765625",
480 27, 19, "7450580596923828125", /* dp 8-9 */
481 };
482
483 static void
mulascii(char * a,int * na,int * dp,int * bp)484 mulascii(char *a, int *na, int *dp, int *bp)
485 {
486 char *p;
487 int d, b;
488 Tab *t;
489
490 d = -*dp;
491 if(d >= (int)(nelem(tab2)))
492 d = (int)(nelem(tab2))-1;
493 t = tab2 + d;
494 b = t->bp;
495 if(memcmp(a, t->cmp, t->siz) < 0)
496 d--;
497 p = a + *na;
498 *bp -= b;
499 *dp += d;
500 *na += d;
501 mulby(a, p+d, p, b);
502 }
503
504 static int
xcmp(char * a,char * b)505 xcmp(char *a, char *b)
506 {
507 int c1, c2;
508
509 while((c1 = *b++)) {
510 c2 = *a++;
511 if(isupper(c2))
512 c2 = tolower(c2);
513 if(c1 != c2)
514 return 1;
515 }
516 return 0;
517 }
518