Как вычислить факториал большого числа?

Уважаемые программисты! Домучив-таки собственные StringToFloat и FloatToString, решил я написать самопальную программку для вычисления n!
Писал на Delphi, используя рекурсию (ногами не бейте, плиз), с целью потом выдрать код в дизассемблере и заюзать его в ассемблерной проге.
Столкнулся со следующей проблемой: используя тип LongInt я могу вычислить максимум 16!, Int64 - 20!, Double - 170!, а дальше - никак
Мне подкинули исходничек на C++, который вычисляет факториалы гораздо больших чисел (10000!, например), используя массив целых чисел в памяти, вот он:

Код (Text):

1. // factorial.cpp : Defines the entry point for the console application.
2. //
3.  
4. #include "stdafx.h"
5. #include <fstream.h>
6. #include <memory.h>
7. #include <stdlib.h>
8.  
9.  
10.  
11. const int MAX_SIZE=99999;//99999;
12.  
13. int sizeEl[3]={0,0,0};//0-numtmp,1- simple,2 - 10, 3- 100, 4- 1000
14.  
15.  
16. int addMass(int mass[],int masstmp[],int stepp,int size1,int size2,int nuli)
17. {
18.     if(stepp>0)
19.     {  
20.     int starti;
21.     if(size1>=(nuli+size2)){starti=size1;}else{starti=size2+nuli;};
22.     for(int j=1;j<=stepp;j++)
23. {//for 2 st                
24.     for(int k=nuli;k<=(starti);k++)
25.     {
26. mass[k]=mass[k]+masstmp[k-nuli];
27. if(mass[k]>9)
28. {//st if 1 
29.     mass[k+1]=mass[k+1]+mass[k]/10;mass[k]=mass[k]-(mass[k]/10*10);    
30.     if((k+1)>starti){starti=k+1;if(starti>MAX_SIZE-1){cout<<"Err. Toooo big param!!!"<<endl;cin>>starti;exit(1);return -1;}}
31. }//if 1 end
32.     }
33. }//for 2 end   
34.     return starti; 
35.     }
36.     else{return size1;}
37.    
38. }
39.  
40. int main(int argc, char* argv[])
41. {
42.  
43. cout<<"-----------------------"<<endl;
44. cout<<"|                     |"<<endl;
45. cout<<"|    wf.exe v2.0      |"<<endl;
46. cout<<"|                     |"<<endl;
47. cout<<"| make by:            |"<<endl;
48. cout<<"| _ZaliZo             |"<<endl;
49. cout<<"| chat32x64@yandex.ru |"<<endl;
50. cout<<"| in                  |"<<endl;
51. cout<<"| Ukraine             |"<<endl;
52. cout<<"| Crimea              |"<<endl;
53. cout<<"|                     |"<<endl;
54. cout<<"| xx.09.2006          |"<<endl;
55. cout<<"|  xx.xx.xx           |"<<endl;
56. cout<<"|                     |"<<endl;
57. cout<<"-----------------------"<<endl;
58. int fact;// факториал, который нужно найти(вводится пользователем)
59. int num[MAX_SIZE];//факториал в виде массива
60. int numtmp[MAX_SIZE];//временный массив для хранения промежуточного факториала
61. for(int pp=0;pp<MAX_SIZE;pp++)//цикл инициализации массивов
62. {num[pp]=0;numtmp[pp]=0;}
63. num[0]=1;
64. numtmp[0]=1;
65. cout<<"Enter factorial"<<endl;
66. cin>>fact;
67. cout<<"Start serch factorial of "<<fact<<" ..."<<endl;
68. for(int i=1;i<fact;i++)
69. {//for 1 st
70. cout<<i*100/fact<<"%"<<'\r';// вывод процента выполненых работ
71. for(int tm=0;tm<=sizeEl[1];tm++){numtmp[tm]=num[tm];sizeEl[0]=sizeEl[1];}
72. if(i==1){numtmp[0]=1;}
73. int step10000=i/10000;
74. int step1000=(i-step10000*10000)/1000;
75. int step100=(i-step10000*10000-step1000*1000)/100;
76. int step10=(i-step10000*10000-step1000*1000-step100*100)/10;
77. int step1=i-step10000*10000-step1000*1000-step100*100-step10*10;
78. sizeEl[1]=addMass(num,numtmp,step1,sizeEl[1],sizeEl[0],0);
79. sizeEl[1]=addMass(num,numtmp,step10,sizeEl[1],sizeEl[0],1);
80. sizeEl[1]=addMass(num,numtmp,step100,sizeEl[1],sizeEl[0],2);
81. sizeEl[1]=addMass(num,numtmp,step1000,sizeEl[1],sizeEl[0],3);
82. sizeEl[1]=addMass(num,numtmp,step10000,sizeEl[1],sizeEl[0],3);
83. }//for 1 end
84. if((sizeEl[1]==-1)||(sizeEl[4]==-1))
85. { cout<<"Size of element is too big, Sorry... ["<<sizeEl<<"]"<<endl;cin>>fact;}
86. else
87. {
88. cout<<endl<<endl;
89. cout<<"write to file"<<endl;
90. ofstream of("data.txt");
91. for(int t=sizeEl[1];t>=0;t--)//вывод часла на консоль/файл
92. {/*cout<<num[t];*/of<<num[t];}
93. of.close();
94. }
95. cout<<endl<<endl;
96. cout<<"End operation"<<endl;
97. cin>>fact;
98.     return 0;
99. }

Вот только яконкретно запутался, дизассемблируя его и не понял, как там всё делается Удручает и то, что стандартный виндовский калькулятор вычисляет те же значения раза в три быстрее...
Пробовал сотворить такую приблуду самостоятельно - столкнулся с необходимостью множить память на память, а mul и imul требуют операнд-регистр...
Помогите алгоритмом и советами по реализации, пожалуйста

 

Калькулятор не высчитывает все знаение факториала, а только первые 32 символа, а потом число знаков после запятой....
Там используются не только основные регистры процессора но и вероятнее всего ФПУ/ММХ

Множить память на память используя основные регистры и цикл - это прямая дорога к зависанию...
стоит ли говорить, что imul eax,eax,eax не прокатит....
Вторым множителем должна быть константа...
Остается только умножение циклом.... И вот тут.... если подсчитать, сколько нужно циклов молотить.... )))))))))
-----------------------------------------------------------------------------------

Вот, почитай тут
http://www.uni-vologda.ac.ru/students/it10/OpisanieGlav.html
Там довольно дзенно перечисленны комманды ммх, фпу, ссе, и несколько примеров.


Кстати, за ассемблер. В данном случае никаких переводов в строки и обратно не требуется. Здесь обычная работа с числами и массивом.

Если дойдут руки, днем сброшу код...

 

Такс... вот... лови...
В регистре st(1) ,будет результат факторизации.
Верхний предел 1500!
Компилятор масм

Делаешь батник main.bat, суешь этот код, потом запускаешь батник.
Получается ехе на выходе.

А там играйся сколько хочешь

Код (Text):

1. ;@Echo off
2. ;goto make
3. ;-----------------------------------------------------------------
4.  
5. .586p
6. .mmx
7. .model flat, stdcall
8. option casemap:none
9. .data
10.  
11. Result  dq 1
12. cm  dq 1
13. Num dq 0
14.  
15.  
16. .code
17. start:
18.     mov eax,1000
19.     call Fact
20.     ret          
21.  
22. ;==========================================================
23. ;       Fact
24. ;==========================================================
25. ; ВХОД:
26. ;   eax - число, которое нужно возвести в степень
27. ; ВЫХОД:
28. ;   FPU [ST(1)] - результат
29. ;==========================================================
30. Fact proc
31.     mov dword ptr [Num],eax
32.     finit
33.     fild qword ptr [cm]
34.     fild qword ptr [Result]
35.     fild qword ptr [Num]
36. @@:
37.     fmul st(1),st(0)    ; m = m*(m-1)
38.     fsub st(0),st(2)    ; n = n-1
39.     dec eax
40.     jnz @B
41.     ret
42. Fact endp
43. end start
44.  
45. ;-----------------------------------------------------------------
46. :make
47. \masm32\bin\ml /c /coff main.bat
48. \masm32\bin\Link /SUBSYSTEM:WINDOWS /SECTION:.text,EWR /OPT:NOREF main.obj
49. del main.obj
50. pause
51. exit

Если нужно перегнать значения регистров в память, наиболее просто заюзать fsave
Если нужно прогнать факториал под 10000 и более - достаточно просто начальное число разложить на множители, - чистая математика, так сказать... и юзай далее эту процедурку...

 

-

 

если на делфи пишешь посмотри в сторону использования модуля FGInt.

 

И вот еще что...
Рекурсия в данном случае просто НЕДОПУСТИМА
Потому как вызов процедуры + параметр + счетчик = 12 байт.
А т.к. ты хочешь юзать факториал больших чисел, больше 10000, то тут прямая дорога к переполнению стека....

 

W4FhLF

Нифига...
2000! не тянет. вся программа падает с криками о переполнении стека.

 

По поводу исходников на си - сразу видно, писал извращенец...
Смысл там простой...

Когда результат достигает поределенного значения, происходит "сокращение"
Т.е. подсчитывается число нулей, загоняется в n, и далее весь результат делится на 10^n
И происходит сложение m = m+n.

Результат представляется в виде
Result = a * (10^m);
При чем сначала преобразовывается inttostr(a), а потом inttostr(m)
Смысл такой. Хотя там как-то криво он реализован...

Некое округление, так сказать.... При чем чем больше m, тем менее точный результат.
И если такой элементарный код автор растянул на столько строк, то в топку такой листинг.

 

использовать длинную арифметику, как в №1 на си, хоть и довольно долго, зато надежно и число ограничивается только объемом допустимой памяти

 

FreeManCPM

А так же продолжительностью жизни автора

 

был когда то у меня старенький АМД под 1 гигагерцег. и читал я тогда МК. и там проскочила стотья про вычисление факториала.. исход был поскальный. я его загнал в делфи. вобщем 3000! около пол-минуты.

 

возьми либу GMP - и прилинкуй статически к проге на дельфях - или в виде ДЛЛ (еще проще, если со статикой есть проблемы) - и делай факториал вплоть до MAX_ULONG - то есть, до 2^32 - 1 - там функция есть специально для этого. Если тебе мало MAX_ULONG - делай в той же GMP но руками (тогда длина будет ограничена только размерами памяти). Хотя если тебе нужны аппроксимированные значения, можно гамма функцию вычислить (по крайней мере порядок результата будешь знать)

 

nitrotoluol

Серьёзно? Не повезло тебе значит, а у меня вот 10000! с помощью этого модуля меньше чем за секунду вычисляется

PS Athlon 3500+

 

- не серьёзно я когда-то давно примерно то что nitrotoluol в #8 написал, реализовал на программируемом МК (если ничего не путаю ~15 операций в секнду) потестировал на "разумных" числах - Ок - зарядил "хорошее" число подождал часа ~2,5 и принялся считать а сколько же нужно времени - убедился, что замечание:

очень верное

 

Всмысле мне не веришь или процессоры с частотой пару гигагерц в диковинку?

 

W4FhLF
Y_Mur

Вы говорите о разных вещах.

Факториал фкториалу рознь. Есть алгоритмы упрощенного нахождения порядка факториала, и алгоритмы нахождения факториала с определенной стпепенью погрешности.

Так вот... В первом случае даже мой код #3 несложно переделать на нахождение порядка факториала 10000!. если изначально разбить число на множители. Но толку от такого факториала - ноль. Согласитесь, кому нужно значение с погрешностью в 70%. А именно столько составляет неточность при упрощенном нахождении порядка...
Время затраченное на такой расчет приблизительно около 2х минут.

Чем выше точность - тем выше время, необходимое для расчета.
При погрешности, стремящейся к 0, время нахождения чудовищно... Еще есть табличные методы нахождения, т.е. известны некоторые значения, от которых в последствии и отталкиваются, а не находят их... Так что спор без выяснения предмета спора - думаю не уместен

 

nitrotoluol

Не совсем понимаю, что ты имеешь ввиду под порядком факториала, что-то вроде такого:
10000! = 2.8462596809170545189064132121199e+35659
?

Если да, то с помощью FGInt я находил именно сам факториал с точностью вплоть до 40000 разряда, дальше идут нули.

 

W4FhLF
Нет в смысле на 2ГГц это неинтересно - толи дело МК ~15 оп/сек )) (это типа шутка была)
nitrotoluol
Может не совсем внимательно глянул твой код, но я тогда в МК я просто разделил мантиссу и порядок на два МК регистра (там практически такой же регистровый стек как в FPU). Когда мантисса зашкаливала за заданный предел - делил её на 1000, прибавлял к регистру порядка +3, когда этого переставало хватать увеличивал делитель, правда не уверен что МК тогда за 2.5 часа всё таки дошёл до этого чудного момента
Ну и под конец приводил регистры с мантиссой порядком в удобочитаемый вид
Накапливающаяся погрешность в этом случае конечно есть но есно далеко но ~70%

 

W4FhLF
Y_Mur
Хм.... интересно.... значит у меня пробелы со знаниями..... сейчас поэксперементирую.... ))

 

Имхо при больших значениях аргумента лучше воспользоваться Гамма-функцией Эйлера
x! = Г(x+1); Апроксимирующую формулу для Г(х) даю. На асме и дельфи.
Просто вставь в код соответсвующие секции и все должно заработать.