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

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