Fast factorial calculation - PrimeSwing

Having stumbled upon this article recently, I realized that the methods of calculating factorial, which are different from the banal multiplication of consecutive numbers, are rarely mentioned. Need to fix this situation.
I propose to consider the "asymptotically fastest" factorial calculation algorithm!
To begin, let me remind you that 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 to

$\ 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 to calculate values efficiently.

$n \ wr$ .

2. Simple swinging factorial multipliers

Denote

$l_p (n \ wr)$ as a power of prime

$p$ in primary decomposition

$n \ wr$ . Then the following formula will be true:

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

Evidence

We use the Legendre theorem on factorial simple factors :

\ 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}

$\ 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 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 ones in the binary representation of the 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 the 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, the calculation

$n!$ it has the same complexity ( Schönhage-Strassen’s algorithm is used in the assessment, hence the interesting complexity; the proof is 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 turns out to be the fastest in practice.

UPDATE: as noted in this commentary , here I was mistaken, multiplying the numbers from 1 to n is more difficult.

Links 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<BigInteger> 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<BigInteger> 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; }

Source: https://habr.com/ru/post/323770/

All Articles

Fast factorial calculation - PrimeSwing

1. Factorial decomposition

2. Simple swinging factorial multipliers

3. The complexity of the algorithm

Links and implementation

More articles: