Fast Factorial Computation - PrimeSwing

    Having recently stumbled upon this article , I realized that ways of calculating factorial are rarely mentioned, other than the banal multiplication of consecutive numbers. It is necessary to correct this situation.
    I propose to consider the "asymptotically fastest" factorial algorithm!
    To begin with, I recall that the factorial n is the product of all natural numbers from 1 to n ($ n! = 1 \ cdot2 \ cdot \ ldots \ cdot n $), wherein $ 0! = 1 $;

    1. Factorial decomposition


    We introduce a function called swinging factorial as follows:

    $ n \ wr = \ frac {n!} {\ lfloor n / 2 \ rfloor! ^ 2} $


    This fraction will always be an integer for a simple reason - it is a multiple of the central binomial coefficient $ \ binom {n} {\ lfloor n / 2 \ rfloor} $which is equal

    $ \ binom {n} {\ lfloor n / 2 \ rfloor} = \ frac {n!} {\ lfloor n / 2 \ rfloor! \ cdot \ lceil n / 2 \ rceil!} $


    Expanding the definition of swinging factorial , we get a new recurrent factorial formula:

    $ n! = \ lfloor n / 2 \ rfloor! ^ 2 \ cdot n \ wr $


    It will be especially good if we learn how to efficiently calculate values $ n \ wr $.

    2. Simple factors swinging factorial


    We denote $ l_p (n \ wr) $ as a prime power $ p $ in primary decomposition $ n \ wr $. Then the following formula will be valid:

    $ l_p (n \ wr) = \ sum_ {k \ geqslant1} \ lfloor \ frac {n} {p ^ k} \ rfloor \: mod \: 2 $


    Evidence
    We use the Legendre theorem on prime factors of the factorial :

    $ \ begin {array} {ccl} l_p (n! / \ lfloor n / 2 \ rfloor! ^ 2) & = & l_p (n!) - 2l_p (\ lfloor n / 2 \ rfloor!) \\ & = & \ sum_ {k \ geqslant1} \ lfloor n / p ^ k \ rfloor-2 \ sum_ {k \ geqslant1} \ lfloor \ lfloor n / 2 \ rfloor / p ^ k \ rfloor \\ & = & \ sum_ {k \ geqslant1 } (\ lfloor n / p ^ k \ rfloor-2 \ lfloor \ lfloor n / p ^ k \ rfloor / 2 \ rfloor) \\\ end {array} $


    For the last expression, we use the fact that $ j-2 \ lfloor j / 2 \ rfloor = j \: mod \: 2 $, and we get the formula we need.

    Consequently, $ l_p (n \ wr) \ leqslant log_p (n) $ and $ p ^ {l_p (n \ wr)} \ leqslant n $. If a$ p $ odd then $ l_p (p ^ a \ wr) = a $. Other special cases:

    $$ display $$ \ begin {array} {lrrl} (a) & \ lfloor n / 2 \ rfloor & <p \ leqslant & n & \ Rightarrow & l_p (n \ wr) = 1 \\ (b) & \ lfloor n / 3 \ rfloor & <p \ leqslant & \ lfloor n / 2 \ rfloor & \ Rightarrow & l_p (n \ wr) = 0 \\ (c) & \ sqrt {n} & <p \ leqslant & \ lfloor n / 3 \ rfloor & \ Rightarrow & l_p (n \ wr) = \ lfloor n / p \ rfloor \: mod \: 2 \\ (d) & 2 & <p \ leqslant & \ sqrt {n} & \ Rightarrow & l_p (n \ wr) <log_2 (n) \\ (e) & & p = & 2 & \ Rightarrow & l_p (n \ wr) = \ sigma_2 (\ lfloor n / 2 \ rfloor) \\ \ end {array } $$ display $$


    $ \ sigma_2 (n) $ here means the number of units in binary representation of a number $ n $. All these facts can be used for additional optimization in the code. I will not give evidence, if you wish, you can easily get it yourself.

    Now, knowing the degrees of all prime divisors$ n \ wr $, we have a way to calculate swinging factorial :

    $ n \ wr = \ prod_ {p \ leqslant n} p ^ {l_p (n \ wr)} $


    3. The complexity of the algorithm


    It can be shown that the calculation $ n \ wr $ has difficulty $ O (n (log \: n) ^ 2log \: log \: n) $. Oddly enough, computing$ n!  $it has the same complexity ( the Schoenhage-Strassen algorithm is used in the assessment , hence such an interesting laboriousness; evidence by reference at the end of the article).

    Despite the fact that formally multiplying numbers from 1 to n has the same complexity, the PrimeSwing algorithm in practice is the fastest.

    UPDATE: as was noted in this comment , here I was mistaken, multiplying numbers from 1 to n has a lot of complexity.

    References and implementation



    Java implementation
    // main function
    public static BigInteger factorial(int n) {
        return factorial(n, primes(n));
    }
    // recursive function with shared primes array
    private static BigInteger factorial(int n, int[] primes) {
        if (n < 2) return BigInteger.ONE;
        BigInteger f = factorial(n / 2, primes);
        BigInteger ps = primeSwing(n, primes);
        return f.multiply(f).multiply(ps);
    }
    // swinging factorial function
    private static BigInteger primeSwing(int n, int[] primes) {
        List multipliers = new ArrayList<>();
        for (int i = 0; i < primes.length && primes[i] <= n; i++) {
            int prime = primes[i];
            BigInteger bigPrime = BigInteger.valueOf(prime);
            BigInteger p = BigInteger.ONE;
            int q = n;
            while (q != 0) {
                q = q / prime;
                if (q % 2 == 1) {
                    p = p.multiply(bigPrime);
                }
            }
            if (!p.equals(BigInteger.ONE)) {
                multipliers.add(p);
            }
        }
        return product(multipliers, 0, multipliers.size() - 1);
    }
    // fast product for the list of numbers
    private static BigInteger product(List multipliers, int i, int j) {
        if (i > j) return BigInteger.ONE;
        if (i == j) return multipliers.get(i);
        int k = (i + j) >>> 1;
        return product(multipliers, i, k).multiply(product(multipliers, k + 1, j));
    }
    // Eratosthenes sieve
    private static int[] primes(int upTo) {
        upTo++;
        if (upTo >= 0 && upTo < 3) {
            return new int[]{};
        }
        int length = upTo >>> 1;
        boolean sieve_bool[] = new boolean[length];
        for (int i = 1, iterations = (int) Math.sqrt(length - 1); i < iterations; i++) {
            if (!sieve_bool[i]) {
                for (int step = 2 * i + 1, j = i * (step + 1); j < length; j += step) {
                    sieve_bool[j] = true;
                }
            }
        }
        int not_primes = 0;
        for (boolean not_prime : sieve_bool) {
            if (not_prime) not_primes++;
        }
        int sieve_int[] = new int[length - not_primes];
        sieve_int[0] = 2;
        for (int i = 1, j = 1; i < length; i++) {
            if (!sieve_bool[i]) {
                sieve_int[j++] = 2 * i + 1;
            }
        }
        return sieve_int;
    }
    


    Also popular now: