Fast Factorial Computation - PrimeSwing

ibessonov April 29, 2017 at 00:03

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

$inline$ ;

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

$inline$ 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

$inline$ 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

$inline$ . 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

$inline$ 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;
}

Tags: