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