diff --git a/c-programming/floatscan_experiment.c b/c-programming/floatscan_experiment.c new file mode 100644 index 0000000..8daa280 --- /dev/null +++ b/c-programming/floatscan_experiment.c @@ -0,0 +1,533 @@ +/* + * This was an experiment to learn about how libc turn floating point values + * into actual values. + * + * Expects the floating point string to be in stdin. + * + * The original code is from musl libc: version 1.2.4, src/internal/floatscan.c + */ + +#include +#include +#include +#include +#include +#include +#include + +/* wrapping code --- beginning */ +#define shgetc(f) fgetc(f) +#define shlim(...) +#define shunget(...) + +#define FLOAT_TYPE 0 +#define DOUBLE_TYPE 1 +#define LONG_DOUBLE_TYPE 2 +/* wrapping code --- end */ + +/* the original musl libc code --- beginning */ +#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 + +#define LD_B1B_DIG 2 +#define LD_B1B_MAX 9007199, 254740991 +#define KMAX 128 + +#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 + +#define LD_B1B_DIG 3 +#define LD_B1B_MAX 18, 446744073, 709551615 +#define KMAX 2048 + +#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 + +#define LD_B1B_DIG 4 +#define LD_B1B_MAX 10384593, 717069655, 257060992, 658440191 +#define KMAX 2048 + +#else +#error Unsupported long double representation +#endif + +#define MASK (KMAX-1) + +static long long scanexp(FILE *f, int pok) +{ + int c; + int x; + long long y; + int neg = 0; + + c = shgetc(f); + if (c=='+' || c=='-') { + neg = (c=='-'); + c = shgetc(f); + if (c-'0'>=10U && pok) shunget(f); + } + if (c-'0'>=10U) { + shunget(f); + return LLONG_MIN; + } + for (x=0; c-'0'<10U && x=0) { + shunget(f); + } + if (!gotdig) { + errno = EINVAL; + shlim(f, 0); + return 0; + } + + /* Handle zero specially to avoid nasty special cases later */ + if (!x[0]) return sign * 0.0; + + /* Optimize small integers (w/no exponent) and over/under-flow */ + if (lrp==dc && dc<10 && (bits>30 || x[0]>>bits==0)) + return sign * (long double)x[0]; + if (lrp > -emin/2) { + errno = ERANGE; + return sign * LDBL_MAX * LDBL_MAX; + } + if (lrp < emin-2*LDBL_MANT_DIG) { + errno = ERANGE; + return sign * LDBL_MIN * LDBL_MIN; + } + + /* Align incomplete final B1B digit */ + if (j) { + for (; j<9; j++) x[k]*=10; + k++; + j=0; + } + + a = 0; + z = k; + e2 = 0; + rp = lrp; + + /* Optimize small to mid-size integers (even in exp. notation) */ + if (lnz<9 && lnz<=rp && rp < 18) { + if (rp == 9) return sign * (long double)x[0]; + if (rp < 9) return sign * (long double)x[0] / p10s[8-rp]; + int bitlim = bits-3*(int)(rp-9); + if (bitlim>30 || x[0]>>bitlim==0) + return sign * (long double)x[0] * p10s[rp-10]; + } + + /* Drop trailing zeros */ + for (; !x[z-1]; z--); + + /* Align radix point to B1B digit boundary */ + if (rp % 9) { + int rpm9 = rp>=0 ? rp%9 : rp%9+9; + int p10 = p10s[8-rpm9]; + 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; + } + + /* Upscale until desired number of bits are left of radix point */ + while (rp < 9*LD_B1B_DIG || (rp == 9*LD_B1B_DIG && x[a] 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; + a = (a-1 & MASK); + if (a == z) { + z = (z-1 & MASK); + x[z-1 & MASK] |= x[z]; + } + x[a] = carry; + } + } + + /* Downscale until exactly number of bits are left of radix point */ + 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); + i--; + rp -= 9; + } + } + if (carry) { + if ((z+1 & MASK) != a) { + x[z] = carry; + z = (z+1 & MASK); + } else x[z-1 & MASK] |= 1; + } + } + + /* Assemble desired bits into floating point variable */ + for (y=i=0; i LDBL_MANT_DIG+e2-emin) { + bits = LDBL_MANT_DIG+e2-emin; + if (bits<0) bits=0; + denormal = 1; + } + + /* Calculate bias term to force rounding, move out lower bits */ + 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; + } + + /* Process tail of decimal input so it can affect rounding */ + 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; + } + if (LDBL_MANT_DIG-bits >= 2 && !fmodl(frac, 1)) + frac++; + } + + y += frac; + y -= bias; + + if ((e2+LDBL_MANT_DIG & INT_MAX) > emax-5) { + if (fabsl(y) >= 2/LDBL_EPSILON) { + if (denormal && bits==LDBL_MANT_DIG+e2-emin) + denormal = 0; + y *= 0.5; + e2++; + } + if (e2+LDBL_MANT_DIG>emax || (denormal && frac)) + errno = ERANGE; + } + + return scalbnl(y, e2); +} + +static long double hexfloat(FILE *f, int bits, int emin, int sign, int pok) +{ + uint32_t x = 0; + long double y = 0; + long double scale = 1; + long double bias = 0; + int gottail = 0, gotrad = 0, gotdig = 0; + long long rp = 0; + long long dc = 0; + long long e2 = 0; + int d; + int c; + + c = shgetc(f); + + /* Skip leading zeros */ + for (; c=='0'; c = shgetc(f)) gotdig = 1; + + if (c=='.') { + gotrad = 1; + c = shgetc(f); + /* Count zeros after the radix point before significand */ + for (rp=0; c=='0'; c = shgetc(f), rp--) gotdig = 1; + } + + for (; c-'0'<10U || (c|32)-'a'<6U || c=='.'; c = shgetc(f)) { + if (c=='.') { + if (gotrad) break; + rp = dc; + gotrad = 1; + } else { + gotdig = 1; + if (c > '9') d = (c|32)+10-'a'; + else d = c-'0'; + if (dc<8) { + x = x*16 + d; + } else if (dc < LDBL_MANT_DIG/4+1) { + y += d*(scale/=16); + } else if (d && !gottail) { + y += 0.5*scale; + gottail = 1; + } + dc++; + } + } + if (!gotdig) { + shunget(f); + if (pok) { + shunget(f); + if (gotrad) shunget(f); + } else { + shlim(f, 0); + } + return sign * 0.0; + } + if (!gotrad) rp = dc; + while (dc<8) x *= 16, dc++; + if ((c|32)=='p') { + e2 = scanexp(f, pok); + if (e2 == LLONG_MIN) { + 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) { + 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) { + x += x + 1; + y += y - 1; + } else { + x += x; + y += y; + } + e2--; + } + + if (bits > 32+e2-emin) { + bits = 32+e2-emin; + if (bits<0) bits=0; + } + + if (bits < LDBL_MANT_DIG) + bias = copysignl(scalbn(1, 32+LDBL_MANT_DIG-bits-1), sign); + + if (bits<32 && y && !(x&1)) x++, y=0; + + y = bias + sign*(long double)x + sign*y; + y -= bias; + + if (!y) errno = ERANGE; + + return scalbnl(y, e2); +} + +long double __floatscan(FILE *f, int prec, int pok) +{ + int sign = 1; + size_t i; + int bits; + int emin; + int c; + + switch (prec) { + case 0: + bits = FLT_MANT_DIG; + emin = FLT_MIN_EXP-bits; + break; + case 1: + bits = DBL_MANT_DIG; + emin = DBL_MIN_EXP-bits; + break; + case 2: + bits = LDBL_MANT_DIG; + emin = LDBL_MIN_EXP-bits; + break; + default: + return 0; + } + + while (isspace((c=shgetc(f)))); + + if (c=='+' || c=='-') { + sign -= 2*(c=='-'); + c = shgetc(f); + } + + for (i=0; i<8 && (c|32)=="infinity"[i]; i++) + if (i<7) c = shgetc(f); + if (i==3 || i==8 || (i>3 && pok)) { + if (i!=8) { + shunget(f); + if (pok) for (; i>3; i--) shunget(f); + } + return sign * INFINITY; + } + if (!i) for (i=0; i<3 && (c|32)=="nan"[i]; i++) + if (i<2) c = shgetc(f); + if (i==3) { + if (shgetc(f) != '(') { + shunget(f); + return NAN; + } + for (i=1; ; i++) { + c = shgetc(f); + if (c-'0'<10U || c-'A'<26U || c-'a'<26U || c=='_') + continue; + if (c==')') return NAN; + shunget(f); + if (!pok) { + errno = EINVAL; + shlim(f, 0); + return 0; + } + while (i--) shunget(f); + return NAN; + } + return NAN; + } + + if (i) { + shunget(f); + errno = EINVAL; + shlim(f, 0); + return 0; + } + + if (c=='0') { + c = shgetc(f); + if ((c|32) == 'x') + return hexfloat(f, bits, emin, sign, pok); + shunget(f); + c = '0'; + } + + return decfloat(f, c, bits, emin, sign, pok); +} +/* the original musl libc code --- end */ + +/* test code --- beginning */ +int main(void) { + printf("LD_B1B_DIG = %d\n", LD_B1B_DIG); + printf("%La\n", __floatscan(stdin, LONG_DOUBLE_TYPE, 1)); + return 0; +} +/* test code --- end */