1#include <stdint.h>
2#include <stdio.h>
3#include <math.h>
4#include <float.h>
5#include <limits.h>
6#include <errno.h>
7#include <ctype.h>
8
9#include "shgetc.h"
10#include "floatscan.h"
11
12#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
13
14#define LD_B1B_DIG 2
15#define LD_B1B_MAX 9007199, 254740991
16#define KMAX 128
17
18#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
19
20#define LD_B1B_DIG 3
21#define LD_B1B_MAX 18, 446744073, 709551615
22#define KMAX 2048
23
24#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
25
26#define LD_B1B_DIG 4
27#define LD_B1B_MAX 10384593, 717069655, 257060992, 658440191
28#define KMAX 2048
29
30#else
31#error Unsupported long double representation
32#endif
33
34#define MASK (KMAX-1)
35
36#define CONCAT2(x,y) x ## y
37#define CONCAT(x,y) CONCAT2(x,y)
38
39static long long scanexp(FILE *f, int pok)
40{
41	int c;
42	int x;
43	long long y;
44	int neg = 0;
45
46	c = shgetc(f);
47	if (c=='+' || c=='-') {
48		neg = (c=='-');
49		c = shgetc(f);
50		if (c-'0'>=10U && pok) shunget(f);
51	}
52	if (c-'0'>=10U) {
53		shunget(f);
54		return LLONG_MIN;
55	}
56	for (x=0; c-'0'<10U && x<INT_MAX/10; c = shgetc(f))
57		x = 10*x + c-'0';
58	for (y=x; c-'0'<10U && y<LLONG_MAX/100; c = shgetc(f))
59		y = 10*y + c-'0';
60	for (; c-'0'<10U; c = shgetc(f));
61	shunget(f);
62	return neg ? -y : y;
63}
64
65
66static long double decfloat(FILE *f, int c, int bits, int emin, int sign, int pok)
67{
68	uint32_t x[KMAX];
69	static const uint32_t th[] = { LD_B1B_MAX };
70	int i, j, k, a, z;
71	long long lrp=0, dc=0;
72	long long e10=0;
73	int lnz = 0;
74	int gotdig = 0, gotrad = 0;
75	int rp;
76	int e2;
77	int emax = -emin-bits+3;
78	int denormal = 0;
79	long double y;
80	long double frac=0;
81	long double bias=0;
82	static const int p10s[] = { 10, 100, 1000, 10000,
83		100000, 1000000, 10000000, 100000000 };
84
85	j=0;
86	k=0;
87
88	/* Don't let leading zeros consume buffer space */
89	for (; c=='0'; c = shgetc(f)) gotdig=1;
90	if (c=='.') {
91		gotrad = 1;
92		for (c = shgetc(f); c=='0'; c = shgetc(f)) gotdig=1, lrp--;
93	}
94
95	x[0] = 0;
96	for (; c-'0'<10U || c=='.'; c = shgetc(f)) {
97		if (c == '.') {
98			if (gotrad) break;
99			gotrad = 1;
100			lrp = dc;
101		} else if (k < KMAX-3) {
102			dc++;
103			if (c!='0') lnz = dc;
104			if (j) x[k] = x[k]*10 + c-'0';
105			else x[k] = c-'0';
106			if (++j==9) {
107				k++;
108				j=0;
109			}
110			gotdig=1;
111		} else {
112			dc++;
113			if (c!='0') {
114				lnz = (KMAX-4)*9;
115				x[KMAX-4] |= 1;
116			}
117		}
118	}
119	if (!gotrad) lrp=dc;
120
121	if (gotdig && (c|32)=='e') {
122		e10 = scanexp(f, pok);
123		if (e10 == LLONG_MIN) {
124			if (pok) {
125				shunget(f);
126			} else {
127				shlim(f, 0);
128				return 0;
129			}
130			e10 = 0;
131		}
132		lrp += e10;
133	} else if (c>=0) {
134		shunget(f);
135	}
136	if (!gotdig) {
137		errno = EINVAL;
138		shlim(f, 0);
139		return 0;
140	}
141
142	/* Handle zero specially to avoid nasty special cases later */
143	if (!x[0]) return sign * 0.0;
144
145	/* Optimize small integers (w/no exponent) and over/under-flow */
146	if (lrp==dc && dc<10 && (bits>30 || x[0]>>bits==0))
147		return sign * (long double)x[0];
148	if (lrp > -emin/2) {
149		errno = ERANGE;
150		return sign * LDBL_MAX * LDBL_MAX;
151	}
152	if (lrp < emin-2*LDBL_MANT_DIG) {
153		errno = ERANGE;
154		return sign * LDBL_MIN * LDBL_MIN;
155	}
156
157	/* Align incomplete final B1B digit */
158	if (j) {
159		for (; j<9; j++) x[k]*=10;
160		k++;
161		j=0;
162	}
163
164	a = 0;
165	z = k;
166	e2 = 0;
167	rp = lrp;
168
169	/* Optimize small to mid-size integers (even in exp. notation) */
170	if (lnz<9 && lnz<=rp && rp < 18) {
171		if (rp == 9) return sign * (long double)x[0];
172		if (rp < 9) return sign * (long double)x[0] / p10s[8-rp];
173		int bitlim = bits-3*(int)(rp-9);
174		if (bitlim>30 || x[0]>>bitlim==0)
175			return sign * (long double)x[0] * p10s[rp-10];
176	}
177
178	/* Drop trailing zeros */
179	for (; !x[z-1]; z--);
180
181	/* Align radix point to B1B digit boundary */
182	if (rp % 9) {
183		int rpm9 = rp>=0 ? rp%9 : rp%9+9;
184		int p10 = p10s[8-rpm9];
185		uint32_t carry = 0;
186		for (k=a; k!=z; k++) {
187			uint32_t tmp = x[k] % p10;
188			x[k] = x[k]/p10 + carry;
189			carry = 1000000000/p10 * tmp;
190			if (k==a && !x[k]) {
191				a = (a+1 & MASK);
192				rp -= 9;
193			}
194		}
195		if (carry) x[z++] = carry;
196		rp += 9-rpm9;
197	}
198
199	/* Upscale until desired number of bits are left of radix point */
200	while (rp < 9*LD_B1B_DIG || (rp == 9*LD_B1B_DIG && x[a]<th[0])) {
201		uint32_t carry = 0;
202		e2 -= 29;
203		for (k=(z-1 & MASK); ; k=(k-1 & MASK)) {
204			uint64_t tmp = ((uint64_t)x[k] << 29) + carry;
205			if (tmp > 1000000000) {
206				carry = tmp / 1000000000;
207				x[k] = tmp % 1000000000;
208			} else {
209				carry = 0;
210				x[k] = tmp;
211			}
212			if (k==(z-1 & MASK) && k!=a && !x[k]) z = k;
213			if (k==a) break;
214		}
215		if (carry) {
216			rp += 9;
217			a = (a-1 & MASK);
218			if (a == z) {
219				z = (z-1 & MASK);
220				x[z-1 & MASK] |= x[z];
221			}
222			x[a] = carry;
223		}
224	}
225
226	/* Downscale until exactly number of bits are left of radix point */
227	for (;;) {
228		uint32_t carry = 0;
229		int sh = 1;
230		for (i=0; i<LD_B1B_DIG; i++) {
231			k = (a+i & MASK);
232			if (k == z || x[k] < th[i]) {
233				i=LD_B1B_DIG;
234				break;
235			}
236			if (x[a+i & MASK] > th[i]) break;
237		}
238		if (i==LD_B1B_DIG && rp==9*LD_B1B_DIG) break;
239		/* FIXME: find a way to compute optimal sh */
240		if (rp > 9+9*LD_B1B_DIG) sh = 9;
241		e2 += sh;
242		for (k=a; k!=z; k=(k+1 & MASK)) {
243			uint32_t tmp = x[k] & (1<<sh)-1;
244			x[k] = (x[k]>>sh) + carry;
245			carry = (1000000000>>sh) * tmp;
246			if (k==a && !x[k]) {
247				a = (a+1 & MASK);
248				i--;
249				rp -= 9;
250			}
251		}
252		if (carry) {
253			if ((z+1 & MASK) != a) {
254				x[z] = carry;
255				z = (z+1 & MASK);
256			} else x[z-1 & MASK] |= 1;
257		}
258	}
259
260	/* Assemble desired bits into floating point variable */
261	for (y=i=0; i<LD_B1B_DIG; i++) {
262		if ((a+i & MASK)==z) x[(z=(z+1 & MASK))-1] = 0;
263		y = 1000000000.0L * y + x[a+i & MASK];
264	}
265
266	y *= sign;
267
268	/* Limit precision for denormal results */
269	if (bits > LDBL_MANT_DIG+e2-emin) {
270		bits = LDBL_MANT_DIG+e2-emin;
271		if (bits<0) bits=0;
272		denormal = 1;
273	}
274
275	/* Calculate bias term to force rounding, move out lower bits */
276	if (bits < LDBL_MANT_DIG) {
277		bias = copysignl(scalbn(1, 2*LDBL_MANT_DIG-bits-1), y);
278		frac = fmodl(y, scalbn(1, LDBL_MANT_DIG-bits));
279		y -= frac;
280		y += bias;
281	}
282
283	/* Process tail of decimal input so it can affect rounding */
284	if ((a+i & MASK) != z) {
285		uint32_t t = x[a+i & MASK];
286		if (t < 500000000 && (t || (a+i+1 & MASK) != z))
287			frac += 0.25*sign;
288		else if (t > 500000000)
289			frac += 0.75*sign;
290		else if (t == 500000000) {
291			if ((a+i+1 & MASK) == z)
292				frac += 0.5*sign;
293			else
294				frac += 0.75*sign;
295		}
296		if (LDBL_MANT_DIG-bits >= 2 && !fmodl(frac, 1))
297			frac++;
298	}
299
300	y += frac;
301	y -= bias;
302
303	if ((e2+LDBL_MANT_DIG & INT_MAX) > emax-5) {
304		if (fabs(y) >= CONCAT(0x1p, LDBL_MANT_DIG)) {
305			if (denormal && bits==LDBL_MANT_DIG+e2-emin)
306				denormal = 0;
307			y *= 0.5;
308			e2++;
309		}
310		if (e2+LDBL_MANT_DIG>emax || (denormal && frac))
311			errno = ERANGE;
312	}
313
314	return scalbnl(y, e2);
315}
316
317static long double hexfloat(FILE *f, int bits, int emin, int sign, int pok)
318{
319	uint32_t x = 0;
320	long double y = 0;
321	long double scale = 1;
322	long double bias = 0;
323	int gottail = 0, gotrad = 0, gotdig = 0;
324	long long rp = 0;
325	long long dc = 0;
326	long long e2 = 0;
327	int d;
328	int c;
329
330	c = shgetc(f);
331
332	/* Skip leading zeros */
333	for (; c=='0'; c = shgetc(f)) gotdig = 1;
334
335	if (c=='.') {
336		gotrad = 1;
337		c = shgetc(f);
338		/* Count zeros after the radix point before significand */
339		for (rp=0; c=='0'; c = shgetc(f), rp--) gotdig = 1;
340	}
341
342	for (; c-'0'<10U || (c|32)-'a'<6U || c=='.'; c = shgetc(f)) {
343		if (c=='.') {
344			if (gotrad) break;
345			rp = dc;
346			gotrad = 1;
347		} else {
348			gotdig = 1;
349			if (c > '9') d = (c|32)+10-'a';
350			else d = c-'0';
351			if (dc<8) {
352				x = x*16 + d;
353			} else if (dc < LDBL_MANT_DIG/4+1) {
354				y += d*(scale/=16);
355			} else if (d && !gottail) {
356				y += 0.5*scale;
357				gottail = 1;
358			}
359			dc++;
360		}
361	}
362	if (!gotdig) {
363		shunget(f);
364		if (pok) {
365			shunget(f);
366			if (gotrad) shunget(f);
367		} else {
368			shlim(f, 0);
369		}
370		return sign * 0.0;
371	}
372	if (!gotrad) rp = dc;
373	while (dc<8) x *= 16, dc++;
374	if ((c|32)=='p') {
375		e2 = scanexp(f, pok);
376		if (e2 == LLONG_MIN) {
377			if (pok) {
378				shunget(f);
379			} else {
380				shlim(f, 0);
381				return 0;
382			}
383			e2 = 0;
384		}
385	} else {
386		shunget(f);
387	}
388	e2 += 4*rp - 32;
389
390	if (!x) return sign * 0.0;
391	if (e2 > -emin) {
392		errno = ERANGE;
393		return sign * LDBL_MAX * LDBL_MAX;
394	}
395	if (e2 < emin-2*LDBL_MANT_DIG) {
396		errno = ERANGE;
397		return sign * LDBL_MIN * LDBL_MIN;
398	}
399
400	while (x < 0x80000000) {
401		if (y>=0.5) {
402			x += x + 1;
403			y += y - 1;
404		} else {
405			x += x;
406			y += y;
407		}
408		e2--;
409	}
410
411	if (bits > 32+e2-emin) {
412		bits = 32+e2-emin;
413		if (bits<0) bits=0;
414	}
415
416	if (bits < LDBL_MANT_DIG)
417		bias = copysignl(scalbn(1, 32+LDBL_MANT_DIG-bits-1), sign);
418
419	if (bits<32 && y && !(x&1)) x++, y=0;
420
421	y = bias + sign*(long double)x + sign*y;
422	y -= bias;
423
424	if (!y) errno = ERANGE;
425
426	return scalbnl(y, e2);
427}
428
429long double __floatscan(FILE *f, int prec, int pok)
430{
431	int sign = 1;
432	size_t i;
433	int bits;
434	int emin;
435	int c;
436
437	switch (prec) {
438	case 0:
439		bits = FLT_MANT_DIG;
440		emin = FLT_MIN_EXP-bits;
441		break;
442	case 1:
443		bits = DBL_MANT_DIG;
444		emin = DBL_MIN_EXP-bits;
445		break;
446	case 2:
447		bits = LDBL_MANT_DIG;
448		emin = LDBL_MIN_EXP-bits;
449		break;
450	default:
451		return 0;
452	}
453
454	while (isspace((c=shgetc(f))));
455
456	if (c=='+' || c=='-') {
457		sign -= 2*(c=='-');
458		c = shgetc(f);
459	}
460
461	for (i=0; i<8 && (c|32)=="infinity"[i]; i++)
462		if (i<7) c = shgetc(f);
463	if (i==3 || i==8 || (i>3 && pok)) {
464		if (i!=8) {
465			shunget(f);
466			if (pok) for (; i>3; i--) shunget(f);
467		}
468		return sign * INFINITY;
469	}
470	if (!i) for (i=0; i<3 && (c|32)=="nan"[i]; i++)
471		if (i<2) c = shgetc(f);
472	if (i==3) {
473		if (shgetc(f) != '(') {
474			shunget(f);
475			return NAN;
476		}
477		for (i=1; ; i++) {
478			c = shgetc(f);
479			if (c-'0'<10U || c-'A'<26U || c-'a'<26U || c=='_')
480				continue;
481			if (c==')') return NAN;
482			shunget(f);
483			if (!pok) {
484				errno = EINVAL;
485				shlim(f, 0);
486				return 0;
487			}
488			while (i--) shunget(f);
489			return NAN;
490		}
491		return NAN;
492	}
493
494	if (i) {
495		shunget(f);
496		errno = EINVAL;
497		shlim(f, 0);
498		return 0;
499	}
500
501	if (c=='0') {
502		c = shgetc(f);
503		if ((c|32) == 'x')
504			return hexfloat(f, bits, emin, sign, pok);
505		shunget(f);
506		c = '0';
507	}
508
509	return decfloat(f, c, bits, emin, sign, pok);
510}
511