summaryrefslogtreecommitdiffstats
path: root/fpconv.c
diff options
context:
space:
mode:
Diffstat (limited to 'fpconv.c')
-rw-r--r--fpconv.c289
1 files changed, 289 insertions, 0 deletions
diff --git a/fpconv.c b/fpconv.c
new file mode 100644
index 0000000..92e6f85
--- /dev/null
+++ b/fpconv.c
@@ -0,0 +1,289 @@
+/*
+ * 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 <stdio.h>
+#include <inttypes.h>
+#include <stdlib.h>
+#include <time.h>
+#include <sys/time.h>
+#include <ctype.h>
+#include <string.h>
+#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;
+}