33,008
社区成员
发帖
与我相关
我的任务
分享擂台：计算阶乘的 18 位非零尾数
阶乘不用解释，阶乘的 18 位非零尾数需要简单解释一下。阶乘的十进制结果中末尾都会有很多个数字 0 (5以下的除外)，去掉这些连续的 0 之后最后 18 位数字即是这里所说的 18 位非零尾数。之所以称这为“非零”，是因为最后一个数字肯定不是 0。举两个例子，10! = 3,628,800，其 18 位非零尾数为 36,288；24! = 620,448,401,733,239,439,360,000，其 18 位非零尾数为 044,840,173,323,943,936
目前我写的最新版本，对于 10,000 位的整数，在我的机器上(P4 1.8A + 256M, 超频到 2.4)可以在 1 秒内计算其阶乘的 18 位非零尾数。现在对于我使用的算法来讲，速度已经算是达到极至了。故来此设立擂台。
程序中虽然还有可以再进行优化的地方，但都不值得去做了，因为对速度的提高根本起不了多大的作用。
相关程序和更多的介绍可以访问 http://www.ysiw.info/T18/T18.html
欢迎高手打擂，欢迎有兴趣的人参与讨论。
能给出更好的算法者，我可以再送 2000 分。
...全文
请发表友善的回复…
发表回复
liangbch 2005-09-01
- 打赏
- 举报
对于n!，如果不要求计算精确值(实际上，在实际应用中，大多不需要精确值)，用stirling公式依然可以可以求出 很大的 n(如10^10000)阶乘。
mathe 2005-09-01
- 打赏
- 举报
怎么变成算n!的帖子了. 计算n!其实没有多少意义,最到计算到(10^10)!就非常不错了.
这里是计算n!最后几个非零的数,其中n可以非常非常巨大,比如10^10000等
这里是计算n!最后几个非零的数,其中n可以非常非常巨大,比如10^10000等
oo 2005-08-23
- 打赏
- 举报
mark
liangbch 2005-08-18
- 打赏
- 举报
如果你对计算阶乘感兴趣，请看我的两个帖子：
http://61.186.252.131/Expert/topic/2267/2267097.xml?temp=.4246332
http://61.186.252.131/Expert/topic/2427/2427925.xml?temp=.3127558
和gxq的两个帖子和一篇文章：
http://community.csdn.net/Expert/topic/3441/3441530.xml
http://community.csdn.net/Expert/topic/3197/3197332.xml
http://blog.csdn.net/gxqcn/archive/2004/06/22/18810.aspx
目前我看到的一个计算大数阶乘最快的程序是刘楚雄写的一个程序，可从http://download.enet.com.cn/html/033432005081102.html 下载，
其核心代码中FFT来自http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html,
使用FFT 或者 FNT 的大数乘法运算是目前理论界最快的算法了。
我也准备利用前面的FFT代码编写一个大数阶乘程序，看看能否比刘楚雄更快。
最后,希望我们共同进步。
http://61.186.252.131/Expert/topic/2267/2267097.xml?temp=.4246332
http://61.186.252.131/Expert/topic/2427/2427925.xml?temp=.3127558
和gxq的两个帖子和一篇文章：
http://community.csdn.net/Expert/topic/3441/3441530.xml
http://community.csdn.net/Expert/topic/3197/3197332.xml
http://blog.csdn.net/gxqcn/archive/2004/06/22/18810.aspx
目前我看到的一个计算大数阶乘最快的程序是刘楚雄写的一个程序，可从http://download.enet.com.cn/html/033432005081102.html 下载，
其核心代码中FFT来自http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html,
使用FFT 或者 FNT 的大数乘法运算是目前理论界最快的算法了。
我也准备利用前面的FFT代码编写一个大数阶乘程序，看看能否比刘楚雄更快。
最后,希望我们共同进步。
junmayang 2005-08-17
- 打赏
- 举报
楼上的有道理，我的程序还很粗略，在内存分配上欠细，开始的时候固定死了，所以n会有上限，这是第一个需要改进的。空间的占用没有考虑过多，这个不是关键，考虑优化一下，你的评估的方法值得借鉴。运算速度只有用汇编代码优化可能会快一点，不过提升的幅度应该很有限。改进以后再发代码，呵呵
总的来说代码上的优化的空间已经不大了，要获得性能的大幅提高，只能从数学方面了，因该优化的空间很大。
总的来说代码上的优化的空间已经不大了，要获得性能的大幅提高，只能从数学方面了，因该优化的空间很大。
liangbch 2005-08-17
- 打赏
- 举报
to junmayang(笨猪):
再说说占用的存储空间问题：你的程序每个__int64 表示12位十进制数，即每个byte表示1.5位数。我的程序每个DWORD表示9位数据，即byte表示2.25位数，从这方面说，你的方案不是最优。
再说说占用的存储空间问题：你的程序每个__int64 表示12位十进制数，即每个byte表示1.5位数。我的程序每个DWORD表示9位数据，即byte表示2.25位数，从这方面说，你的方案不是最优。
liangbch 2005-08-17
- 打赏
- 举报
这个程序在一开始便分配一块固定大小的内存，按当前的值，当n的值达到2042268，将导致写内存错。
下面给出我编的一个程序，应该比楼上的更简单。因为核心算法使用汇编代码，所以速度比楼上的更快。这个程序的另一个特点是，在内存的使用上按需分配，对n值没有硬性限制。
#include "stdio.h"
#include "windows.h"
#include "math.h"
#define RADIX 1000000000
#define PI 3.1415926535897932384626433832795
#define E 2.7182818284590452353602874713527
#ifndef DWORD
typedef unsigned long DWORD;
#endif
DWORD rad=RADIX;
int calcResultLen(DWORD n) //ÀûÓÃstirling¹«Ê½¼ÆËã n!µÄ³¤¶È
{
if (n<=1)
return 1;
return log10(2*PI*(double)n)/2 + n * log10((double)n/E)+1;
}
DWORD calcCore(DWORD *p,DWORD n,DWORD c)
{
_asm
{
push ebx
;
mov ebx,p
xor edx,edx
mov eax,[ebx]
mul n
add eax,c
adc edx,0
div rad
mov [ebx],edx
;
pop ebx
}
}
void Num2Str(DWORD n, char *strEnd)
// Êýת»¯Îª´®£¬²»ÓÃsprintf µÄÀíÓÉÊÇΪÁËÌá¸ßËÙ¶È
{
char *p;
int i;
for ( i=0,p=strEnd;i<9;i++)
{
*p-- =( n % 10 +'0');
n /=10;
}
}
DWORD* calc( DWORD n, DWORD *pBuff,int buffLen)
{
DWORD *p; //lowest unit
DWORD *pEnd=pBuff+buffLen-1; //highest unit
DWORD c; //c carry number
//---------------
*pEnd=1;
while (n>1)
{
for (c=0,p=pBuff+buffLen-1;p>=pEnd;p--)
c=calcCore(p,n,c);
if (c>0)
{
pEnd--; *pEnd=c;
}
n--;
}
return pEnd;
}
int main(int argc, char* argv[])
{
DWORD n,buffLen,time;
DWORD *pEnd,*p;
DWORD *pBuff=NULL;
char *pRes=NULL;
char *pChar;
FILE *fp=NULL;
char fileName[16];
//------------------------
while (true)
{
printf("please In put a number, 0 exit:\nn=?");
scanf("%d",&n);
if (n<1)
goto End;
sprintf(fileName,"%d!.txt",n);
buffLen=(calcResultLen(n)+8) / 9;
pBuff= new DWORD[buffLen];
pRes= new char[buffLen * 9];
fp=fopen(fileName,"wb");
if (pBuff==NULL || pRes==NULL || fp==NULL)
goto Exit;
//---------------------
time=GetTickCount();
pEnd=calc(n,pBuff,buffLen);
time= GetTickCount()-time;
printf("CAL TIME:%d ms\n",time);
//-----------------------------
time=GetTickCount();
for (p=pEnd,pChar=pRes+8;p<pBuff+buffLen;p++)
{
Num2Str(*p,pChar);
pChar+=9;
}
pChar=pRes;
while ( *pChar=='0' ) pChar++; //ÂÔ¹ýÇ°µ¼Áã
fwrite(pChar,1,pRes+buffLen*9-pChar,fp); //Ò»´ÎдÈëÔËËã½á¹û
time= GetTickCount()-time;
printf("OUT TIME:%d ms\n\n",time);
Exit:
if ( pBuff!=NULL)
{ delete[] pBuff; pBuff=NULL; }
if (pRes!=NULL)
{ delete[] pRes; pRes=NULL; }
if (fp!=NULL)
{ fclose(fp); fp=NULL; }
}
End:
return 0;
}
下面给出我编的一个程序，应该比楼上的更简单。因为核心算法使用汇编代码，所以速度比楼上的更快。这个程序的另一个特点是，在内存的使用上按需分配，对n值没有硬性限制。
#include "stdio.h"
#include "windows.h"
#include "math.h"
#define RADIX 1000000000
#define PI 3.1415926535897932384626433832795
#define E 2.7182818284590452353602874713527
#ifndef DWORD
typedef unsigned long DWORD;
#endif
DWORD rad=RADIX;
int calcResultLen(DWORD n) //ÀûÓÃstirling¹«Ê½¼ÆËã n!µÄ³¤¶È
{
if (n<=1)
return 1;
return log10(2*PI*(double)n)/2 + n * log10((double)n/E)+1;
}
DWORD calcCore(DWORD *p,DWORD n,DWORD c)
{
_asm
{
push ebx
;
mov ebx,p
xor edx,edx
mov eax,[ebx]
mul n
add eax,c
adc edx,0
div rad
mov [ebx],edx
;
pop ebx
}
}
void Num2Str(DWORD n, char *strEnd)
// Êýת»¯Îª´®£¬²»ÓÃsprintf µÄÀíÓÉÊÇΪÁËÌá¸ßËÙ¶È
{
char *p;
int i;
for ( i=0,p=strEnd;i<9;i++)
{
*p-- =( n % 10 +'0');
n /=10;
}
}
DWORD* calc( DWORD n, DWORD *pBuff,int buffLen)
{
DWORD *p; //lowest unit
DWORD *pEnd=pBuff+buffLen-1; //highest unit
DWORD c; //c carry number
//---------------
*pEnd=1;
while (n>1)
{
for (c=0,p=pBuff+buffLen-1;p>=pEnd;p--)
c=calcCore(p,n,c);
if (c>0)
{
pEnd--; *pEnd=c;
}
n--;
}
return pEnd;
}
int main(int argc, char* argv[])
{
DWORD n,buffLen,time;
DWORD *pEnd,*p;
DWORD *pBuff=NULL;
char *pRes=NULL;
char *pChar;
FILE *fp=NULL;
char fileName[16];
//------------------------
while (true)
{
printf("please In put a number, 0 exit:\nn=?");
scanf("%d",&n);
if (n<1)
goto End;
sprintf(fileName,"%d!.txt",n);
buffLen=(calcResultLen(n)+8) / 9;
pBuff= new DWORD[buffLen];
pRes= new char[buffLen * 9];
fp=fopen(fileName,"wb");
if (pBuff==NULL || pRes==NULL || fp==NULL)
goto Exit;
//---------------------
time=GetTickCount();
pEnd=calc(n,pBuff,buffLen);
time= GetTickCount()-time;
printf("CAL TIME:%d ms\n",time);
//-----------------------------
time=GetTickCount();
for (p=pEnd,pChar=pRes+8;p<pBuff+buffLen;p++)
{
Num2Str(*p,pChar);
pChar+=9;
}
pChar=pRes;
while ( *pChar=='0' ) pChar++; //ÂÔ¹ýÇ°µ¼Áã
fwrite(pChar,1,pRes+buffLen*9-pChar,fp); //Ò»´ÎдÈëÔËËã½á¹û
time= GetTickCount()-time;
printf("OUT TIME:%d ms\n\n",time);
Exit:
if ( pBuff!=NULL)
{ delete[] pBuff; pBuff=NULL; }
if (pRes!=NULL)
{ delete[] pRes; pRes=NULL; }
if (fp!=NULL)
{ fclose(fp); fp=NULL; }
}
End:
return 0;
}
WoodJohn 2005-08-17
- 打赏
- 举报
真高人也
junmayang 2005-08-17
- 打赏
- 举报
大部分的代码处理输入和输出了，呵呵
junmayang 2005-08-17
- 打赏
- 举报
复杂？用于计算的代码才二十多行，我觉得存储结构还是可取的，占的内存也很少
liangbch 2005-08-17
- 打赏
- 举报
代码写得太复杂了。
junmayang 2005-08-17
- 打赏
- 举报
用于计算的代码实际很少，基本思路：采用10^12进制。
junmayang 2005-08-17
- 打赏
- 举报
发布我的初级算法，没有做过数学上的优化，完全是代码计算，P42.8G，10000！耗时1s。
在VC60编译通过，水平有限，程序比较拙劣，不要笑我啊，呵呵
------------------------------------------------------------------------------
#include "stdio.h"
#include "stdlib.h"
#include "memory.h"
#include "string.h"
#include <windows.h>
#pragma warning(disable : 4035)
inline unsigned __int64 GetCycleCount(void)
{
_asm _emit 0x0F
_asm _emit 0x31
}
class KTimer
{
unsigned __int64 m_startcycle;
public:
unsigned __int64 m_overhead;
KTimer(void)
{
m_overhead = 0;
Start();
m_overhead = Stop();
}
void Start(void)
{
m_startcycle = GetCycleCount();
}
unsigned __int64 Stop(void)
{
return GetCycleCount() - m_startcycle - m_overhead;
}
};
#define FORWARD 1000000000000
void main () {
int hir_length = 0;
int zero_num = 0;
int i = 0;
int j = 0;
int k = 0;
int N;
char x[100];
char *s = x;
char flag;
double f;
int H_INT_SIZE = 1000000;
__int64 *hir_value;
__int64 *p;
__int64 temp = 0;
__int64 temp1 = 0;
__int64 tm1, tm2;
__int64 cpuFre;
KTimer timer;
timer.Start();
printf("Calculate sytem frequency...\n");
Sleep(1000);
tm1 = timer.Stop();
cpuFre = tm1;
printf("CPU FREQUENCY : %I64dHZ\n\n", cpuFre);
hir_value = (__int64 *)malloc(H_INT_SIZE * sizeof(__int64));
while (1) {
printf("please In put a number, 0 exit:\nn=?");
gets(x);
i = 0;
s = x;
while (*s) {
s++;
i++;
}
s = x;
N = 0;
if (*s == '0') {
if (*(s + 1) != 0) {
continue;
}
else {
return;
}
}
while (i) {
i--;
j = (*s - 48);
if (j > 9 || j < 0) {
N = -1;
break;
}
for (k = 0; k < i; k++) {
j *= 10;
}
N += j;
s++;
}
if (N ==0) {
return;
}
else if (N == -1) {
continue;
}
memset(hir_value, 0, H_INT_SIZE * sizeof(__int64));
p = hir_value;
*p = 1;
hir_length = 1;
zero_num = 0;
timer.Start();
for (i = 1; i <= N; i++) {
p = hir_value + zero_num;
while (!*p) {
p++;
zero_num++;
}
temp = 0;
for(j = zero_num; j < hir_length; j++, p++) {
temp1 = hir_value[j] * i + temp;
temp = temp1 / FORWARD;
*p = temp1 - temp * FORWARD;
}
if (temp) {
*p = temp;
hir_length++;
}
}
tm1 = timer.Stop();
timer.Start();
temp = hir_value[hir_length - 1];
int len_r = 1;
while (temp >= 10) {
temp = temp / 10;
len_r++;
}
len_r += 12 * (hir_length - 1);
FILE *fp = fopen(strcat(x, "!.txt"), "w");
flag = 1;
fprintf(fp, " %12I64d", hir_value[hir_length - 1]);
p = hir_value + hir_length - 2;
for (i = hir_length - 2; i >= 0; i--, p--) {
if (flag == 5) {
fprintf(fp, "\n ");
flag = 0;
}
else {
fprintf(fp, " ");
}
flag++;
temp = *p;
j = 1;
while (temp >= 10) {
j++;
temp = temp / 10;
}
for (k = 0; k < 12 - j; k++) {
fprintf(fp, "0");
}
fprintf(fp, "%I64d", *p);
}
fprintf(fp, "\n");
fprintf(fp, "\nLength : %d\n", len_r);
tm2 = timer.Stop();
f = ((double)tm1) / cpuFre;
fprintf(fp, "\nCAL TIME : %fs\n", f);
printf("cal time : %fs\n", f);
f = ((double)tm2) / cpuFre;
fprintf(fp, "OUT TIME : %fs\n", f);
printf("out time : %fs\n", f);
fclose(fp);
}
if (hir_value) {
free(hir_value);
hir_value = NULL;
}
return;
}
在VC60编译通过，水平有限，程序比较拙劣，不要笑我啊，呵呵
------------------------------------------------------------------------------
#include "stdio.h"
#include "stdlib.h"
#include "memory.h"
#include "string.h"
#include <windows.h>
#pragma warning(disable : 4035)
inline unsigned __int64 GetCycleCount(void)
{
_asm _emit 0x0F
_asm _emit 0x31
}
class KTimer
{
unsigned __int64 m_startcycle;
public:
unsigned __int64 m_overhead;
KTimer(void)
{
m_overhead = 0;
Start();
m_overhead = Stop();
}
void Start(void)
{
m_startcycle = GetCycleCount();
}
unsigned __int64 Stop(void)
{
return GetCycleCount() - m_startcycle - m_overhead;
}
};
#define FORWARD 1000000000000
void main () {
int hir_length = 0;
int zero_num = 0;
int i = 0;
int j = 0;
int k = 0;
int N;
char x[100];
char *s = x;
char flag;
double f;
int H_INT_SIZE = 1000000;
__int64 *hir_value;
__int64 *p;
__int64 temp = 0;
__int64 temp1 = 0;
__int64 tm1, tm2;
__int64 cpuFre;
KTimer timer;
timer.Start();
printf("Calculate sytem frequency...\n");
Sleep(1000);
tm1 = timer.Stop();
cpuFre = tm1;
printf("CPU FREQUENCY : %I64dHZ\n\n", cpuFre);
hir_value = (__int64 *)malloc(H_INT_SIZE * sizeof(__int64));
while (1) {
printf("please In put a number, 0 exit:\nn=?");
gets(x);
i = 0;
s = x;
while (*s) {
s++;
i++;
}
s = x;
N = 0;
if (*s == '0') {
if (*(s + 1) != 0) {
continue;
}
else {
return;
}
}
while (i) {
i--;
j = (*s - 48);
if (j > 9 || j < 0) {
N = -1;
break;
}
for (k = 0; k < i; k++) {
j *= 10;
}
N += j;
s++;
}
if (N ==0) {
return;
}
else if (N == -1) {
continue;
}
memset(hir_value, 0, H_INT_SIZE * sizeof(__int64));
p = hir_value;
*p = 1;
hir_length = 1;
zero_num = 0;
timer.Start();
for (i = 1; i <= N; i++) {
p = hir_value + zero_num;
while (!*p) {
p++;
zero_num++;
}
temp = 0;
for(j = zero_num; j < hir_length; j++, p++) {
temp1 = hir_value[j] * i + temp;
temp = temp1 / FORWARD;
*p = temp1 - temp * FORWARD;
}
if (temp) {
*p = temp;
hir_length++;
}
}
tm1 = timer.Stop();
timer.Start();
temp = hir_value[hir_length - 1];
int len_r = 1;
while (temp >= 10) {
temp = temp / 10;
len_r++;
}
len_r += 12 * (hir_length - 1);
FILE *fp = fopen(strcat(x, "!.txt"), "w");
flag = 1;
fprintf(fp, " %12I64d", hir_value[hir_length - 1]);
p = hir_value + hir_length - 2;
for (i = hir_length - 2; i >= 0; i--, p--) {
if (flag == 5) {
fprintf(fp, "\n ");
flag = 0;
}
else {
fprintf(fp, " ");
}
flag++;
temp = *p;
j = 1;
while (temp >= 10) {
j++;
temp = temp / 10;
}
for (k = 0; k < 12 - j; k++) {
fprintf(fp, "0");
}
fprintf(fp, "%I64d", *p);
}
fprintf(fp, "\n");
fprintf(fp, "\nLength : %d\n", len_r);
tm2 = timer.Stop();
f = ((double)tm1) / cpuFre;
fprintf(fp, "\nCAL TIME : %fs\n", f);
printf("cal time : %fs\n", f);
f = ((double)tm2) / cpuFre;
fprintf(fp, "OUT TIME : %fs\n", f);
printf("out time : %fs\n", f);
fclose(fp);
}
if (hir_value) {
free(hir_value);
hir_value = NULL;
}
return;
}
wuyi8808 2005-08-15
- 打赏
- 举报
向 mathe() 学习.
junmayang 2005-08-15
- 打赏
- 举报
我的程序确实是初级算法，本人数学水平有限，数学上没有经过优化，只是代码上初步进行了优化。
不过我的程序理论上没有上限，只要内存足够大。
不过我的程序理论上没有上限，只要内存足够大。
liangbch 2005-08-11
- 打赏
- 举报
to junmayang(笨猪):
我在贴子(http://community.csdn.net/Expert/topic/4154/4154011.xml)发表的一个执行文件仅201字节的程序就可以赶上你的程序。该程序可计算1-31999的阶乘，运行结果输出到txt文件，详情见该贴。
我在贴子(http://community.csdn.net/Expert/topic/4154/4154011.xml)发表的一个执行文件仅201字节的程序就可以赶上你的程序。该程序可计算1-31999的阶乘，运行结果输出到txt文件，详情见该贴。
liangbch 2005-08-11
- 打赏
- 举报
to junmayang(笨猪) ;
计算1万的阶乘用时1秒，这个速度并不快，仅相当于我的初级算法。我的高级算法和gxq的程序可比这个你的程序快上百倍。
计算1万的阶乘用时1秒，这个速度并不快，仅相当于我的初级算法。我的高级算法和gxq的程序可比这个你的程序快上百倍。
junmayang 2005-08-11
- 打赏
- 举报
楼主，我的程序把10000的阶乘全算出来才要1秒
gxqcn 2005-07-11
- 打赏
- 举报
为防止CSDN“私吞”帖子，我已将本贴及本版4106353号贴收录在了我的个人主页中：
http://maths.diy.myrice.com/florilegium/
为保持原滋原味，收录时是从原贴直接按保存得到的，未作任何删改编辑。
（本贴对应网页为：http://maths.diy.myrice.com/florilegium/csdn/FactorialTail.htm）
另，在我发布的 HugeCalc V5.0.0.1 中，已包含该问题的相关文档、源代码及相关程序，
下载页为：http://maths.myrice.com/software.htm#02
http://maths.diy.myrice.com/florilegium/
为保持原滋原味，收录时是从原贴直接按保存得到的，未作任何删改编辑。
（本贴对应网页为：http://maths.diy.myrice.com/florilegium/csdn/FactorialTail.htm）
另，在我发布的 HugeCalc V5.0.0.1 中，已包含该问题的相关文档、源代码及相关程序，
下载页为：http://maths.myrice.com/software.htm#02
gxqcn 2005-07-11
- 打赏
- 举报
to xstring(麻雀)：
楼主的 http://www.ysiw.info/T18/T18.html 近期怎么无法访问了？
另外，楼主程序当 N ≡ a (mod 10^18) , 1<=a<=8 运行错误的 bug 解掉否？
楼主的 http://www.ysiw.info/T18/T18.html 近期怎么无法访问了？
另外，楼主程序当 N ≡ a (mod 10^18) , 1<=a<=8 运行错误的 bug 解掉否？
加载更多回复（59）
