/* * fpconv.c * * Convert a decimal floating-point number to IEEE with correct * rounding in all cases. * * This code is written for correctness and portability, not * speed. By leveraging the FPU on the host computer you can often * get much faster results. * * The basic algorithm involves treating the decimal floating-point * number as a fraction multiplied by a power of two: p/q * 2^n * where q is always a power of 5. */ #include #include #include #include #include #include #include #include "pbn.h" #ifdef __GNUC__ #define likely(x) __builtin_expect(!!(x), 1) #define unlikely(x) __builtin_expect(!!(x), 0) #else #define likely(x) (!!(x)) #define unlikely(x) (!!(x)) #endif /* XXX: consider caching powers of five, or at least 5^2^n */ static struct pbn *power_of_five(unsigned int pw) { static struct pbn **powers = NULL; static unsigned int max_power = 0; struct pbn *p; unsigned int pwx, lowbit; if (unlikely(pw >= max_power)) { struct pbn **npowers; pwx = (pw+1) << 1; npowers = realloc(powers, pwx * sizeof(struct pbn *)); if (!npowers) return NULL; memset(npowers + max_power, 0, (pwx - max_power) * sizeof(struct pbn *)); if (!powers) { /* The first time through... */ npowers[0] = pbn_int(1); npowers[1] = pbn_int(5); } powers = npowers; max_power = pwx; } else if (likely(powers[pw])) { return pbn_ref(powers[pw]); } lowbit = pw & -pw; if (lowbit == pw) { /* A power of two */ struct pbn *t = power_of_five(pw >> 1); p = pbn_mul(t, pbn_ref(t)); } else { /* Not a power of two */ p = pbn_mul(power_of_five(lowbit), power_of_five(pw & ~lowbit)); } printf("5^%-4u : ", pw); pbn_dump(stdout, p); putchar('\n'); powers[pw] = pbn_ref(p); return p; } struct fpgen { struct pbn *mnt; /* Mantissa */ struct pbn *num; /* Remainder numerator */ struct pbn *den; /* Remainder denominator */ unsigned int pw2; /* Power of two */ }; enum fp_type { FP_INVALID, FP_ZERO, FP_INT, FP_FRAC, }; #if PBN_LIMB_BITS == 32 #define LIMB_DIGITS 9 #define LIMB_DIGITS_VAL UINT32_C(1000000000) #elif PBN_LIMB_BITS == 64 #define LIMB_DIGITS 19 #define LIMB_DIGITS_VAL UINT64_C(10000000000000000000) #else #error "Need to define LIMB_DIGITS for PBN_LIMB_BITS" #endif static enum fp_type analyze_string(const char *str, struct fpgen *fpg) { const char *p, *dot, *first_dig, *last_dig; char c; long exponent, exponent2, exponent5; int tz; int since_dot, digits, zeroes; int fl; pbn_limb_t lv; struct pbn *num, *den; memset(fpg, 0, sizeof *fpg); dot = first_dig = last_dig = NULL; exponent = digits = zeroes = since_dot = 0; for (p = str; (c = *p); p++) { switch (c) { case 'e': case 'E': exponent += strtol(p+1, NULL, 10); goto done; case '.': if (!dot) dot = p; else return FP_INVALID; break; case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9': if (!first_dig) { first_dig = p; zeroes = 0; } last_dig = p; digits += zeroes + 1; zeroes = 0; exponent = dot ? --since_dot : 0; break; case '0': if (dot) --since_dot; else exponent++; zeroes++; break; case '_': break; default: return FP_INVALID; } } done: if (unlikely(!first_dig)) { printf("Zero mantissa\n"); return FP_ZERO; } fl = ((digits-1) % LIMB_DIGITS) + 1; lv = 0; num = NULL; for (p = first_dig; p <= last_dig; p++) { c = *p; putchar(c); if (c >= '0' && c <= '9') { lv = (lv * 10) + (c - '0'); if (!--fl) { if (num) num = pbn_mulus_add(num, LIMB_DIGITS_VAL, lv); else num = pbn_uint(lv); lv = 0; fl = LIMB_DIGITS; } } } printf("e%+ld (%d digits)\n", exponent, digits); exponent2 = exponent5 = exponent; tz = pbn_ctz(num); exponent2 += tz; num = pbn_shr(num, tz); pbn_dump(stdout, num); printf(" * 2^%ld * 5^%ld\n", exponent2, exponent5); if (exponent5 < 0) { fpg->num = num; fpg->den = den = power_of_five(-exponent5); printf("num: "); pbn_dump(stdout, num); printf("\nden: "); pbn_dump(stdout, den); putchar('\n'); return FP_FRAC; } else { if (exponent5 > 0) num = pbn_mul(num, power_of_five(exponent5)); fpg->mnt = num; printf("mnt: "); pbn_dump(stdout, num); putchar('\n'); return FP_INT; } } /* * mbits is the number of true mantissa bits, excluding * the implicit 1. */ static int adjust_number(struct fpgen *fpg, int mbits) { struct pbn *num = fpg->num; struct pbn *den = fpg->den; int xbits, rem; if (den) { /* * Excess quotient bits: we need at least 2, plus we have a 1-bit * uncertainty in the quotient length */ mbits += 3; xbits = num->bits - den->bits + 1; if (xbits < mbits) { num = pbn_shl(num, mbits - xbits); fpg->pw2 -= mbits - xbits; } rem = pbn_div(&fpg->mnt, &fpg->num, num, pbn_ref(den)); printf("mnt: "); pbn_dump(stdout, fpg->mnt); printf(" [%d]\nexp: %d\nnum: ", rem, fpg->pw2); pbn_dump(stdout, fpg->num); printf("\nden: "); pbn_dump(stdout, den); putchar('\n'); } /* Left-shift the mantissa so the 1. falls in its own limb */ xbits = (1 - fpg->mnt->bits) & (PBN_LIMB_BITS-1); if (xbits) fpg->mnt = pbn_shl(fpg->mnt, xbits); printf("mnt: "); pbn_dump(stdout, fpg->mnt); putchar('\n'); /* * XXX: The output generator needs to know what bits in mnt are * actually significant, which depends on xbits among other things */ } int main(int argc, char *argv[]) { int i; struct fpgen fpg; enum fp_type ft; for (i = 1; i < argc; i++) { ft = analyze_string(argv[i], &fpg); if (ft >= FP_INT) adjust_number(&fpg, 52); } return 0; }