public class MillerRabin {

    public static void main(String[] args) {
        long t0, t1;
        t0 = System.nanoTime();
        boolean b = !isComposite(479001599);
        boolean c = !isComposite(456789012);
        t1 = System.nanoTime();
        System.out.println(t1 - t0);
        System.out.println(b + " " + c);
    }

    /**
     * <p>Miller-Rabin 测试某一个数是否是合数</p>
     *
     * @param n 需要测试的数
     * @return true: 该数为合数；false: 该数为素数
     */
    public static boolean isComposite(int n) {
        if (n < 2) {
            throw new IllegalArgumentException("number must greater than or equals 2");
        }
        // 排除 2、3、5、7 以加速测试
        if (n == 2 || n == 3 || n == 5 || n == 7) {
            return false;
        }

        // 偶数
        if ((n & 1) == 0) {
            return true;
        }

        // 排除 3、5、7 的倍数，以加速测试
        if (n % 3 == 0) {
            return true;
        }
        if (n % 5 == 0) {
            return true;
        }
        if (n % 7 == 0) {
            return true;
        }

        // 寻找 s 和 d 以满足 n = 2^s * d + 1
        int s = 0, d = n - 1;
        while ((d & 1) == 0) {
            d >>= 1;
            s++;
        }

        // 对于各种数值需要进行 Miller-Rabin 基准测试的素数值
        // 参考：http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Deterministic_variants_of_the_test
        // if n < 1,373,653, it is enough to test a = 2 and 3;
        // if n < 9,080,191, it is enough to test a = 31 and 73;
        // if n < 4,759,123,141, it is enough to test a = 2, 7, and 61;
        // if n < 2,152,302,898,747, it is enough to test a = 2, 3, 5, 7, and 11;
        // if n < 3,474,749,660,383, it is enough to test a = 2, 3, 5, 7, 11, and 13;
        // if n < 341,550,071,728,321, it is enough to test a = 2, 3, 5, 7, 11, 13, and 17.

        if (n < 1373653) {
            if (loopMillerRabin(s, d, n, 2, 3)) {
                return true;
            }
        } else if (n < 9080191) {
            if (loopMillerRabin(s, d, n, 31, 73)) {
                return true;
            }
        } else {
            // 4,759,123,141 已经超过 int 的最大值，因此大于等于 9080191 就采用 4,759,123,141 的基准测试
            if (loopMillerRabin(s, d, n, 2, 7, 61)) {
                return true;
            }
        }
        return false;
    }

    /**
     * <p>循环 Miller-Rabin 测试</p>
     *
     * @param s n = 2^s * d + 1 中的 s 值
     * @param d n = 2^s * d + 1 中的 d 值
     * @param n 需要测试的数
     * @param t 测试的基准素数
     */
    private static boolean loopMillerRabin(int s, int d, int n, int... t) {
        for (int i = 0; i < t.length; i++) {
            if (testMillerRabin(t[i], s, d, n)) {
                return true;
            }
        }
        return false;
    }

    /**
     * <p>Miller-Rabin 基本测试</p>
     *
     * @param a 素性测试基准素数
     * @param s n = 2^s * d + 1 中的 s 值
     * @param d n = 2^s * d + 1 中的 d 值
     * @param n 需要测试的数
     * @return 测试某一数是否是合数。true: 该数是合数；false: 该数可能是素数。若返回 false
     * 需要进行多基准的联合测试才能判断该数确实是素数
     */
    private static boolean testMillerRabin(int a, int s, int d, int n) {
        if (montgomery(a, d, n) != 1) {
            int e = 1;
            for (int i = 0; i < s; i++) {
                if (montgomery(a, d * e, n) + 1 == n) {
                    return false;
                }
                e <<= 1;
            }
            return true;
        }
        return false;
    }

    /**
     * <p>使用 Montgomery 算法计算 (base ^ exp) % mod 的值。</p>
     * 
     * <p>由于 Java 中 int 的运算速度远远大于 long 的运算速度，因此该算法需要改进！</p>
     *
     * @param base 基数
     * @param exp 指数
     * @param mod 模数
     */
    private static int montgomery(int base, int exp, int mod) {
        if (base > 46340 || mod > 46340) {
            long temp = 1;
            long prod = base % mod;
            while (exp > 1) {
                if ((exp & 1) != 0) {
                    temp = (temp * prod) % mod;
                }
                prod = (prod * prod) % mod;
                exp >>= 1;
            }
            return (int) ((temp * prod) % mod);
        } else {
            int temp = 1;
            int prod = base % mod;
            while (exp > 1) {
                if ((exp & 1) != 0) {
                    temp = (temp * prod) % mod;
                }
                prod = (prod * prod) % mod;
                exp >>= 1;
            }
            return (temp * prod) % mod;
        }
    }

    /**
     * <p>
     * 根据复化辛普生法则计算高斯 Li 函数。Li(x) 近似于素数个数函数 π(x)
     * </p>
     * 
     * @param x 数值范围
     * @return 该值以内的素数个数的估计值
     */
    private static double gaussLi(int x) {
        int n = x;
        double s = 2;
        double h = (x - s) / n;
        double a = f(s + h / 2);
        double b = 0;
        for (int k = 1; k < n; k++) {
            a += f(s + k * h + h / 2);
            b += f(s + k * h);
        }
        return (h / 6) * (f(s) + 4 * a + 2 * b + f(x));
    }

    /**
     * <p>
     * Guass Li(x) 积分函数项
     * </p>
     * 
     * @param x
     * @return
     */
    private static double f(double x) {
        return 1 / Math.log(x);
    }
}