Предлагаю рассмотреть «асимптотически наиболее быстрый» алгоритм вычисления факториала!
Для начала напомню, что факториал n — это произведение всех натуральных чисел от 1 до n (), при этом ;
1. Декомпозиция факториала
Введём функцию, именуемую swinging factorial, следующим образом:
Данная дробь всегда будет целым числом по простой причине — она кратна центральному биномиальному коэффициенту , который равен
Разворачивая определение swinging factorial, мы получим новую рекуррентную формулу факториала:
Она будет особенно хороша, если мы научимся эффективно вычислять значения .
2. Простые множители swinging factorial
Обозначим как степень простого числа в примарном разложении . Тогда будет справедлива следующая формула:
Для последнего выражения воспользуемся тем фактом, что , и получим нужную нам формулу.
Как следствие, и . Если нечётно, то . Другие частные случаи:
здесь означает количество единиц в двоичном представлении числа . Все эти факты могут быть использованы для дополнительной оптимизации в коде. Доказательства я приводить не буду, при желании вы легко сможете получить их самостоятельно.
Теперь, зная степени всех простых делителей , у нас есть способ вычисления swinging factorial:
3. Трудоёмкость алгоритма
Можно показать, что вычисление имеет сложность . Как ни странно, вычисление имеет ту же сложность (в оценке используется алгоритм Шёнхаге-Штрассена, отсюда и такая интересная трудоёмкость; доказательства по ссылке в конце статьи).
Несмотря на то, что формально перемножение чисел от 1 до n имеет ту же трудоёмкость, алгоритм PrimeSwing на практике оказывается самым быстрым.
UPDATE: как было замечено в этом комментарии, тут я ошибся, перемножение чисел от 1 до n имеет большую трудоёмкость.
Ссылки и реализация
- страница с различными алгоритмами вычисления факториала;
- детальное описание алгоритма из статьи (и не только).
// 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;
}
Комментарии (21)
choupa
29.04.2017 11:31+1Мне кажется, самый практичный метод вычисления факториала
array[0, 1, 2, 6, 24, 120, ...]
Простите за серость, а зачем нужно вычислять большие факториалы? Криптография? Комбинаторика? Ряды Тейлора? Собеседование в Гугл?
ibessonov
29.04.2017 11:35Боюсь, что это вопрос не ко мне. Вероятнее всего различные алгоритмы вычисления факториала имеют соревновательный характер.
Хранить массив может быть очень затратно по памяти, особенно если нужны значения порядка 1000000! (вдруг они кому-то нужны). На практике я такого не встречал и сам всегда хранил кэшированные значения в обычном массиве (во время решения задач по программированию).
alexeykuzmin0
03.05.2017 11:14Факториалы чисел порядка миллионов иногда встречаются в олимпиадном программировании в задачах на комбинаторику, где после вывода формулы решения оказывается, что нужно считать большие биномиальные коэффициенты по какому-нибудь простому модулю. Но даже в этом случае проще сначала посчитать все прямые и обратные факториалы, а потом уже из них собирать то, что нужно, более сложный алгоритм вычисления факториалов мне лично ни разу не пригодился.
mikhaelkh
29.04.2017 11:46+1Несмотря на то, что формально перемножение чисел от 1 до n имеет ту же трудоёмкость, алгоритм PrimeSwing на практике оказывается самым быстрым.
Предположим, что разделим числа на две половины, рекурсивно вычислим произведение и перемножим половины. Сложность у этого метода ~M(n*log(n))*log(n), что хуже ~M(n*log(n)) у PrimeSwing. Или как предлагается перемножать?ibessonov
29.04.2017 12:12Извиняюсь, действительно у меня в этом моменте ошибка. Вычисление с помощью PrimeSwing имеет трудоёмкость как у простого перемножение двух чисел порядка n!, а не вычисления n! стандартным способом. Сейчас внесу обновление в статью.
Спасибо!
koldyr
29.04.2017 11:53Стоит отметить, что статья имеет практическую направленность, но согласно вики:
Метод Шёнхаге — Штрассена считался асимптотически быстрейшим методом с 1971 до 2007 годы, пока не был заявлен новый метод с лучшей оценкой сложности умножения.[3] На практике метод Шёнхаге — Штрассена начинает превосходить более ранние классические методы, такие как умножение Карацубы и алгоритм Тоома — Кука (обобщение метода Карацубы), начиная с целых чисел порядка 10^10000-10^40000.
Таким образом становится актуальным вопрос о константе в оценке.ibessonov
29.04.2017 12:23Интересно вышло, алгоритм PrimeSwing был опубликован в 2008-м, автор уже мог знать о более оптимальном способе умножения чисел.
Если я верно понял, автор производил вычисления факториалов примерно до 1000000, то есть результаты были не настолько большими, чтобы использовать Алгоритма Фюрера (если верить википедии). А так да, теоретическую оценку наверняка можно улучшить.
На практике я не заморачивался и использовал то умножение, какое было реализовано за меня в BigInteger.
Спасибо за замечание.alexeykuzmin0
03.05.2017 11:16В BigInteger (по крайней мере раньше было) реализовано умножение столбиком (за квадрат). И не забудьте еще расходы времени на перевод из десятичной в двоичную СС.
vlanko
04.05.2017 11:51В Java используются более быстрые алгоритмы. А перевод для вывода — если печатать в консоль — да, займет время :)
malishich
29.04.2017 16:39Стать интересная, не знал что и здесь простые числа применяются. На практике мне не встречались задачи где нужно было вычислять просто факториалы чисел больше 20-30 (например в uint64_t влезает 20, а вот 21! уже нет, зато в double влезает, хоть и с погрешностью). Встречались отношения (дроби) больших факториалов (комбинаторные сочетания Cmn для близких чисел и размещения Amn для сильно разных чисел), порядка 300-т. Но в этих задачах я решал просто разложением на простые множители всех чисел произведения в числителе и знаменателе, подсчёту количества таких множителей в числителе и знаменателе, сокращению, и дальнейшему вычислению отношения (там обычно много чего сокращалось и части дроби «влазили» с малой погрешностью в тип double, но чаще вообще в uint64_t).
MikailBag
29.04.2017 22:08Зачем считать для цешек факториалы, если есть треугольник паскаля?
Я могу оценить время на генерацию N строк как O(N^3).
В Вашем случае (N = 300) — пара минут максимум.alexeykuzmin0
03.05.2017 11:19Зачем использовать треугольник Паскаля, если есть замечательное соотношение C(n, k) = C(n — 1, k — 1) * (n / k)? Вычисляется, в отличие от треугольника Паскаля, за O(k) операций над числами, а из операций потребуется только умножение и деление на короткое число (причем деление всегда гарантированно нацело будет).
mikhaelkh
29.04.2017 21:41В этом алгоритме мы по сути выделяем большой множитель, являющийся квадратом. То есть A[k] = A[k+1]^2 * B[k+1], A[0] = n!, для вычисления A[k] мы вычисляем A[k+1], B[k+1]. Но мы можем выделить ещё больший множитель, оставляя в B[k+1] только произведение простых, входящие в A[k] нечётных степенях. Кажется, это должно работать ещё быстрее, хотя разница должна быть невелика.
nikitadanilov
30.04.2017 02:26«в примарном разложении» — должно быть «в разложении на простые множители». «Примарное разложение» это немного другое (см. теорему Ласкера).
koldyr
С какой сложностью нам достается это знание?
ibessonov
Простых чисел от 1 до n всего O(n / log n), вычисление lp происходит за O(log n).
Сложностью же перемножения простых чисел в соответствующих степенях будет O(n (log n)2 log log n) — это уже не очевидная оценка, её доказательство есть по второй ссылке.
koldyr
Асимптотику распределения простых чисел я знаю.
Я правильно понимаю, что придется хранить или генерировать таблицу простых чисел от 2 до n/2?
Какова трудоемкость детерминированного алгоритма теста на простоту?
ibessonov
Не совсем, понадобится таблица простых чисел от 2 до n.
Решето Эратосфена составляется за O(n * log log n), это меньше, чем время вычисления факториала. Отдельного способа проверки на простоту алгоритм не требует.
koldyr
Прекрасно. Как-то я упустил, что сложность решета ниже сложности вычисления факториала.