[Из песочницы] Два способа быстрого вычисления факториала

Понятие факториала известно всем. Это функция, вычисляющая произведение последовательных натуральных чисел от 1 до N включительно: N! = 1×2 * 3 *… * N. Факториал — быстрорастущая функция, уже для небольших значений N значение N! имеет много значащих цифр.Попробуем реализовать эту функцию на языке программирования. Очевидно, нам понадобиться язык, поддерживающий длинную арифметику. Я воспользуюсь C#, но с таким же успехом можно взять Java или Python.

Итак, простейшая реализация (назовем ее наивной) получается прямо из определения факториала:

static BigInteger FactNaive (int n) { BigInteger r = 1; for (int i = 2; i <= n; ++i) r *= i; return r; } На моей машине эта реализация работает примерно 1,7 секунд для N=50000.Далее я покажу два алгоритма, которые работают более чем в два раза быстрее наивной реализации.Первый алгоритм основан на том соображении, что длинные числа примерно одинаковой длины умножать эффективнее, чем длинное число умножать на короткое (как в наивной реализации). То есть нам нужно добиться, чтобы при вычислении факториала множители постоянно были примерно одинаковой длины.

Пусть нам нужно найти произведение последовательных чисел от L до R, обозначим его как P (L, R). Разделим интервал от L до R пополам и посчитаем P (L, R) как P (L, M) * P (M + 1, R), где M находится посередине между L и R, M = (L + R) / 2. Заметим, что множители будут примерно одинаковой длины. Аналогично разобьем P (L, M) и P (M + 1, R). Будем производить эту операцию, пока в каждом интервале останется не более двух множителей. Очевидно, что P (L, R) = L, если L и R равны, и P (L, R) = L * R, если L и R отличаются на единицу. Чтобы найти N! нужно посчитать P (2, N).

Посмотрим, как будет работать наш алгоритм для N=10, найдем P (2, 10):

P (2, 10)P (2, 6) * P (7, 10)(P (2, 4) * P (5, 6)) * (P (7, 8) * P (9, 10))((P (2, 3) * P (4)) * P (5, 6)) * (P (7, 8) * P (9, 10))(((2×3) * (4)) * (5×6)) * ((7×8) * (9×10))((6×4) * 30) * (56×90)(24×30) * (5 040)720×5 0403 628 800

Получается своеобразное дерево, где множители находятся в узлах, а результат получается в корнеДерево вычисления факториала

Реализуем описанный алгоритм:

static BigInteger ProdTree (int l, int r) { if (l > r) return 1; if (l == r) return l; if (r — l == 1) return (BigInteger)l * r; int m = (l + r) / 2; return ProdTree (l, m) * ProdTree (m + 1, r); } static BigInteger FactTree (int n) { if (n < 0) return 0; if (n == 0) return 1; if (n == 1 || n == 2) return n; return ProdTree(2, n); } Для N=50000 факториал вычисляется за 0,8 секунд, что вдвое быстрее, чем в наивной реализации.Второй способ быстрого вычисления использует разложение факториала на простые множители. Очевидно, что в разложении N! участвуют только простые множители от 2 до N. Попробуем посчитать, сколько раз простой множитель K содержится в N!, то есть узнаем степень множителя K в разложении. Каждый K-ый член произведения 1 * 2 * 3 *… * N увеличивает показатель на единицу, то есть показатель степени будет равен N / K. Но каждый K2-ый член увеличивает степень еще на единицу, то есть показатель становится N / K + N / K2. Аналогично для K3, K4 и так далее. В итоге получим, что показатель степени при простом множителе K будет равен N / K + N / K2 + N / K3 + N / K4 +…

Для наглядности посчитаем, сколько раз двойка содержится в 10! Двойку дает каждый второй множитель (2, 4, 6, 8 и 10), всего таких множителей 10 / 2 = 5. Каждый четвертый дает четверку (22), всего таких множителей 10 / 4 = 2 (4 и 8). Каждый восьмой дает восьмерку (23), такой множитель всего один 10 / 8 = 1 (8). Шестнадцать (24) и более уже не дает ни один множитель, значит, подсчет можно завершать. Суммируя, получим, что показатель степени при двойке в разложении 10! на простые множители будет равен 10 / 2 + 10 / 4 + 10 / 8 = 5 + 2 + 1 = 8.

Если действовать таким же образом, можно найти показатели при 3, 5 и 7 в разложении 10!, после чего остается только вычислить значение произведения:

10! = 28×34 * 52×71 = 3 628 800

Осталось найти простые числа от 2 до N, для этого можно использовать решето Эратосфена:

static BigInteger FactFactor (int n) { if (n < 0) return 0; if (n == 0) return 1; if (n == 1 || n == 2) return n; bool[] u = new bool[n + 1]; // маркеры для решета Эратосфена List> p = new List>(); // множители и их показатели степеней for (int i = 2; i <= n; ++i) if (!u[i]) // если i - очередное простое число { // считаем показатель степени в разложении int k = n / i; int c = 0; while (k > 0) { c += k; k /= i; } // запоминаем множитель и его показатель степени p.Add (new Tuple(i, c)); // просеиваем составные числа через решето int j = 2; while (i * j <= n) { u[i * j] = true; ++j; } } // вычисляем факториал BigInteger r = 1; for (int i = p.Count() - 1; i >= 0; --i) r *= BigInteger.Pow (p[i].Item1, p[i].Item2); return r; } Эта реализация также тратит примерно 0,8 секунд на вычисление 50000!

© Habrahabr.ru