Experimatation floating-point

/**
 * Experimentations autour de nombre flottants en faisant varier le nombre de
 * bits de la mantisse.
 *
 * @author Samuel DEVULDER Avril 2016
 */
public class FP extends Number implements Comparable<FP> {
    // http://www.quinapalus.com/efunc.html

    private static final long serialVersionUID = 1L;
    byte sign, exp;
    long mantissa;

    static final int BITS = 15;

    static long ONE_MSK = 1l << BITS;

    static FP ZERO = new FP(0, 0, 0);
    static FP LN2 = new FP(1, -1, (long) (Math.log(2) * ONE_MSK * 2 + .5));
    static FP ONE = new FP(1, 0, 1 << BITS);
    static FP SQRT2 = new FP(1, 0, (long) (Math.sqrt(2) * ONE_MSK + .5));
    static FP TWO = new FP(1, 1, 1 << BITS);
    static FP LN10 = new FP(1, 1, (long) (Math.log(10) * ONE_MSK / 2 + .5));
    static FP E = new FP(1, 1, (long) (Math.exp(1) * ONE_MSK / 2 + .5));
    static FP PI = new FP(1, 1, (long) (Math.PI * ONE_MSK / 2 + .5));
    static FP TEN = new FP(1, 3, 10 << (BITS - 3));
    static FP MAX = new FP(1, 127, ONE_MSK * 2 - 1);

    public static void main(String[] args) {
        System.out.println("LN2=" + LN2);
        System.out.println("ONE=" + ONE);
        System.out.println("SQRT2=" + SQRT2);
        System.out.println("TWO=" + TWO);
        System.out.println("LN10=" + LN10);
        System.out.println("PI=" + PI);
        System.out.println("TEN=" + TEN);
        System.out.println("MAX=" + MAX);
        System.out.println("ONE+SQRT2=" + ONE.add(SQRT2));
        System.out.println("ONE-SQRT2=" + ONE.sub(SQRT2));
        System.out.println("SQRT2-SQRT2=" + SQRT2.sub(SQRT2));
        System.out.println("SQRT2*SQRT2=" + SQRT2.mul(SQRT2));
        System.out.println("SQRT2/SQRT2=" + SQRT2.div(SQRT2));
        System.out.println("ONE/MAX=" + ONE.div(MAX));
        System.out.println("parse1=" + new FP("-314.1592653589E-2"));
        System.out.println("parse2=" + new FP("3.141592653589"));

        System.out.println("ln(2)=" + TWO.ln());
        System.out.println("ln(10)=" + TEN.ln());

        System.out.println("log10(2)=" + TWO.log10());
        System.out.println("log10(10)=" + TEN.log10());

        System.out.println("sqrt(2)=" + TWO.sqrt());
        System.out.println("sqrt(10)=" + TEN.sqrt());

        System.out.println("exp(-2)=" + TWO.neg().exp());
        System.out.println("exp(-1)=" + ONE.neg().exp());
        System.out.println("exp(0)=" + ZERO.exp());
        System.out.println("exp(1)=" + ONE.exp());
        System.out.println("exp(2)=" + TWO.exp());
        System.out.println("exp(ln2)=" + LN2.exp());
        System.out.println("exp(ln10)=" + LN10.exp());
    }

    FP(int sign, int exp, long mantissa) {
        build(sign, exp, mantissa);
    }

    protected int sign() {
        return sign;
    }

    protected int exponent() {
        return exp;
    }

    protected long mantissa() {
        return mantissa;
    }

    public FP(double d) {
        int sgn = d < 0 ? -1 : d > 0 ? 1 : 0;
        int exp = Math.getExponent(d);
        int man = (int) Math.scalb(sgn < 0 ? -d : d, BITS - exp);
        build(man == 0 ? 0 : sgn, exp, man);
    }

    protected FP build(int sign, int exp, long mantissa) {
        this.sign = (byte) sign;
        this.exp = (byte) exp;
        this.mantissa = mantissa & (ONE_MSK * 2 - 1);
        if (sign != 0 && (mantissa & ONE_MSK) == 0)
            throw new RuntimeException("Bad mantissa "
                    + Long.toHexString(this.mantissa));
        return this;
    }

    public FP(String s_) {
        int i = 0;
        char[] s = s_.toCharArray();

        // read sign
        int sgn = 1;
        while (i < s.length && (s[i] == '-' || s[i] == '+')) {
            if (s[i++] == '-')
                sgn = -sgn;
        }

        // read int part
        FP v = ZERO;
        while (i < s.length && s[i] >= '0' && s[i] <= '9') {
            v = v.mul(TEN).add(s[i++] - '0');
        }

        // read decimal
        if (i < s.length && s[i] == '.') {
            FP div = ONE;
            ++i;
            while (i < s.length && s[i] >= '0' && s[i] <= '9') {
                div = div.div(TEN);
                v = v.add(div.mul(s[i++] - '0'));
            }
        }

        // read exponent
        if (i < s.length && (s[i] == 'e' || s[i] == 'E')) {
            ++i;

            int sgn2 = 1;
            while (i < s.length && (s[i] == '-' || s[i] == '+')) {
                if (s[i++] == '-')
                    sgn2 = -sgn2;
            }

            FP e = ONE;
            while (i < s.length && s[i] >= '0' && s[i] <= '9') {
                e = e.pow(10).mul(pow_tab[s[i++] - '0']);
            }

            v = sgn2 < 0 ? v.div(e) : v.mul(e);
        }

        if (sgn < 0)
            v = v.neg();

        build(v.sign(), v.exponent(), v.mantissa());
    }

    static FP[] pow_tab = { TEN.pow(0), TEN.pow(1), TEN.pow(2), TEN.pow(3),
            TEN.pow(4), TEN.pow(5), TEN.pow(6), TEN.pow(7), TEN.pow(8),
            TEN.pow(9) };

    @Override
    public boolean equals(Object obj) {
        return obj instanceof FP && compareTo((FP) obj) == 0;
    }

    @Override
    public int hashCode() {
        if (sign == 0)
            return 0;
        return (int) ((sign * 31 + mantissa) * 31 + exp);
    }

    @Override
    public double doubleValue() {
        if (sign() == 0)
            return 0;
        double t = mantissa() * Math.pow(2, exponent() - BITS);
        return sign() < 0 ? -t : t;
    }

    @Override
    public float floatValue() {
        return (float) doubleValue();
    }

    @Override
    public int intValue() {
        return (int) longValue();
    }

    @Override
    public long longValue() {
        int sgn = sign();
        if (sgn == 0)
            return 0;
        long val = mantissa() >>> (BITS - exponent());
        return sgn < 0 ? -val : val;
    }

    public FP floor() {
        int exp = exponent();
        if (exp < 0)
            return ZERO;
        if (exp > BITS)
            return this;
        long m = (2l << BITS) - ((1l << BITS) >>> exp);
        return new FP(sign(), exp, mantissa() & m);
    }

    static FP PRECISION;

    static {
        FP x = ONE, y = new FP(1, -BITS, 1 << BITS);
        while (x.compareTo(y) > 0) {
            PRECISION = x;
            x = x.div(TEN);
        }
    }

    public String toString() {
        // return Double.toString(doubleValue());

        // http://stackoverflow.com/questions/2302969/how-to-implement-char-ftoafloat-num-without-sprintf-library-function-i
        if (sign() == 0)
            return "0";
        StringBuilder sb = new StringBuilder();

        boolean neg = sign() < 0;
        FP x = this.abs();
        // if(x.mantissa!=0xFFFF) ++x.mantissa;
        int m = 0, m1 = 0;
        // while(TEN.pow(m+1).compareTo(x)<=0) ++m;
        m = x.log10().intValue();
        boolean useExp = (m >= 5 || (neg && m >= 3) || m <= -5);
        if (neg)
            sb.append('-');
        // set up for scientific notation
        if (useExp) {
            if (m < 0)
                --m;
            x = x.mul(TEN.pow(-m));
            m1 = m;
            m = 0;
        }
        if (m < 1.0) {
            m = 0;
        }
        // convert the number
        while (x.compareTo(PRECISION) > 0 || m >= 0) {
            FP weight = TEN.pow(m);
            if (weight.tst() > 0) {
                int digit = x.div(weight).intValue();
                x = x.sub(weight.mul(digit));
                sb.append((char) ('0' + digit));
            }
            if (m == 0 && x.tst() > 0)
                sb.append('.');
            m--;
        }
        if (useExp) {
            // convert the exponent
            sb.append('e');
            if (m1 > 0) {
                sb.append('+');
            } else {
                sb.append('-');
                m1 = -m1;
            }
            if (m1 > 10) {
                sb.append((char) ('0' + m1 / 10));
                m1 = m1 % 10;
            }
            sb.append((char) ('0' + m1));
        }
        return sb.toString();
    };

    public FP neg() {
        return new FP(-sign(), exponent(), mantissa());
    }

    public FP abs() {
        return new FP(sign() == 0 ? 0 : 1, exponent(), mantissa());
    }

    public FP mulSign(FP other) {
        return new FP((other.sign() & sign()) == 0 ? 0
                : ((other.sign() ^ sign()) | 1), exponent(), mantissa());
    }

    public int tst() {
        return sign();
    }

    public int compareTo(FP other) {
        int d = this.sign() - other.sign();
        if (d != 0) {
            return d;
        }
        if (this.sign() == 0)
            return 0;
        d = this.exponent() - other.exponent();
        if (d != 0)
            return d;
        d = (int) (this.mantissa() - other.mantissa());
        return d;
    }

    public FP mul(FP other) {
        if ((this.sign() & other.sign()) == 0)
            return ZERO;
        int sign = (this.sign() ^ other.sign()) | 1;
        int exp = this.exponent() + other.exponent();
        long mant = (this.mantissa() * other.mantissa()) >>> BITS;
        if ((mant & (2l << BITS)) != 0) {
            mant >>>= 1;
            ++exp;
        }
        return new FP(mant == 0 ? 0 : sign, exp, mant);
    }

    public FP mul(double v) {
        return mul(new FP(v));
    }

    public FP square() {
        return mul(this);
    }

    public FP pow(int n) {
        boolean neg = n < 0;
        if (neg)
            n = -n;
        FP x = this, y = ONE;
        while (n != 0) {
            if ((n & 1) != 0) {
                y = y.mul(x);
            }
            x = x.square();
            n >>>= 1;
        }
        return neg ? y.inv() : y;
    }

    public FP div(FP other) {
        if (other.sign() == 0)
            return MAX.mulSign(this).mulSign(other);
        if (this.sign() == 0)
            return ZERO;
        int sign = (this.sign() ^ other.sign()) | 1;
        int exp = this.exponent() - other.exponent();
        long mant = (this.mantissa() << (1 + BITS)) / other.mantissa();
        if ((mant & (2l << BITS)) != 0) {
            mant >>>= 1;
        } else {
            --exp;
        }
        return new FP(sign, exp, mant);
    }

    public FP inv() {
        return ONE.div(this);
    }

    public FP sub(FP other) {
        return add(other.neg());
    }

    public FP add(FP other) {
        if (sign() == 0)
            return other;
        if (other.sign() == 0)
            return this;

        FP small = this;
        FP big = other;
        if (small.exponent() == big.exponent()) {
            if (small.mantissa() > big.mantissa()) {
                small = other;
                big = this;
            }
        } else if (small.exponent() > big.exponent()) {
            small = other;
            big = this;
        }
        int exp = big.exponent();
        long mant = big.mantissa();
        long delta = small.mantissa() >>> (big.exponent() - small.exponent());
        if (big.sign() == small.sign()) {
            mant += delta;
            if ((mant & (2l << BITS)) != 0) {
                mant >>>= 1;
                ++exp;
            }
        } else {
            mant -= delta;
            if (mant == 0)
                return ZERO;
            while ((mant & (1l << BITS)) == 0) {
                mant <<= 1;
                --exp;
            }
        }
        return new FP(big.sign(), exp, mant);
    }

    protected FP add(double v) {
        return add(new FP(v));
    }

    static protected FP[] ln = { new FP(8.690628917915396e-06),
            new FP(0.9992996219306662), new FP(-0.4907434587925273),
            new FP(0.2867076984634543), new FP(-0.1332213692732176),
            new FP(0.03110469698575859) };

    public FP ln() {
        if (sign() <= 0)
            return MAX.neg();
        FP x = new FP(1, 0, mantissa()).sub(ONE), y = ZERO;
        for (int i = ln.length; --i >= 0;)
            y = y.mul(x).add(ln[i]);
        int exp = this.exponent() < 0 ? -this.exponent() : this.exponent();
        while (exp > 0) {
            y = y.add(LN2);
            --exp;
        }
        if (this.exponent() < 0)
            y = y.neg();
        return y;
    }

    public FP log10() {
        return ln().div(LN10);
    }

    public FP sqrt() {
        if (sign() <= 0)
            return ZERO;

        long op = mantissa() << BITS;
        long one = 1l << (BITS + BITS);
        long res = 0;
        while (one != 0) {
            if (op >= res + one) {
                op -= res + one;
                res += one << 1;
            }
            res >>>= 1;
            one >>>= 2;
        }
        /* Do arithmetic rounding to nearest integer */
        if (op > res) {
            res++;
        }

        FP x = new FP(1, this.exponent() >> 1, res);
        if ((this.exponent() & 1) != 0)
            x = x.mul(SQRT2);
        return x;
    }

    public FP exp() {
        // exp(2^x * y) = exp(y)^(2^x) = exp(2^0) * exp(2^-1)

        if (sign() == 0)
            return ONE;

        // exp(y)
        FP x = ONE;
        if (false) {
            FP y = E;
            long m = mantissa(), msk = ONE_MSK;
            while (m != 0 && !ONE.equals(y)) {
                if ((m & msk) != 0) {
                    x = x.mul(y);
                    m ^= msk;
                }
                y = y.sqrt();
                msk >>>= 1;
            }
        } else {
            // optim
            long z = mantissa(), y = ONE_MSK, t, k;
            int exp = 0;
            for (k = 1; (t = (long) (Math.log((k + 1.0) / k) * ONE_MSK + 0.5)) != 0; k <<= 1) {
                t = z - t;
                if (t >= 0) {
                    z = t;
                    y += y / k;
                    if ((y & (ONE_MSK * 2)) != 0) {
                        y >>>= 1;
                        ++exp;
                    }
                }
            }
            x = new FP(1, exp, y);
        }

        // exp(2^exponent)
        int exp = exponent();
        while (exp > 0 && !MAX.equals(x)) {
            x = x.square();
            --exp;
        }
        while (exp < 0 && !ONE.equals(x)) {
            x = x.sqrt();
            ++exp;
        }

        return sign() < 0 ? x.inv() : x;
    }
}