Fib.java

public class Fib {
    public static void main(String[] args) {
        long[] js = new long[] {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,50,100,300,1000,10000,1000000,10000000L};
        for(long j : js) {
            long st = System.nanoTime();
            String result = fib(j);
            long tm = System.nanoTime() - st;
            System.out.println("fib("+j+") == "+result+";time="+tm/1000000000+"s+"+(tm/1000000)%1000+"ms");
            System.gc();
        }
    }
    public static String fib(long num) {
        //fib(num): num
        //fib(num) ~= (1+rt(5)/2)^num
        //log 10000 fib(num) ~= log(1+rt(5)/2)*num/log(10000)
        //lg 10000 fib(num) ~= log(1+rt(5)/2)*num/log(10000)*log(2)
        FFTBI bi = new FFTBI((int)(num/20+1000));
        BigIntM a = bi.newInst(0);
        BigIntM b = bi.newInst(1);
        BigIntM c = bi.newInst(0);
        BigIntM d = bi.newInst(0);
        int lgnum = 2; while((1L<<lgnum)<num)lgnum++;
        for(int i = lgnum; i >= 0; i--) {
            if(((num>>i)&1)==0) {
                //(a,b) = ((2*fn1-fn) * fn, fn1*fn1 + fn*fn)
                //(a,b) <- ((2*b-a) * a, b*b + a*a)
                c.copyM(a);
                c.mulM(a);
                d.copyM(b);
                d.mulM(b);
                d.addM(c);

                c.copyM(b);
                c.addM(b);
                c.diffM(a);
                c.mulM(a);

                BigIntM t;
                t=a;a=c;c=t;
                t=b;b=d;d=t;
            } else {
                //(a,b) = (fn*fn+fn1*fn1, 2*fn1*fn +fn1*fn1)
                //(a,b) <- (a*a+b*b, (2*a+b)*b)
                c.copyM(a);
                c.mulM(a);
                d.copyM(b);
                d.mulM(b);
                c.addM(d);

                d.copyM(a);
                d.addM(a);
                d.addM(b);
                d.mulM(b);

                BigIntM t;
                t=a;a=c;c=t;
                t=b;b=d;d=t;
            }
        }
        return a.toString();
    }
}

abstract class BigIntM {
    abstract void copyM(BigIntM t);
    abstract void addM(BigIntM t);
    abstract void diffM(BigIntM t);
    abstract void mulM(BigIntM t);
}

class FFTBI {
    private int n;
    private int lgn;
    private int[] brev;
    private double[] tbl_r;
    private double[] tbl_i;
    private double[] swp_r;
    private double[] swp_i;

    private double[] mul0_r;
    private double[] mul0_i;
    private double[] mul1_r;
    private double[] mul1_i;
    public BigIntM newInst(int num) {
        return new BigIntMI(num);
    }
    private class BigIntMI extends BigIntM {
        private FFTBI parent;
        //10000-based(Decimal)
        private int[] data;
        BigIntMI() {
            parent = FFTBI.this;
            data = new int[n];
        }
        //0<=num<100000000
        BigIntMI(int num) {
            parent = FFTBI.this;
            data = new int[n];
            data[0] = (num%10000);
            data[1] = (num/10000)%10000;
        }
        void copyM(BigIntM t) {
            BigIntMI ti = (BigIntMI)t;
            if(ti.parent!=parent)throw new RuntimeException("parent not matched");
            System.arraycopy(ti.data,0,data,0,n);
        }
        void addM(BigIntM t) {
            BigIntMI ti = (BigIntMI)t;
            if(ti.parent!=parent)throw new RuntimeException("parent not matched");
            int cl=0;
            for(int i = 0; i < n; i++) {
                data[i] += ti.data[i]+cl;
                cl = data[i]/10000;
                data[i] = data[i]%10000;
            }
        }
        void diffM(BigIntM t) {
            BigIntMI ti = (BigIntMI)t;
            if(ti.parent!=parent)throw new RuntimeException("parent not matched");
            int cl=0;
            for(int i = 0; i < n; i++) {
                data[i] -= ti.data[i]+cl;
                cl = 0;
                if(data[i]<0) {
                    cl = 1;
                    data[i]+=10000;
                }
            }
        }
        void mulM(BigIntM t) {
            BigIntMI ti = (BigIntMI)t;
            if(ti.parent!=parent)throw new RuntimeException("parent not matched");
            synchronized(parent) {
                for(int i = 0; i < n; i++) {
                    mul0_r[i] = data[i];
                    mul0_i[i] = 0;
                    mul1_r[i] = ti.data[i];
                    mul1_i[i] = 0;
                }
                applyFFT(mul0_r, mul0_i, false);
                applyFFT(mul1_r, mul1_i, false);
                for(int i = 0; i < n; ++i) {
                    double c_r = mul0_r[i]*mul1_r[i] - mul0_i[i]*mul1_i[i];
                    double c_i = mul0_r[i]*mul1_i[i] + mul0_i[i]*mul1_r[i];
                    mul0_r[i] = c_r;
                    mul0_i[i] = c_i;
                }
                applyFFT(mul0_r, mul0_i, true);
                long tmp = 0;
                for(int i = 0; i < n; ++i) {
                    tmp += (long)(mul0_r[i]+0.5);
                    data[i] = (int)(tmp%10000);
                    tmp /= 10000;
                }
            }
        }
        public String toString() {
            boolean flag=false;
            StringBuffer buf = new StringBuffer();
            for(int i = n-1; i >= 0; --i) {
                if(flag) {
                    buf.append((char)('0'+(data[i]/1000%10)));
                    buf.append((char)('0'+(data[i]/100%10)));
                    buf.append((char)('0'+(data[i]/10%10)));
                    buf.append((char)('0'+(data[i]/1%10)));
                } else if(data[i]!=0) {
                    flag=true;
                    if(data[i]>=1) {
                        if(data[i]>=10) {
                            if(data[i]>=100) {
                                if(data[i]>=1000) {
                                    buf.append((char)('0'+(data[i]/1000%10)));
                                }
                                buf.append((char)('0'+(data[i]/100%10)));
                            }
                            buf.append((char)('0'+(data[i]/10%10)));
                        }
                        buf.append((char)('0'+(data[i]/1%10)));
                    }
                }
            }
            if(buf.length()==0)buf.append("0");
            return buf.toString();
        }
    }
    private void applyFFT(double[] dat_r, double[] dat_i, boolean inverse) {
        double per = 1.0;
        if(inverse)per /= n;

        for(int j =0; j < n; ++j) {
            swp_r[j] = dat_r[j];
            swp_i[j] = dat_i[j];
        }

        for(int i = lgn-1; i >= 0; --i) {
            int l = 1<<i;
            for(int j = 0; j < n; ++j) {
                if((j&l)==0) {
                    double t0_r = swp_r[j] + tbl_r[(j>>i)] * swp_r[j^l] -
                        tbl_i[(j>>i)] * swp_i[j^l];
                    double t0_i = swp_i[j] + tbl_r[(j>>i)] * swp_i[j^l] +
                        tbl_i[(j>>i)] * swp_r[j^l];
                    double t1_r = swp_r[j] + tbl_r[(j>>i)+1] * swp_r[j^l] -
                        tbl_i[(j>>i)+1] * swp_i[j^l];
                    double t1_i = swp_i[j] + tbl_r[(j>>i)+1] * swp_i[j^l] +
                        tbl_i[(j>>i)+1] * swp_r[j^l];
                    swp_r[j] = t0_r;
                    swp_i[j] = t0_i;
                    swp_r[j^l] = t1_r;
                    swp_i[j^l] = t1_i;
                }
            }
        }
        for(int j = 0; j < n; ++j) {
            dat_r[j] = swp_r[brev[j]]*per;
            dat_i[j] = -swp_i[brev[j]]*per;
        }
    }
    FFTBI(int n) {
        while((1<<lgn)<n)lgn++;
        this.n=n=(1<<lgn);
        brev = new int[n];
        tbl_r = new double[n];
        tbl_i = new double[n];
        swp_r = new double[n];
        swp_i = new double[n];
        mul0_r = new double[n];
        mul0_i = new double[n];
        mul1_r = new double[n];
        mul1_i = new double[n];
        //create bit-reverse table.
        for(int j = 0; j < lgn; j++) {
            for(int i = 0; i < n; i++) {
                brev[i|(1<<j)] |= (1<<(lgn-1-j));
            }
        }
        //create exp table.
        for(int i = 0; i < n; ++i) {
            tbl_r[brev[i]] = Math.cos(2*Math.PI*i/n);
            tbl_i[brev[i]] = Math.sin(2*Math.PI*i/n);
        }
    }
}