Наткнувшись недавно на эту статью, я понял, что редко упоминаются способы вычисления факториала, отличные от банального перемножения последовательных чисел. Нужно эту ситуацию исправить.
Предлагаю рассмотреть «асимптотически наиболее быстрый» алгоритм вычисления факториала!
Для начала напомню, что факториал n — это произведение всех натуральных чисел от 1 до n (
), при этом
;
Введём функцию, именуемую swinging factorial, следующим образом:
Данная дробь всегда будет целым числом по простой причине — она кратна центральному биномиальному коэффициенту
, который равен
Разворачивая определение swinging factorial, мы получим новую рекуррентную формулу факториала:
Она будет особенно хороша, если мы научимся эффективно вычислять значения
.
Обозначим
как степень простого числа
в примарном разложении
. Тогда будет справедлива следующая формула:
Как следствие,
и
. Если
нечётно, то
. Другие частные случаи:
здесь означает количество единиц в двоичном представлении числа
. Все эти факты могут быть использованы для дополнительной оптимизации в коде. Доказательства я приводить не буду, при желании вы легко сможете получить их самостоятельно.
Теперь, зная степени всех простых делителей
, у нас есть способ вычисления swinging factorial:
Можно показать, что вычисление
имеет сложность
. Как ни странно, вычисление
имеет ту же сложность (в оценке используется алгоритм Шёнхаге-Штрассена, отсюда и такая интересная трудоёмкость; доказательства по ссылке в конце статьи).
Несмотря на то, что формально перемножение чисел от 1 до n имеет ту же трудоёмкость, алгоритм PrimeSwing на практике оказывается самым быстрым.
UPDATE: как было замечено в этом комментарии, тут я ошибся, перемножение чисел от 1 до n имеет большую трудоёмкость.
Предлагаю рассмотреть «асимптотически наиболее быстрый» алгоритм вычисления факториала!
Для начала напомню, что факториал n — это произведение всех натуральных чисел от 1 до n (
1. Декомпозиция факториала
Введём функцию, именуемую swinging factorial, следующим образом:
Данная дробь всегда будет целым числом по простой причине — она кратна центральному биномиальному коэффициенту
Разворачивая определение swinging factorial, мы получим новую рекуррентную формулу факториала:
Она будет особенно хороша, если мы научимся эффективно вычислять значения
2. Простые множители swinging factorial
Обозначим
Доказательство
Воспользуемся теоремой Лежандра о простых множителях факториала:
Для последнего выражения воспользуемся тем фактом, что
, и получим нужную нам формулу.
Для последнего выражения воспользуемся тем фактом, что
Как следствие,
Теперь, зная степени всех простых делителей
3. Трудоёмкость алгоритма
Можно показать, что вычисление
Несмотря на то, что формально перемножение чисел от 1 до n имеет ту же трудоёмкость, алгоритм PrimeSwing на практике оказывается самым быстрым.
UPDATE: как было замечено в этом комментарии, тут я ошибся, перемножение чисел от 1 до n имеет большую трудоёмкость.
Ссылки и реализация
- страница с различными алгоритмами вычисления факториала;
- детальное описание алгоритма из статьи (и не только).
Реализация на Java
// 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;
}