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