Tuesday, February 12, 2008

Как посчитать очень большой факториал?

Не так давно обратил внимание на то, как быстро Maple считает факториалы(встроенный функционал). 500000! - огромное число - он посчитал за пару секунд. Стало интересно - как же это реализовано? В итоге было написано несколько тестовых кусочков кода на ЯП, которые поддерживают большие числа. Первым таким языком был Erlang. Код весьма прост:
[UPD]
fact3(X)->fact3(X,1).

fact3(0,Acc)->Acc;
fact3(X,Acc)->fact3(X-1,Acc*X).
[/UPD]
Erlang соптимизирует рекурсию в цикл и в итоге этот код будет достаточно оптимальным. [UPD]На сравнительно небольших числах он работает отменно. Например, 50000! он посчитал за 6 секунд(особой точности здесь не требуется, поэтому я не старался получить максимально точные данные или тестировать в абсолютно чистой среде). Но 100000! у меня он вообще считать отказался - упал по нехватке памяти.[/UPD]
Вторым языком, который был опробован, стал PHP. Наивный код на PHP на больших числах не стесняясь выдавал INF, но тов. Josser(не знаю, куда поставить ссылку) довёл код до ума:

$targ = 50000;
$res = 1;

$startime = time();
for ($i=1; $i<=$targ; $i++) {
echo $i."\r\n";
$res = bcmul($res, $i);
}
$endtime = time();

echo $res."\r\n";

echo $endtime-$startime."\r\n";

Также Josser и протестировал его:
30000! - 99 секунды
50000! - 242 секунды

Как видим, работает оно медленней, чем вариант на Erlang, но хотя бы работает на 50000!.
Следующим языком стал Python. Код также прост, как и в предыдущих случаях:


def fact(x):
s=1
for i in range(2,x+1):
s*=i
return s

Протестировав, получаем:
30000! - 1.64100003242
50000! - 7.54700016975
100000! - 40.75

Как видим, Python смотрится просто героем на фоне Erlang и PHP. Однако, с увеличением чисел понимаем, что подобраться к заветному 500000! таким образом не получится.
Следующим протестированным языком стал Haskell. Я не особо знаком с этим языком, поэтому не могу похвастать и абсолютно правильным кодом. Главное - работает:

import Text.Printf
import Control.Exception
import System.CPUTime

time :: IO t -> IO t
time a = do
start <- getCPUTime
v <- a
end <- getCPUTime
let diff = (fromIntegral (end - start)) / (10^12)
printf "Computation time: %0.3f sec\n" (diff :: Double)
return v

fac n = if n == 0 then 1 else n * fac (n-1)

main = do
putStrLn "Starting..."
time $ fac 100000 `seq` return ()
putStrLn "Done."

Результаты тестирования:
30000! - 4.828 sec
50000! - 14.531 sec
100000! - 70.625 sec
Получше, чем Erlang и PHP, но хуже, чем Python. Также я вспомнил свои познания из области математики и попробовал организовать вычисление факториала через гамма-функцию. Код для вычисления логарифма гамма-функции:

cof :: [Double]
cof = [76.18009172947146,-86.50532032941677,24.01409824083091,-1.231739572450155,0.001208650973866179,-0.000005395239384953]

ser :: Double
ser = 1.000000000190015

gammaln :: Double -> Double
gammaln xx = let tmp' = (xx+5.5) - (xx+0.5)*log(xx+5.5)
ser' = foldl (+) ser $ map (\(y,c) -> c/(xx+y)) $ zip [1..] cof
in -tmp' + log(2.5066282746310005 * ser' / xx) where

Используется потом так:
exp (gammaln 20001)
Так считается 20000! Но, к сожалению, здесь Haskell прямо ответил - Infinity.

В итоге, после того, как мы опробовали кучу языков программирования, можно сделать вывод, что ни один из них изначально не подходит для решения заветной задачи - нахождения 500000!. Но раз стандартные средства не работают, время поискать нестандартное решение. И оно сразу находится - библиотека GNU MP или GMP. Кому хочется посмотреть, на что способна эта библиотека - попробуйте здесь. Также можно попробовать скачать и установить - даже без тонкой настройки всё равно получаем очень неплохие результаты.
Ну, и финальный аккорд - нахождение искомого 500000! при помощи GNU MP:

The result of executing 500000! is:

computation took 1418 ms
result is about 2632341 digits, not printing it


Итоговый результат: хотите посчитать очень большой факториал? Используйте GNU MP![UPD](который и используется в Хаскеле(?) )[/UPD]

[UPD]
Спасибо lionet, который заметил чушь, которую я назвал хвостовой рекурсией раньше, чем я ее исправил. Также спасибо yorool-gui за развитие Хаскелевского варианта. Я не буду исправлять текущий код - вы можете посмотреть на более оптимальные варианты решения задачи здесь.
[/UPD][UPD]
Кстати, в GMP используется функция mpz_fac_ui для того, чтобы подсчитать факториал. Т.е., это отдельная функция. И именно она даёт такие хорошие результаты. Тогда становится понятным, почему Haskell c GMP не даёт таких результатов, как просто использование данной функции.[/UPD]

18 comments:

  1. Комментарий почему-то пропал (ошибка или премодерация - непонятно). Попробую еще раз: ответ вот здесь - http://yorool-gui.livejournal.com/175800.html

    ReplyDelete
  2. Будешь смеяться, но Haskell и использует библиотеку GMP.

    ReplyDelete
  3. Вариант с Эрлангом левый - в этом случае здесь НЕТ хвостовой рекурсии. У меня вариант для 50000! делается за 30 секунд.

    Вариант с хвостовой рекурсией был бы таков:

    fac(N) -> fac(N, 1).
    fac(0, Acc) -> Acc;
    fac(N, Acc) -> fac(N - 1, Acc * N).

    ReplyDelete
  4. 2lionet:
    Действительно, с Erlang весьма нехорошо получилось - написал ерунду, теперь вот с утра разгребать. Правильный вариант у меня:
    fact3(X)->fact3(X,1).
    fact3(0,Acc)->Acc;
    fact3(X,Acc)->fact3(X-1,Acc*X).

    У меня считает за 6 секунд 50000!. Но 100000! падает по нехватке памяти.

    2yorool_gui: спасибо за комментарии по Хаскелю. Язык новый, неведомый, а так как я умудряюсь ошибиться и в более-менее знакомом Эрланге(и не успеваю исправить, пока не заметили:) ), то можно не удивляться проблемам с Хаскелем.

    2lionet:Также любопытно узнать, раз Хаскель использует GMP, то каким образом получаются такие разительные результаты?

    P.S. За все исправления спасибо - сейчас исправлю.

    ReplyDelete
  5. А где исходный код GMP-вычислялки?

    ReplyDelete
  6. Присоединяюсь к общемус сумасшедсвию.

    Код на плюсах который за 12 секунд посчитал 100000!, и за 325 секунд 500000!. +черти сколько времени на питоне конвертации из HEX в DEC.
    http://dobrokot.nm.ru/tempora/factorial.rar

    использовался Intel C++ Compiler, так как у него лучше всего с 64-битной арифметикой из того, что у меня под рукой.

    ReplyDelete
  7. Стало интересно - как же это реализовано? В итоге было написано несколько тестовых кусочков кода на ЯП, которые поддерживают большие числа.

    Дело в том, что это сравнение попы с пальцем, вы уж простите меня за прямоту аналогии. Maple быстро считает факториалы в первую очередь потому что он быстро перемножает большие числа (применяет что-нибудь типа FFT). И именно это должно быть вам интересно, а не алгоритм вычисления факториала, который сводится к простому циклу
    Когда же вы пишите подобный код, то кроме всего прочего вы в результатах измерения видите различные издержки связанные с ленивостью, вызовом функций и т.д.

    ReplyDelete
  8. Кстати в википедии есть кое-что про быстрое вычисление факториалов.

    ReplyDelete
  9. >Maple быстро считает факториалы в первую очередь потому что он быстро перемножает большие числа (применяет что-нибудь типа FFT). И именно это должно быть вам интересно, а не алгоритм вычисления факториала, который сводится к простому циклу

    Если вы не заметили, алгоритм вычисления факториала здесь и не рассматривается. Здесь рассматриваются скоростные характеристики и возможности тех алгоритмов работы с большими числами, которые заложены в ЯП. Т.е., я не даю описаний алгоритмов, я лишь предлагаю свои практические выкладки, ибо теоретических сведений действительно можно найти множество.

    ReplyDelete
  10. Если вы не заметили, то я привел соображения, по которым Ваш тест не может рассматриваться, как тест библиотек для вычислений с большими числами: а именно обилие вызовов влечёт за собой соотвественные накладные расходы, такие расходы в языках работающих на VM могут быть совершенно не приемлемыми, например, в случае когда стэк виртуальной машины размещается не в реальном стэке, а в динамически выделяемой памяти, и оптимизация хвостовой рекурсии не проводится.

    Кроме-того нигде не дана постановка задачи, а потому и возникает впечатление, что вы меряете непонятно что.

    ReplyDelete
  11. Постановка задачи здесь кроется в небольшом введении. Цель данной статьи - это отыскать способ подсчёта больших факториалов. Ну, и как разумный программист, поиск я начинаю с уже существующих велосипедов.

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

    В конце концов, приведённые числа(время подсчёта) не являются окончательными - они лишь призваны показать, насколько применим тот или иной "простой" подход к решению поставленной задачи.

    ReplyDelete
  12. Ну вы сами себе противоречите =)
    В первом своём ответе мне вы сказали, что алгоритм вычисления факториала здесь не рассматривается, теперь говорите, что цель данной статьи - это отыскать способ подсчёта больших факториалов.

    Вообщем это я так, занудствую =)

    ReplyDelete
  13. А попробуйте такой вариант:

    factorial n = product [ product [x, (x+k)..n] | x <- [1..k] ] where k = 14

    На моём стареньком Celeron-1.7 эта функция вычисляет 100 000! за 3.5 сек, а 200 000! за 10.4 сек, тогда как простой вариант product [1..n] соответственно за 27.3 и 117.0 сек... :о) т.е. от 8 до 11 раз быстрее ;о)

    Варьируя k можно подбирать наибольшее быстродействие. У меня оптимум: k=12..16

    ReplyDelete
  14. 2Prokhorov:
    Очень интересный вариант. При k=155 получаем 100000! за 1.125 на моей машине и 300000! за 10.266 сек. Если вы не против, я добавлю ваш вариант в статью как наилучший с указанием авторства?

    ReplyDelete
  15. очень странно, что эрланг отожрал память, с аккумулятором и хвостовой - не должен
    наблюдаю профайлером за f(100000) - не ест
    считал 300 секунд

    итого
    f(50000) - 62 секунды (6 у тебя)
    f(100000) - 292 секунды (значит примерно 30 должно быть у тебя)

    на всякий случай:
    -module(f).
    -export([f/1]).

    f(X) -> f(X, 1).

    f(0, Acc) -> Acc;
    f(X, Acc) -> f(X-1, Acc * X).

    ReplyDelete
  16. @Phil Pirozhkov
    А что у Вас за ОСь?

    ReplyDelete
  17. А попробуйте эту программу на "Питоне"

    def fact(n,e=None):
    if e is None:
    e=2**(int(math.log(n)/math.log(2)))
    if n<2:
    return 1
    elif e==0:
    return (n*(n-1))
    elif e>0:
    return (fact(n,(e/2))) * (fact((n-(e*2)),(e/2)))
    500000! вычисляет за минуту.
    Подобная программа на "Haskell" 1 000 000! вычисляет за 8 секунд.

    ReplyDelete
  18. На 3 Питоне:

    import math
    math.factorial(500000)

    За 22 секунды...

    ReplyDelete