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
11  * ANY 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 "plan9.h"
21 #include "fmt.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]; 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++) {
525 		c2 = *a++;
526 		if(isupper(c2))
527 			c2 = tolower(c2);
528 		if(c1 != c2)
529 			return 1;
530 	}
531 	return 0;
532 }