Skip to main content
Post Made Community Wiki by C. K. Young
Source Link
14mRh4X0r
  • 348
  • 1
  • 5

Java - 125.15 (21,400,000 / 171,000)

Also shamelessly copied from Peter Luschny's Github repo (thanks @semi-extrinsic) and licensed under the MIT license, this uses the "prime factorization nested squaring" algorithm as proposed by Albert Schönhage et al. (according to Luschny's factorial algorithms description page).

I slightly adapted the algorithm to use Java's BigInteger and to not use a lookup table for n < 20.

Compiled with gcj, which uses GMP for its BigInteger implementation, and ran on Linux 3.12.4 (Gentoo), on a Core i7 4700MQ at 2.40GHz

import java.math.BigInteger;

public class PrimeSieveFactorialSchoenhage {

    private static int[] primeList, multiList;

    public static BigInteger factorial(int n) {
        int log2n = 31 - Integer.numberOfLeadingZeros(n);
        int piN = log2n < 2 ? 1 : 2 + (15 * n) / (8 * (log2n - 1));

        primeList = new int[piN];
        multiList = new int[piN];

        int len = primeFactors(n);
        return nestedSquare(len).shiftLeft(n - Integer.bitCount(n));
    }

    private static BigInteger nestedSquare(int len) {
        if (len == 0) {
            return BigInteger.ONE;
        }

        int i = 0, mult = multiList[0];

        while (mult > 1) {
            if ((mult & 1) == 1) { // is mult odd ?
                primeList[len++] = primeList[i];
            }

            multiList[i++] = mult / 2;
            mult = multiList[i];
        }
        BigInteger ns = nestedSquare(i);
        if (len <= i) {
            return ns.multiply(ns);
        }

        return product(primeList, i, len - i).multiply(ns.multiply(ns));
    }

    private static BigInteger product(int[] a, int start, int length) {
        if (length == 0) {
            return BigInteger.ONE;
        }

        int len = (length + 1) / 2;
        long[] b = new long[len];

        int i, j, k;

        for (k = 0, i = start, j = start + length - 1; i < j; i++, k++, j--) {
            b[k] = a[i] * (long) a[j];
        }

        if (i == j) {
            b[k++] = a[j];
        }

        return recProduct(b, 0, k - 1);
    }

    private static BigInteger recProduct(long[] s, int n, int m) {
        if (n > m) {
            return BigInteger.ONE;
        }
        if (n == m) {
            return BigInteger.valueOf(s[n]);
        }
        int k = (n + m) >> 1;
        return recProduct(s, n, k).multiply(recProduct(s, k + 1, m));
    }

    private static int primeFactors(int n) {
        int[] primes = new int[n < 17 ? 6 : (int) Math.floor(n / (Math.log(n) - 1.5))];
        int numPrimes = makePrimeList(n, primes);

        int maxBound = n / 2, count = 0;

        int start = indexOf(primes, 2, 0, numPrimes - 1);
        int end = indexOf(primes, n, start, numPrimes);

        for (int i = start; i < end; i++) {
            int prime = primes[i];
            int m = prime > maxBound ? 1 : 0;

            if (prime <= maxBound) {
                int q = n;
                while (q >= prime) {
                    m += q /= prime;
                }
            }

            primeList[count] = prime;
            multiList[count++] = m;
        }
        return count;
    }

    private static int indexOf(final int[] data, int value, int low, int high) {
        while (low < high) {
            int mid = (low + high) >>> 1;

            if (data[mid] < value) {
                low = mid + 1;
            } else {
                high = mid;
            }
        }

        if (low >= data.length) {
            return low;
        }

        if (data[low] == value) {
            low++;
        }

        return low;
    }

    private static int makePrimeList(int n, int[] prime) {
        boolean[] composite = new boolean[n / 3];

        sieveOfEratosthenes(composite);

        boolean toggle = false;
        int p = 5, i = 0, j = 2;

        prime[0] = 2;
        prime[1] = 3;

        while (p <= n) {
            if (!composite[i++]) {
                prime[j++] = p;
            }
            // -- never mind, it's ok.
            p += (toggle = !toggle) ? 2 : 4;
        }

        return j; // number of primes
    }

    private static void sieveOfEratosthenes(final boolean[] composite) {
        int d1 = 8;
        int d2 = 8;
        int p1 = 3;
        int p2 = 7;
        int s1 = 7;
        int s2 = 3;
        int n = 0;
        int len = composite.length;
        boolean toggle = false;

        while (s1 < len) { // -- scan sieve
            if (!composite[n++]) { // -- if a prime is found, cancel its multiples
                int inc = p1 + p2;

                for (int k = s1; k < len; k += inc) {
                    composite[k] = true;
                }

                for (int k = s1 + s2; k < len; k += inc) {
                    composite[k] = true;
                }
            }

            if (toggle = !toggle) { // Never mind, it's ok.
                s1 += d2;
                d1 += 16;
                p1 += 2;
                p2 += 2;
                s2 = p2;
            } else {
                s1 += d1;
                d2 += 8;
                p1 += 2;
                p2 += 6;
                s2 = p1;
            }
        }
    }

    public static void main(String[] args) {
        int n = Integer.parseInt(args[0]);
        long nanos = System.nanoTime();
        BigInteger fact = factorial(n);
        nanos = System.nanoTime() - nanos;
        // Commented out because it takes ages to print
        //System.out.println(fact);
        System.out.println(nanos / 1e9);
    }
}