optimize floatscan downscaler to skip results that won't be needed
[musl] / src / internal / floatscan.c
index 15ad5e1..b231329 100644 (file)
@@ -3,7 +3,9 @@
 #include <math.h>
 #include <float.h>
 #include <limits.h>
+#include <errno.h>
 
+#include "shgetc.h"
 #include "floatscan.h"
 
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
 #define MASK (KMAX-1)
 
 
-#if 1
-#include "stdio_impl.h"
-#undef ungetc
-#define ungetc(c,f) ((f)->rpos--,(c))
-#undef getc
-#define getc getc_unlocked
-#endif
-
-
-static long long scanexp(FILE *f, off_t *pcnt)
+static long long scanexp(FILE *f, int pok)
 {
        int c;
        int x;
        long long y;
        int neg = 0;
        
-       *pcnt += (c=getc(f))>=0;
+       c = shgetc(f);
        if (c=='+' || c=='-') {
                neg = (c=='-');
-               *pcnt += (c=getc(f))>=0;
-               if (c-'0'>=10U) {
-                       if (c>=0) {
-                               ungetc(c, f);
-                               --*pcnt;
-                       }
-                       return LLONG_MIN;
-               }
+               c = shgetc(f);
+               if (c-'0'>=10U && pok) shunget(f);
        }
-       for (x=0; c-'0'<10U && x<INT_MAX/10; *pcnt += (c=getc(f))>=0)
+       if (c-'0'>=10U) {
+               shunget(f);
+               return LLONG_MIN;
+       }
+       for (x=0; c-'0'<10U && x<INT_MAX/10; c = shgetc(f))
                x = 10*x + c-'0';
-       for (y=x; c-'0'<10U && x<LLONG_MAX/10; *pcnt += (c=getc(f))>=0)
+       for (y=x; c-'0'<10U && x<LLONG_MAX/100; c = shgetc(f))
                y = 10*y + c-'0';
-       for (; c-'0'<10U; *pcnt += (c=getc(f))>=0);
-       if (c>=0) {
-               ungetc(c, f);
-               --*pcnt;
-       }
+       for (; c-'0'<10U; c = shgetc(f));
+       shunget(f);
        return neg ? -y : y;
 }
 
 
-static long double decfloat(FILE *f, int c, int bits, int emin, int sign, int pok, off_t *pcnt)
+static long double decfloat(FILE *f, int c, int bits, int emin, int sign, int pok)
 {
        uint32_t x[KMAX];
        static const uint32_t th[] = { LD_B1B_MAX };
        int i, j, k, a, z;
        long long lrp=-1, dc=0;
+       long long e10=0;
        int gotdig = 0;
        int rp;
-       int e10=0;
        int e2;
        long double y;
        long double frac=0;
@@ -81,17 +69,15 @@ static long double decfloat(FILE *f, int c, int bits, int emin, int sign, int po
        j=0;
        k=0;
 
-       if (c<0) *pcnt += (c=getc(f))>=0;
-
        /* Don't let leading zeros consume buffer space */
-       for (; c=='0'; *pcnt += (c=getc(f))>=0) gotdig=1;
+       for (; c=='0'; c = shgetc(f)) gotdig=1;
 
        x[0] = 0;
-       for (; c-'0'<10U || c=='.'; *pcnt += (c=getc(f))>=0) {
+       for (; c-'0'<10U || c=='.'; c = shgetc(f)) {
                if (c == '.') {
                        if (lrp!=-1) break;
                        lrp = dc;
-               } else if (k < KMAX) {
+               } else if (k < KMAX-2) {
                        dc++;
                        if (j) x[k] = x[k]*10 + c-'0';
                        else x[k] = c-'0';
@@ -102,27 +88,29 @@ static long double decfloat(FILE *f, int c, int bits, int emin, int sign, int po
                        gotdig=1;
                } else {
                        dc++;
-                       x[KMAX-1] |= c-'0';
+                       if (c!='0') x[KMAX-3] |= 1;
                }
        }
        if (lrp==-1) lrp=dc;
 
        if (gotdig && (c|32)=='e') {
-               e10 = scanexp(f, pcnt);
+               e10 = scanexp(f, pok);
                if (e10 == LLONG_MIN) {
-                       if (!pok) {
-                               *pcnt = 0;
+                       if (pok) {
+                               shunget(f);
+                       } else {
+                               shlim(f, 0);
                                return 0;
                        }
                        e10 = 0;
                }
                lrp += e10;
        } else if (c>=0) {
-               ungetc(c, f);
-               --*pcnt;
+               shunget(f);
        }
        if (!gotdig) {
-               *pcnt = 0;
+               errno = EINVAL;
+               shlim(f, 0);
                return 0;
        }
 
@@ -130,10 +118,14 @@ static long double decfloat(FILE *f, int c, int bits, int emin, int sign, int po
                return sign * 0.0;
        if (lrp==dc && (!k || (k==1 && !j)) && (bits>30 || x[0]>>bits==0))
                return sign * (long double)x[0];
-       if (lrp > -emin/2)
+       if (lrp > -emin/2) {
+               errno = ERANGE;
                return sign * LDBL_MAX * LDBL_MAX;
-       if (lrp < emin-2*LDBL_MANT_DIG)
+       }
+       if (lrp < emin-2*LDBL_MANT_DIG) {
+               errno = ERANGE;
                return sign * LDBL_MIN * LDBL_MIN;
+       }
 
        if (k<KMAX && j) {
                for (; j<9; j++) x[k]*=10;
@@ -146,7 +138,28 @@ static long double decfloat(FILE *f, int c, int bits, int emin, int sign, int po
        e2 = 0;
        rp = lrp;
 
-       while (rp < 18+9*LD_B1B_DIG) {
+       if (rp % 9) {
+               static const int p10s[] = {
+                       100000000, 10000000, 1000000, 100000,
+                       10000, 1000, 100, 10
+               };
+               int rpm9 = rp>=0 ? rp%9 : rp%9+9;
+               int p10 = p10s[rpm9-1];
+               uint32_t carry = 0;
+               for (k=a; k!=z; k++) {
+                       uint32_t tmp = x[k] % p10;
+                       x[k] = x[k]/p10 + carry;
+                       carry = 1000000000/p10 * tmp;
+                       if (k==a && !x[k]) {
+                               a = (a+1 & MASK);
+                               rp -= 9;
+                       }
+               }
+               if (carry) x[z++] = carry;
+               rp += 9-rpm9;
+       }
+
+       while (rp < 9*LD_B1B_DIG || (rp == 9*LD_B1B_DIG && x[0]<th[0])) {
                uint32_t carry = 0;
                e2 -= 29;
                for (k=(z-1 & MASK); ; k=(k-1 & MASK)) {
@@ -172,32 +185,6 @@ static long double decfloat(FILE *f, int c, int bits, int emin, int sign, int po
                }
        }
 
-       if (rp % 9) {
-               static const int p10s[] = {
-                       100000000, 10000000, 1000000, 100000,
-                       10000, 1000, 100, 10
-               };
-               int rpm9 = rp % 9;
-               int p10 = p10s[rpm9-1];
-               uint32_t carry = 0;
-               for (k=a; k!=z; k=(k+1 & MASK)) {
-                       uint32_t tmp = x[k] % p10;
-                       x[k] = x[k]/p10 + carry;
-                       carry = 1000000000/p10 * tmp;
-                       if (k==a && !x[k]) {
-                               a = (a+1 & MASK);
-                               rp -= 9;
-                       }
-               }
-               if (carry) {
-                       if ((z+1 & MASK) != a) {
-                               x[z] = carry;
-                               z = (z+1 & MASK);
-                       } else x[z-1 & MASK] |= 1;
-               }
-               rp += 9-rpm9;
-       }
-
        for (;;) {
                uint32_t carry = 0;
                int sh = 1;
@@ -213,16 +200,17 @@ static long double decfloat(FILE *f, int c, int bits, int emin, int sign, int po
                /* FIXME: find a way to compute optimal sh */
                if (rp > 9+9*LD_B1B_DIG) sh = 9;
                e2 += sh;
-               for (k=a; k!=z; k=(k+1 & MASK)) {
+               for (i=0; (k=(a+i & MASK))!=z && i<LD_B1B_DIG+3; i++) {
                        uint32_t tmp = x[k] & (1<<sh)-1;
                        x[k] = (x[k]>>sh) + carry;
                        carry = (1000000000>>sh) * tmp;
                        if (k==a && !x[k]) {
                                a = (a+1 & MASK);
+                               i--;
                                rp -= 9;
                        }
                }
-               if (carry) {
+               if (carry && k==z) {
                        if ((z+1 & MASK) != a) {
                                x[z] = carry;
                                z = (z+1 & MASK);
@@ -230,8 +218,10 @@ static long double decfloat(FILE *f, int c, int bits, int emin, int sign, int po
                }
        }
 
-       for (y=i=0; i<LD_B1B_DIG && (a+i & MASK)!=z; i++)
+       for (y=i=0; i<LD_B1B_DIG; i++) {
+               if ((a+i & MASK)==z) x[z=(z+1 & MASK)] = 0;
                y = 1000000000.0L * y + x[a+i & MASK];
+       }
 
        y *= sign;
 
@@ -268,10 +258,12 @@ static long double decfloat(FILE *f, int c, int bits, int emin, int sign, int po
 
        y = scalbnl(y, e2);
 
+       if (!y) errno = ERANGE;
+
        return y;
 }
 
-static long double hexfloat(FILE *f, int c, int bits, int emin, int sign, int pok, off_t *pcnt)
+static long double hexfloat(FILE *f, int bits, int emin, int sign, int pok)
 {
        uint32_t x = 0;
        long double y = 0;
@@ -282,20 +274,21 @@ static long double hexfloat(FILE *f, int c, int bits, int emin, int sign, int po
        long long dc = 0;
        long long e2 = 0;
        int d;
+       int c;
 
-       if (c<0) *pcnt += (c=getc(f))>=0;
+       c = shgetc(f);
 
        /* Skip leading zeros */
-       for (; c=='0'; *pcnt += (c=getc(f))>=0) gotdig = 1;
+       for (; c=='0'; c = shgetc(f)) gotdig = 1;
 
        if (c=='.') {
                gotrad = 1;
-               *pcnt += (c=getc(f))>=0;
+               c = shgetc(f);
                /* Count zeros after the radix point before significand */
-               for (rp=0; c=='0'; *pcnt += (c=getc(f))>=0, rp--) gotdig = 1;
+               for (rp=0; c=='0'; c = shgetc(f), rp--) gotdig = 1;
        }
 
-       for (; c-'0'<10U || (c|32)-'a'<6U || c=='.'; *pcnt += (c=getc(f))>=0) {
+       for (; c-'0'<10U || (c|32)-'a'<6U || c=='.'; c = shgetc(f)) {
                if (c=='.') {
                        if (gotrad) break;
                        rp = dc;
@@ -316,31 +309,42 @@ static long double hexfloat(FILE *f, int c, int bits, int emin, int sign, int po
                }
        }
        if (!gotdig) {
-               if (c>=0) {
-                       ungetc(c, f);
-                       --*pcnt;
+               shunget(f);
+               if (pok) {
+                       shunget(f);
+                       if (gotrad) shunget(f);
+               } else {
+                       shlim(f, 0);
                }
-               if (pok) *pcnt -= 1+gotrad; /* uncount the rp, x of 0x */
-               else *pcnt = 0;
                return 0;
        }
        if (!gotrad) rp = dc;
        while (dc<8) x *= 16, dc++;
        if ((c|32)=='p') {
-               e2 = scanexp(f, pcnt);
+               e2 = scanexp(f, pok);
                if (e2 == LLONG_MIN) {
-                       if (!pok) {
-                               *pcnt = 0;
+                       if (pok) {
+                               shunget(f);
+                       } else {
+                               shlim(f, 0);
                                return 0;
                        }
                        e2 = 0;
                }
+       } else {
+               shunget(f);
        }
        e2 += 4*rp - 32;
 
        if (!x) return sign * 0.0;
-       if (e2 > -emin) return sign * LDBL_MAX * LDBL_MAX;
-       if (e2 < emin-2*LDBL_MANT_DIG) return sign * LDBL_MIN * LDBL_MIN;
+       if (e2 > -emin) {
+               errno = ERANGE;
+               return sign * LDBL_MAX * LDBL_MAX;
+       }
+       if (e2 < emin-2*LDBL_MANT_DIG) {
+               errno = ERANGE;
+               return sign * LDBL_MIN * LDBL_MIN;
+       }
 
        while (x < 0x80000000) {
                if (y>=0.5) {
@@ -366,18 +370,18 @@ static long double hexfloat(FILE *f, int c, int bits, int emin, int sign, int po
        y = bias + sign*(long double)x + sign*y;
        y -= bias;
 
+       if (!y) errno = ERANGE;
+
        return scalbnl(y, e2);
 }
 
-long double __floatscan(FILE *f, int c, int prec, int pok, off_t *pcnt)
+long double __floatscan(FILE *f, int c, int prec, int pok)
 {
        int sign = 1;
        int i;
        int bits;
        int emin;
 
-       *pcnt = 0;
-
        switch (prec) {
        case 0:
                bits = 24;
@@ -395,44 +399,41 @@ long double __floatscan(FILE *f, int c, int prec, int pok, off_t *pcnt)
                return 0;
        }
 
-       if (c<0) *pcnt += (c=getc(f))>=0;
+       if (c<0) c = shgetc(f);
 
        if (c=='+' || c=='-') {
                sign -= 2*(c=='-');
-               *pcnt += (c=getc(f))>=0;
+               c = shgetc(f);
        }
 
        for (i=0; i<8 && (c|32)=="infinity"[i]; i++)
-               if (i<7) c = getc(f);
+               if (i<7) c = shgetc(f);
        if (i==3 || i==8 || (i>3 && pok)) {
-               if (i==3 && c>=0) ungetc(c, f);
-               if (i==8) *pcnt += 7;
-               else *pcnt += 2;
+               if (i==3) shunget(f);
+               if (pok) for (; i>3; i--) shunget(f);
+               else shlim(f, 0);
                return sign * INFINITY;
        }
        if (!i) for (i=0; i<3 && (c|32)=="nan"[i]; i++)
-               if (i<3) c = getc(f);
+               if (i<3) c = shgetc(f);
        if (i==3) {
-               *pcnt += 2;
-               return sign>0 ? NAN : -NAN;
+               return NAN;
        }
 
        if (i) {
-               if (c>=0) ungetc(c, f);
-               *pcnt = 0;
+               shunget(f);
+               errno = EINVAL;
+               shlim(f, 0);
                return 0;
        }
 
        if (c=='0') {
-               *pcnt += (c=getc(f))>=0;
+               c = shgetc(f);
                if ((c|32) == 'x')
-                       return hexfloat(f, -1, bits, emin, sign, pok, pcnt);
-               if (c>=0) {
-                       ungetc(c, f);
-                       --*pcnt;
-               }
+                       return hexfloat(f, bits, emin, sign, pok);
+               shunget(f);
                c = '0';
        }
 
-       return decfloat(f, c, bits, emin, sign, pok, pcnt);
+       return decfloat(f, c, bits, emin, sign, pok);
 }