/*
 * 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;
}