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 (), wherein ;
We introduce a function called swinging factorial as follows:
This fraction will always be an integer for a simple reason - it is a multiple of the central binomial coefficient which is equal
Expanding the definition of swinging factorial , we get a new recurrent factorial formula:
It will be especially good if we learn how to efficiently calculate values .
We denote as a prime power in primary decomposition . Then the following formula will be valid:
Consequently, and . If a odd then . Other special cases:
here means the number of units in binary representation of a number . 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, we have a way to calculate swinging factorial :
It can be shown that the calculation has difficulty . Oddly enough, computingit 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.
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 (), wherein ;
1. Factorial decomposition
We introduce a function called swinging factorial as follows:
This fraction will always be an integer for a simple reason - it is a multiple of the central binomial coefficient which is equal
Expanding the definition of swinging factorial , we get a new recurrent factorial formula:
It will be especially good if we learn how to efficiently calculate values .
2. Simple factors swinging factorial
We denote as a prime power in primary decomposition . Then the following formula will be valid:
Evidence
We use the Legendre theorem on prime factors of the factorial :
For the last expression, we use the fact that , and we get the formula we need.
For the last expression, we use the fact that , and we get the formula we need.
Consequently, and . If a odd then . Other special cases:
here means the number of units in binary representation of a number . 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, we have a way to calculate swinging factorial :
3. The complexity of the algorithm
It can be shown that the calculation has difficulty . Oddly enough, computingit 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
- page with various factorial calculation algorithms;
- A detailed description of the algorithm from the article (and not only).
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;
}