#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, -1, -1 #define KMAX 128 #else /* LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 */ #define LD_B1B_DIG 3 #define LD_B1B_MAX 18, 446744073, 709551615, -1 #define KMAX 2048 #endif #define MASK (KMAX-1) int dohex(FILE *f, long double *y, int sign) { return -1; } int parsefloat(FILE *f, long double *yout) { uint32_t x[KMAX]; static const uint32_t th[] = { LD_B1B_MAX }; int c; int i, j, k, a, z; int sign=1; int rp=-1, dc=0; int e10; int e2; long double y; long double frac=0; long double bias=0; int bits=64; int emin=-16445; j=0; k=0; c = getc(f); if (c=='+' || c=='-') { sign -= 2*(c=='-'); c = getc(f); } if (c=='0') { j = 1; c = getc(f); if ((c|32) == 'x') { return dohex(f, yout, sign); } } /* Don't let leading zeros consume buffer space */ for (; c=='0'; c=getc(f)); for (; c-'0'<10U || c=='.'; c=getc(f)) { if (c == '.') { if (rp!=-1) goto mismatch; rp = dc; } else if (k < KMAX) { dc++; x[k] = x[k]*10 + c-'0'; if (++j==9) { k++; j=0; } } else { dc++; x[KMAX-1] |= c-'0'; } } if (rp==-1) rp=dc; if (k=10U) { goto mismatch; } fscanf(f, "%d", &e10); c = getc(f); rp += e10; /* FIXME: overflow */ } a = 0; z = k; e2 = 0; while (rp < 18+9*LD_B1B_DIG) { uint32_t carry = 0; e2 -= 29; for (k=(z-1 & MASK); ; k=(k-1 & MASK)) { uint64_t tmp = ((uint64_t)x[k] << 29) + carry; if (tmp > 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) z = (z+1 & MASK); x[z-1 & MASK] = carry; } 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) z = (z+1 & MASK); x[z-1 & MASK] = carry; } } for (i=a; i!=z; i=(i+1 & MASK)) printf("%.9u ", x[i]); printf("[%d]\n%*s\n", e2, rp+(rp/9)+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) { printf("old frac = %La\n", frac); 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; } } printf("y = %.21Lg\n", scalbnl(y-bias,e2)); printf("frac = %.21Lg\n", scalbnl(frac,e2)); printf("y = %La\n", y); printf("frac = %La\n", frac); printf("bias = %La\n", bias); y += frac; y -= bias; y = scalbnl(y, e2); printf("y = %La\n", y); printf("y ~= %.17Lg\n", y); printf("%.21Lg [%.21Lg, %.21Lg]\n", y, nextafterl(y,-INFINITY), nextafterl(y, INFINITY)); #if 0 printf("%.18g [%.18g, %.18g]\n", (double)y, nextafter(y,-INFINITY), nextafter(y, INFINITY)); printf("%.9g [%.9g, %.9g]\n", (float)y, nextafterf(y,-INFINITY), nextafterf(y, INFINITY)); #endif return 0; mismatch: return -1; } int main() { long double y; //fesetround(FE_DOWNWARD); parsefloat(stdin, &y); }