#define _XOPEN_SOURCE 700 #include #include #include #include #include #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 #define LD_B1B_DIG 2 #define LD_B1B_MAX 9007199, 254740991 #define KMAX 128 #else /* LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 */ #define LD_B1B_DIG 3 #define LD_B1B_MAX 18, 446744073, 709551615 #define KMAX 2048 #endif #define MASK (KMAX-1) #if 0 static inline fast_unget(int c, FILE *f) { if (f->rpos) f->rpos--; } #undef ungetc #define ungetc fast_unget #endif #undef getc #define getc getc_unlocked static long double decfloat(FILE *f, int c, int prec, int sign, int pok, off_t *pcnt) { uint32_t x[KMAX]; static const uint32_t th[] = { LD_B1B_MAX }; int i, j, k, a, z; int rp=-1, dc=0; int e10=0; int e2; long double y; long double frac=0; long double bias=0; int bits=64; int emin; switch (prec) { case 0: bits = 24; emin = -149; break; case 1: bits = 53; emin = -1074; break; case 2: bits = LDBL_MANT_DIG; emin = -16445; break; default: return 0; } 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); x[0] = 0; for (; c-'0'<10U || c=='.'; *pcnt += (c=getc(f))>=0) { if (c == '.') { if (rp!=-1) break; rp = dc; } else if (k < KMAX) { dc++; if (j) x[k] = x[k]*10 + c-'0'; else x[k] = c-'0'; if (++j==9) { k++; j=0; } } else { dc++; x[KMAX-1] |= c-'0'; } } if (rp==-1) rp=dc; if ((c|32)=='e') { c = getc(f); ungetc(c, f); if (c!='+' && c!='-' && c-'0'>=10U) { if (!pok) { *pcnt = 0; return 0; } else { --*pcnt; } } else { int tmp = 0; fscanf(f, "%d%n", &e10, &tmp); *pcnt += tmp; rp += e10; /* FIXME: overflow */ } } else if (c>=0) { ungetc(c, f); --*pcnt; } if (rp==dc && (!k || (k==1 && !j)) && (bits>30 || x[0]>>bits==0)) { y = sign * (long double)x[0]; return y; } if (k 1000000000) { carry = tmp / 1000000000; x[k] = tmp % 1000000000; } else { carry = 0; x[k] = tmp; } if (k==(z-1 & MASK) && k!=a && !x[k]) z = k; if (k==a) break; } if (carry) { rp += 9; if (a == z) { z = (z-1 & MASK); x[z-1 & MASK] |= x[z]; } a = (a-1 & MASK); x[a] = carry; } } 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; for (i=0; i th[i]) break; } if (i==LD_B1B_DIG && rp==9*LD_B1B_DIG) break; /* 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)) { uint32_t tmp = x[k] & (1<>sh) + carry; carry = (1000000000>>sh) * 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; } } for (y=i=0; i LDBL_MANT_DIG+e2-emin) { bits = LDBL_MANT_DIG+e2-emin; if (bits<0) bits=0; } if (bits < LDBL_MANT_DIG) { bias = copysignl(scalbn(1, 2*LDBL_MANT_DIG-bits-1), y); frac = fmodl(y, scalbn(1, LDBL_MANT_DIG-bits)); y -= frac; y += bias; } if ((a+i & MASK) != z) { uint32_t t = x[a+i & MASK]; if (t < 500000000 && (t || (a+i+1 & MASK) != z)) frac += 0.25*sign; else if (t > 500000000) frac += 0.75*sign; else if (t == 500000000) { if ((a+i+1 & MASK) == z) frac += 0.5*sign; else frac += 0.75*sign; } } y += frac; y -= bias; y = scalbnl(y, e2); return y; } static long double hexfloat(FILE *f, int c, int prec, int sign, int pok, off_t *pcnt) { return 0; } long double __floatscan(FILE *f, int c, int prec, int pok, off_t *pcnt) { int sign = 1; int i; if (c<0) *pcnt += (c=getc(f))>=0; if (c=='+' || c=='-') { sign -= 2*(c=='-'); *pcnt += (c=getc(f))>=0; } for (i=0; i<8 && (c|32)=="infinity"[i]; i++) if (i<7) c = getc(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; return sign * INFINITY; } if (!i) for (i=0; i<3 && (c|32)=="nan"[i]; i++) if (i<3) c = getc(f); if (i==3) { *pcnt += 2; return sign>0 ? NAN : -NAN; } if (i || c-'0'>10U) { ungetc(c, f); *pcnt = 0; return -1; } if (c=='0') { *pcnt += (c=getc(f))>=0; if ((c|32) == 'x') return hexfloat(f, -1, prec, sign, pok, pcnt); } return decfloat(f, c, prec, sign, pok, pcnt); } #include int main(int argc, char **argv) { FILE *f = 0; if (argc>1) f = fmemopen(argv[1], strlen(argv[1]), "r"); if (!f) f = stdin; off_t pos = 0; long double y = __floatscan(f, -1, 0, 1, &pos); printf("%.21Lg (read %lld)\n", y, (long long)pos); return 0; }