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

Понятие факториала известно всем. Это функция, вычисляющая произведение последовательных натуральных чисел от 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 040

3 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. Но каждый K²-ый член увеличивает степень еще на единицу, то есть показатель становится N / K + N / K². Аналогично для K³, K⁴ и так далее. В итоге получим, что показатель степени при простом множителе K будет равен N / K + N / K² + N / K³ + N / K⁴ +…

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

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

10! = 2⁸ * 3⁴ * 5² * 7¹ = 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<Tuple<int, int>> p = new List<Tuple<int, int>>(); // множители и их показатели степеней
    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<int, int>(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!

This entry passed through the Full-Text RSS service - if this is your content and you're reading it on someone else's site, please read the FAQ at http://ift.tt/jcXqJW.