import java.util.ArrayList;
import java.util.List;

/**
 * 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;
    int exp;
    long mantissa;

    static final int EXPONENT_BITS = 8;
    static final int MANTISSA_BITS = 15;

    static long ONE_MSK = 1l << MANTISSA_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 << MANTISSA_BITS);
    static FP SQRT2 = new FP(1, 0, (long) (Math.sqrt(2) * ONE_MSK + .5));
    static FP TWO = new FP(1, 1, 1 << MANTISSA_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 << (MANTISSA_BITS - 3));
    static FP MAX = new FP(1, 127, ONE_MSK * 2 - 1);

    static FP[] NUMS = { ZERO, LN2, ONE, SQRT2, TWO, LN10, E, PI, TEN, MAX };

    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("SQRT2 % 1=" + SQRT2.mod(ONE));
        System.out.println("SQRT2*TEN % LN2=" + SQRT2.mul(TEN).mod(LN2));

        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"));

        for (FP x : NUMS) {
            System.out.println("ln(" + x + ")=" + x.ln());
            System.out.println("log10(" + x + ")=" + x.log10());
            System.out.println("sqrt(" + x + ")=" + x.sqrt());
            System.out.println("" + x + "²=" + x.square());
            System.out.println("exp(" + x.neg() + ")=" + x.neg().exp());
            System.out.println("exp(" + x + ")=" + x.exp());
            System.out.println("sin(" + x + ")=" + x.sin());
            System.out.println("cos(" + x + ")=" + x.cos());
            System.out.println("sin²+cos²(" + x + ")=" + x.cos().square().add(x.sin().square()));
            System.out.println("atan(" + x + ")=" + x.atan());
        }

        // for (FP x = ONE; true;) {
        // System.out.println(x);
        // FP y = x.add(x);
        // if (x.equals(y))
        // break;
        // x = y;
        // }

        for (long mant = ONE_MSK; mant < 2 * ONE_MSK; ++mant) {
            FP x = new FP(1, 0, mant);
            x.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);
        long man = (long) Math.scalb(sgn < 0 ? -d : d, MANTISSA_BITS - exp);
        build(man == 0 ? 0 : sgn, exp, man);
    }

    protected FP build(int sign, int exp, long mantissa) {
        this.sign = (byte) sign;

        int EXP_MSK = (1 << (EXPONENT_BITS - 1)) - 1;
        this.exp = (byte) (exp > EXP_MSK ? EXP_MSK : exp < -EXP_MSK ? -EXP_MSK : 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() - MANTISSA_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() >>> (MANTISSA_BITS - exponent());
        return sgn < 0 ? -val : val;
    }

    public FP mul2() {
        return new FP(sign(), exponent() + 1, mantissa());
    }

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

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

    static FP PRECISION;

    static {
        FP x = ONE, y = new FP(1, -MANTISSA_BITS, 1 << MANTISSA_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()) >>> MANTISSA_BITS;
        if ((mant & (2l << MANTISSA_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 + MANTISSA_BITS)) / other.mantissa();
        if ((mant & (2l << MANTISSA_BITS)) != 0) {
            mant >>>= 1;
        } else {
            --exp;
        }
        return new FP(sign, exp, mant);
    }

    public FP mod(FP other) {
        return sub(div(other).floor().mul(other));
    }

    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 << MANTISSA_BITS)) != 0) {
                mant >>>= 1;
                ++exp;
            }
        } else {
            mant -= delta;
            if (mant == 0)
                return ZERO;
            while ((mant & (1l << MANTISSA_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 y = ZERO;
        if (false) {
            // table accurate up to 15 bits
            FP x = new FP(1, 0, mantissa()).sub(ONE);
            for (int i = ln.length; --i >= 0;)
                y = y.mul(x).add(ln[i]);
        } else {
            // exact algorithm log( [1..2[ )
            // multiply mantissa without geting out of the [1..2[
            // interval
            long x = mantissa(), z = exp_tab[0];
            for (int k = 1; k < exp_tab.length; ++k) {
                long t = x + (x >> k);
                if ((t & 2 * ONE_MSK) == 0) {
                    x = t;
                    z -= exp_tab[k];
                }
            }
            int exp = 0;
            if (z != 0)
                while (z < ONE_MSK) {
                    z <<= 1;
                    --exp;
                }
            y = new FP(z == 0 ? 0 : 1, exp, z);
        }
        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() << MANTISSA_BITS;
        long one = 1l << (MANTISSA_BITS + MANTISSA_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;
    }

    static long[] exp_tab;

    static {
        List<Long> array = new ArrayList<Long>();
        for (long k = 1, t = 0; (t = (long) (Math.log((k + 1.0) / k) * ONE_MSK + 0.5)) != 0; k <<= 1)
            array.add(t);
        exp_tab = new long[array.size()];
        for (int i = exp_tab.length; --i >= 0;)
            exp_tab[i] = array.get(i).longValue();
    }

    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;
            int exp = 0;
            for (int k = 0; k < exp_tab.length; ++k) {
                long t = z - exp_tab[k];
                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;
    }

    static protected FP[] sin = { new FP(7.068518675814434e-06), new FP(0.9996898644339388), new FP(0.002193716170952744), new FP(-0.1722388650880223),
            new FP(0.006097383673278212), new FP(0.005721724054854694) };

    public FP sin() {
        boolean neg = sign() < 0;
        FP x = abs().mod(PI.mul2());
        if (x.compareTo(PI) >= 0) {
            x = x.sub(PI.mul2());
            if (x.sign() < 0) {
                neg = !neg;
                x = x.neg();
            }
        }
        if (x.compareTo(PI.div2()) >= 0) {
            x = PI.sub(x);
        }
        FP y = ZERO;
        // table accurate up to 15 bits
        for (int i = sin.length; --i >= 0;)
            y = y.mul(x).add(sin[i]);
        if (neg)
            y = y.neg();
        return y;
    }

    static protected FP[] cos = { new FP(1.000007068162554), new FP(-0.0003343311503654495), new FP(-0.4974330510416236), new FP(-0.007249836999713618),
            new FP(0.05103555758070047), new FP(-0.005721683674851457) };

    public FP cos() {
        // return add(PI.div2()).sin();
        boolean neg = false;
        FP x = abs().mod(PI.mul2());
        if (x.compareTo(PI) >= 0) {
            x = PI.mul2().sub(x);
        }
        if (x.compareTo(PI.div2()) >= 0) {
            x = PI.sub(x);
            neg = true;
        }
        FP y = ZERO;
        // table accurate up to 15 bits
        for (int i = cos.length; --i >= 0;)
            y = y.mul(x).add(cos[i]);
        if (neg)
            y = y.neg();
        return y;
    }

    public FP tan() {
        return sin().div(cos());
    }

    public FP atan() {
        return atan2(ONE);
    }

    static protected FP[] atan = { new FP(-0.0464964749), new FP(0.15931422), new FP(-0.327622764) };

    public FP atan2(FP x) {
        // http://math.stackexchange.com/questions/1098487/atan2-faster-approximation
        FP y = this;
        FP a = x.abs(), b = y.abs();
        if (a.compareTo(b) > 0) {
            a = b;
            b = x.abs();
        }
        FP s = a.square();
        FP r = atan[0].mul(s).add(atan[1]).mul(s).add(atan[2]).mul(s).mul(a).add(a);
        if (y.abs().compareTo(x.abs()) > 0)
            r = PI.div2().sub(r);
        if (x.sign() < 0)
            r = PI.sub(r);
        if (y.sign() < 0)
            r = r.neg();
        return r;
    }

}