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 
scanexp(FILE * f,int pok)39 static 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 
decfloat(FILE * f,int c,int bits,int emin,int sign,int pok)66 static 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 
hexfloat(FILE * f,int bits,int emin,int sign,int pok)317 static 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 
__floatscan(FILE * f,int prec,int pok)429 long 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