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