Существуют ли быстрые алгоритмы получения факториала?

[Dmitriy40]

Думаю для малых значений аргумента используют точное перемножение, а как сильно выходит за разрядную сетку, то используют формулу Муавра-Стирлинга, там если взять несколько членов разложения в ряд точность хорошая вроде бы.

[Eimrine]

Dmitriy40 в сообщении #1126540 писал(а):

как сильно выходит за разрядную сетку

?

[mihaild]

Dmitriy40 в сообщении #1126537 писал(а):

PARI/GP любые факториалы до $10^{16}$ вычисляет мгновенно, вот пример:

Кажется, это посчитано в вещественной арифметике. Может быть какая-нибудь оптимизация формулы Стирлинга (или даже непосредственно она)?

[Aritaborian]

arseniiv в сообщении #1126536 писал(а):

У Mathematica есть часть справки, посвящённая реализации. Может, там и про факториал написано.

Не написано. Гуглёж тоже не особо помог.

[mihaild]

Совсем на коленке написал глупую эвристику на порядок перемножения - на 100к обгоняет наивную реализацию в 30 раз (0.03с против 0.95), на 200к - в 50+ (0.07 против 4.1). На миллионе на моей машине отрабатывает за 0.6с, на 10 миллионах - за 10.6.

Используется синтаксис C++

mpz_class twoway_factorial(const size_t n) {
    std::vector<mpz_class> data;
    size_t r = 1;
    for (size_t i = 1; i <= n; ++i) {
        r *= i;
        if (r > (1L << 63) / n) {
            data.emplace_back(r);
            r = 1;
        }
    }
    data.emplace_back(r);
    while (data.size() > 1) {
        std::vector<mpz_class> new_data;
        new_data.reserve(data.size() / 2 + 1);
        for (size_t i = 0; i < data.size() / 2; ++i) {
            new_data.emplace_back(data[i] * data[data.size() - i - 1]);
        }
        if (data.size() % 2) {
            new_data.emplace_back(data[data.size() / 2]);
        }
        data = new_data;
    }
    return data[0];
}

mpz_fac_ui из gmp на 10 миллионах отрабатывает за 4 секунды, итерируясь по простым числам (https://gmplib.org/manual/Factorial-Algorithm.html).

[Eimrine]

mihaild в сообщении #1126558 писал(а):

Совсем
mpz_fac_ui из gmp на 10 миллионах отрабатывает за 4 секунды, итерируясь по простым числам (https://gmplib.org/manual/Factorial-Algorithm.html).

А остачу кода дадите, чтобы можно было запустить?

Я решил проверить свой захудалый лаптоп вот этим кодом http://www.luschny.de/math/factorial/
- факториала от ста тысяч я так и не дождался, 10000 выводит с заметной задержкой. Алгоритм позиционируется как самый быстрый (нашел по моей ссылке в начале поста), но реализован на явно не самом быстром языке (запускал на Python последней версии).

[arseniiv]

(Оффтоп)

Aritaborian в сообщении #1126557 писал(а):

Не написано. Гуглёж тоже не особо помог.

Эх, подвела!

[mihaild]

Eimrine в сообщении #1126575 писал(а):

А остачу кода дадите, чтобы можно было запустить?

Используется синтаксис C++

#include <gmpxx.h>
#include <iostream>
#include <ctime>
#include <vector>
#include <cassert>

mpz_class naive_factorial(size_t n) {
    mpz_class result(1);
    for (size_t i = 1; i <= n; ++i) {
        result *= i;
    }
    return result;
}

mpz_class twoway_factorial(const size_t n) {
    std::vector<mpz_class> data;
    size_t r = 1;
    for (size_t i = 1; i <= n; ++i) {
        r *= i;
        if (r > (1L << 63) / n) {
            data.emplace_back(r);
            r = 1;
        }
    }
    data.emplace_back(r);
    while (data.size() > 1) {
        std::vector<mpz_class> new_data;
        new_data.reserve(data.size() / 2 + 1);
        for (size_t i = 0; i < data.size() / 2; ++i) {
            new_data.emplace_back(data[i] * data[data.size() - i - 1]);
        }
        if (data.size() % 2) {
            new_data.emplace_back(data[data.size() / 2]);
        }
        data = new_data;
    }
    return data[0];
}

int main() {
    size_t n;
    std::cin >> n;
    clock_t t;
    t = clock();
    mpz_t qq;
    t = clock();
    mpz_init(qq);
    mpz_fac_ui(qq, n);
    std::cout << "gmp: " << static_cast<double>(clock() - t) / CLOCKS_PER_SEC << std::endl;
    auto e = twoway_factorial(n);
    std::cout << "iterative: " << static_cast<double>(clock() - t) / CLOCKS_PER_SEC << std::endl;
    t = clock();
    auto q = naive_factorial(n);
    std::cout << "naive: " << static_cast<double>(clock() - t) / CLOCKS_PER_SEC << std::endl;
    assert(q == e);
    return 0;
}
 

(да, я в курсе, что так мерять производительность нельзя)

Eimrine в сообщении #1126575 писал(а):

факториала от ста тысяч я так и не дождался, 10000 выводит с заметной задержкой

У меня питоновский math.factorial на 100 тысячах за 0.25с отрабатывает.

[g______d]

Насколько я понимаю, алгоритмы умножения больших чисел наиболее эффективно работают, если множители одного и того же порядка. Поэтому можно оптимизировать, разбивая каждый раз на произведение наиболее близких чисел. Где именно выбирать точку разбиения, можно предварительно оценить с помощью формулы Стирлинга (за константу времени в вещественной арифметике).

А, ну кстати да, по простым числам должно быть еще быстрее, потому что степень вхождения вычисляется очень быстро по теореме Лежандра, а большие степени простых можно вычислять за логарифм умножений.

[sqribner48]

Большое спасибо всем. (Пока писал появилось несколько новых сообщений. Однако "не пропадать же добру".)
Не ожидал, что получится такая дискуссия. Я конечно искал в интернете и много нашел разных сайтов с обсуждением этого вопроса. Но меня конкретно интересовало какой алгоритм применяется в Mathematika. В справке доступной мне Mathematika ничего не нашел. Поэтому и обратился на сайт.
Мне кажется, самые лучшие варианты ответов есть на СоХабр (ссылка https://sohabr.net/habr/post/255761/ )
Там были рассмотрены несколько разных вариантов. Первый самый простой (Наивный алгоритм) в цикле. Второй - Алгоритм вычисления деревом, третий - Алгоритм вычисления факторизацией, четвертый - Библиотека GMP. Все подробно и детально расписано, приводятся листинги программ. Очень много есть модификаций алгоритмов в комментариях.

На сайте
http://ru.stackoverflow.com/questions/2479/%D0%A1%D0%B0%D0%BC%D1%8B%D0%B9-%D0%B1%D1%8B%D1%81%D1%82%D1%80%D1%8B%D0%B9-%D1%84%D0%B0%D0%BA%D1%82%D0%BE%D1%80%D0%B8%D0%B0%D0%BB
тоже есть много интересного и много ссылок.

mihaild в сообщении #1126549 писал(а):

Может быть какая-нибудь оптимизация формулы Стирлинга (или даже непосредственно она)?

На этом сайте один из вариантов как раз использует формулу Стирлинга. В результате получается формула

Цитата:

«... с бинарным масштабным коэффициентом, …. которая при правильно выбранной точности вычислений будет точной.»

Там же была указана очень интересная статья On the Complexity Calculating Factorials (по ссылке http://www.cecm.sfu.ca/personal/pborwein/PAPERS/P29.pdf) . В ней предлагают считать факториал за $O(\log \log n m(n \log n))$ (где $M(n)$ — время перемножения двух n-значных чисел).

Eimrine в сообщении #1126538 писал(а):

Луддизм 21-го столетия - соревноваться с компилятором в оптимизации.

Не совсем понял. Причем тут Луддиты? Кажется они ломали машины.
А уважаемый rockclimber вроде ничего не ломал.
Еще раз большое спасибо за ссылки. Буду смотреть и постараюсь разобраться.

P.S. У меня даже появилась одна идейка по рационализации программы простого цикла. За счет нее количество итераций в цикле уменьшается в два раза. Когда полностью ее добью, (если все получится) то выложу для обозрения

[Eimrine]

(Оффтоп)

sqribner48 в сообщении #1126580 писал(а):

Eimrine в сообщении #1126538 писал(а):

Луддизм 21-го столетия - соревноваться с компилятором в оптимизации.

Не совсем понял. Причем тут Луддиты? Кажется они ломали машины.

Мало того, что любой уважающий себя компилятор и так умеет делать эту оптимизацию, так ещё и не каждый высокоуровневый ЯП разрешит такое сделать, и не каждая числовая структура данных умножится от битового сдвига (а не испортится). Задача математика - написать быстрый алгоритм, а думать, на чём будет он исполняться и как будут выглядеть числа на уровне битов это как бы задача другого плана. Если алгоритм со сдвигом будет реализовыван, например, для квантового компьютера, то такая инструкция станет большим тормозом, и нужно будет ещё догадаться деобфусцировать сдвиг до деления на два.

mihaild в сообщении #1126578 писал(а):

У меня питоновский math.factorial на 100 тысячах за 0.25с отрабатывает.

math.factorial считал 50000 секунд 15, но питоновская консоль дико тормозила при скроллировании результата. Точно не замерял, но они оба работают с похожей скоростью, и оба едят только одно ядро. Возможно, для результата имеет значение то, что у меня сейчас используется сильно много свопа, но было видно, что процесс упирался в скорость ядра. Мне кажется, на кластере, где очень много процессоров, или даже на видеокарте, обычное умножение имеет реальные шансы быть быстрее всех из-за потенциально хорошей паралелизуемости.

P.S. К сожалению, мой gcc не запустил ваш код, ругается, что mpz_class, e, q does not have a type и mpz_t, mpz_init, mpz_fac_ui, q, e was not declared in this scope.

[mihaild]

Eimrine в сообщении #1126593 писал(а):

Возможно, для результата имеет значение то, что у меня сейчас используется сильно много свопа

Ну если вычисления в свопе - то, естественно, всё очень медленно. Хотя памяти тут нужно сравнительно немного (в небольшую константу раз больше, чем размер ответа).

Eimrine в сообщении #1126593 писал(а):

обычное умножение имеет реальные шансы быть быстрее всех из-за потенциально простой паралелизуемости

В алгоритме с факторизацией тоже основное время уходит на возведение простых чисел в степень, и на перемножение результатов - что параллелится ничуть не хуже. И даже само по себе умножение длинных чисел уже прекрасно параллелится.

Eimrine в сообщении #1126593 писал(а):

сожалению, мой gcc не запустил ваш код

Что этот код gcc не скомпилировал - неудивительно, т.к. это код на плюсах.

g++ -std=c++11 -O2 -lgmpxx -lgmp factorial.cpp

(и убедитесь, что у вас libgmp есть)

[Aritaborian]

sqribner48 в сообщении #1126580 писал(а):

Но меня конкретно интересовало какой алгоритм применяется в Mathematika. В справке доступной мне Mathematika ничего не нашел.

В разделе Some Notes on Internal Implementation, который имел в виду arseniiv, действительно ничего не говорится о факториале. На будущее можете иметь в виду, что Wolfram Research вообще не особые любители открывать публике внутренности своих программ; отделываются общими словами; впрочем, имеют на это полное право.

[grizzly]

(Оффтоп)

Aritaborian в сообщении #1126602 писал(а):

впрочем, имеют на это полное право.

Я не юрист, но время от времени слышу, что отдельные участники того проекта выражают сомнения по этому поводу.

[sqribner48]

grizzly в сообщении #1126626 писал(а):

В разделе Some Notes on Internal Implementation , который имел в виду arseniiv, действительно ничего не говорится о факториале.

Я в этом разделе нашел только вот это:

Цитата:

Combinatorial Functions
- Most combinatorial functions use sparse caching and recursion.
- $n!$ uses an $O(\log(n) M(n))$ algorithm of Schönhage based on dynamic decomposition to prime powers.

Первую строку я перевел как

Цитата:

Большинство комбинаторных функций использует редкое кэширование и рекурсию.

А что такое редкое кеширование?
Со второй строкой вроде ясно.
$n!$ использует алгоритм Шёнхаге, основанный на динамическом разложении на степени простых чисел.
(Если я правильно перевёл). Об этом писали iifat и g______d в своих постах.