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