:)
我重读这个帖子，期待对另外一个问题程序的改进
我期待这个题有一个新的突破
如果允许小数的话，问题就变得很简单了。
比如 n!=3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446 E1000000000000
那么n=94851898540.1396977006676482198750221184274755718919148371815155857816611145105279808242712122064124653944534478504904850688162913
1# medie2005
做新手任务，凑字数
这个题目真的很有意思!
1# medie2005
难度确实挺大
我曾发过帖子说2的幂能以任何数字开始
这个问题和那帖子上的问题类似
但更难做了
无心人 发表于 2008-12-4 11:38 http://bbs.emath.ac.cn/images/common/back.gif
Matlab?
62# 无心人
今天仔细考虑了这个题目。我认为无需使用大数库。无心人的程序中调用了大数库，而且用到对数，故速度较慢。其实这个题目虽然需要做高精度计算，但只用到小整数乘以大数的运算，完全可以自己实现。我刚刚写了一个程序，未用到任何大数库。这个程序从1开始，一直查找到2^32-1,在我的电脑上只需要3分钟左右。
说明：
1. 我在程序中设定的精度是4个DWORD,每个DWORD表示9位10进制数，故最大精度不超过36位，若需更高精度，则只需改动一处即可。
2. 目前的代码最多可计算到2^32-1，若需要扩大搜索范围，需要实现一个64bit数和大数相乘的算法。
3. 我的程序依次计算每个32位整数n,求n!的科学计数法表示和pi的差，若其差小于当前的最小值，则输出n和n!.
下面是我的程序的输出结果：
20!=2.432902008176640000*10^18
23!=2.5852016738884976640000*10^22
28!=3.04888344611713860501504000000*10^29
154!=3.0897696138473508879585646636*10^271
332!=3.1034430734498676144426322732*10^694
612!=3.1345284031431087019726722580*10^1441
2703!=3.1398094010536875630936179148846*10^8104
8945!=3.139993064736104560624185498*10^31464
17254!=3.14084642621024495446191444792*10^65612
28726!=3.14106463595708329690285209393642766*10^115595
50583!=3.1415379577347882738422653077918*10^215977
527600!=3.14155838836242816120923043886*10^2789957
780018!=3.1415890794236473559982729593497*10^4257193
1812538!=3.141589535267260763834817102767*10^10556211
5385973!=3.141591234607732809306850985*10^33915312
5573998!=3.14159238184880918674222498012222406*10^35182367
52604972!=3.1415925665166406146782547740*10^383318353
65630447!=3.141592602756036723755404377401952*10^484537182
187024036!=3.141592617297502030148476882134*10^1465820139
410390476!=3.141592622630656164649753699858562*10^3356543814
464736214!=3.14159264915921444226889681595*10^3826132373
688395641!=3.141592651962886290907283002629*10^5784962808
It took 188 secs
#include "stdlib.h"
#include "stdio.h"
#include "string.h"
#include "windows.h"
typedef __int64INT64;
typedef unsigned long DWORD;
#define DEC_RADIX 1000000000
#define PRECISION_DIV9 4 //10进制数表示的数字个数除以9
const DWORD PI_VALUE[]={
{3,141592653},
{31,415926535},
{314,159265358},
{3141,592653589},
{31415,926535897},
{314159,265358979},
{3141592,653589793},
{31415926,535897932},
{314159265,358979323},
};
double scales[]=
{
1.00,
10.00,
100.00,
1000.00,
10000.00,
100000.00,
1000000.00,
10000000.00,
100000000.00,
};
int Uint32BitCount(DWORD n)
{
_asm
{
bsr eax,n
inc eax
}
}
//******************************************************
DWORD DEC_Array_Mul_DWORD_ALU(DWORD *a, size_t len, DWORD num)
//******************************************************
{
static DWORD _s_TEN9=DEC_RADIX;
//---------------------------
_asm
{
push esi
mov ecx,len
mov esi,a
lea esi,
xor ecx,ecx //carry
test esi,3 //the a must be 4byte align
jnzthisExit
jmp cmp000
loop_start:
mov eax,
mul dword ptr num
add eax,ecx
adc edx,0
div dword ptr _s_TEN9
mov ,edx
mov ecx,eax
sub esi,4
cmp000:
cmp esi,a
jae loop_start
thisExit:
mov eax,ecx
pop esi
}
}
double check_diff(DWORD num[],int shift)
{
int t,carry;
DWORD tArr;
tArr=PI_VALUE;
tArr=PI_VALUE;
t=tArr-num;
if (t<0)
{
carry=1;
t+=DEC_RADIX;
}
else
carry=0;
tArr=t;
t=tArr-num-carry;
tArr=t;
if (t<0)
{
return 999;
}
else
{
double diff;
diff =tArr + tArr * 0.000000001;
diff /= scales;
return diff;
}
}
void print_number(DWORD n,DWORD res[],int len,double exp)
{
char buff;
char *p;
int i,s;
INT64 e;
p=buff;
p+=sprintf(p,"%d",res);
s=strlen(buff);
if (len>1)
{
for (i=1;i<len;i++)
p+=sprintf(p,"%09d",res[ i ]);
}
printf("\n%u!=",n);
for (i=0;i<strlen(buff);i++)
{
printf("%c",buff[ i ]);
if (i==0)
printf(".");
}
e=(INT64)(exp+(len-1)*9+s-1);
printf("*10^%I64d\n",e);
}
void search_n()
{
DWORD i,j;
DWORD carry;
int shift;
int bc,fac_len; // fac(i)的结果的数组的长度
double diff,min_diff;
double res_exp;
DWORD fac; //fac(i)的结果的数组
DWORD time;
time=GetTickCount();
min_diff=1;
res_exp=0;
fac=1;
fac_len=1;
i=1;
//while (i<DEC_RADIX)
while (i!=0)
{
carry=DEC_Array_Mul_DWORD_ALU(fac,fac_len,i);
if (carry>=DEC_RADIX)
{
for (j=PRECISION_DIV9-2;j>=2;j--)
fac=fac;
fac=carry / DEC_RADIX;
fac=carry % DEC_RADIX;
if ( fac_len+2 <= PRECISION_DIV9)
fac_len+=2;
else if (fac_len+1 ==PRECISION_DIV9)
{
fac_len=PRECISION_DIV9;
res_exp+=9;
}
else
res_exp+=18;
}
else if (carry>0)
{
for (j=PRECISION_DIV9-1;j>=1;j--)
fac=fac;
fac=carry;
if (fac_len<PRECISION_DIV9)
fac_len++;
else
res_exp+=9;
}
bc=Uint32BitCount(fac);
switch (bc)
{
case 2:
shift=0;
diff=check_diff(fac,shift);
if (diff<min_diff)
{
min_diff=diff;
if (fac_len>=2)
print_number(i,fac,fac_len,res_exp);
}
break;
case 5:
shift=1;
diff=check_diff(fac,shift);
if (diff<min_diff)
{
min_diff=diff;
if (fac_len>=2)
print_number(i,fac,fac_len,res_exp);
}
break;
case 9:
shift=2;
diff=check_diff(fac,shift);
if (diff<min_diff)
{
min_diff=diff;
if (fac_len>=2)
print_number(i,fac,fac_len,res_exp);
}
break;
case 12:
shift=3;
diff=check_diff(fac,shift);
if (diff<min_diff)
{
min_diff=diff;
if (fac_len>=2)
print_number(i,fac,fac_len,res_exp);
}
break;
case 15:
shift=4;
diff=check_diff(fac,shift);
if (diff<min_diff)
{
min_diff=diff;
if (fac_len>=2)
print_number(i,fac,fac_len,res_exp);
}
break;
case 19:
shift=5;
diff=check_diff(fac,shift);
if (diff<min_diff)
{
min_diff=diff;
if (fac_len>=2)
print_number(i,fac,fac_len,res_exp);
}
break;
case 22:
shift=6;
diff=check_diff(fac,shift);
if (diff<min_diff)
{
min_diff=diff;
if (fac_len>=2)
print_number(i,fac,fac_len,res_exp);
}
break;
case 25:
shift=7;
diff=check_diff(fac,shift);
if (diff<min_diff)
{
min_diff=diff;
if (fac_len>=2)
print_number(i,fac,fac_len,res_exp);
}
break;
case 29:
shift=8;
diff=check_diff(fac,shift);
if (diff<min_diff)
{
min_diff=diff;
if (fac_len>=2)
print_number(i,fac,fac_len,res_exp);
}
break;
default:
break;
}
i++;
}
time=(GetTickCount()-time);
printf("It took %d secs\n",time/1000);
}
int main(int argc, char* argv[])
{
search_n();
return 0;
}
刚才仔细看了一下无心人的代码，明白了其所用的算法。其主要算法如下，
从start开始，逐个计算每个整数n的阶乘的常用对数的小数部分，记作log(fac(n)), 和PI的常用对数的小数部分比较，符合条件则输出。
当n是1000的整数倍时，计算base=n, lf0=log10(fac(base))的小数部分，
若n不是base的整数倍时，则计算从base+1到n之间所有数的累乘积的自然对数的小数部分lf1, 则 n!的对数的小数部分为 lf= lf0+lf1。
无心人的算法主要受制于 c语言库中(或者CPU浮点指令)中的双精度数的对数函数精度的制约，很难提高精度。
受无心人的启发，我打算编写一个不使用求双精度对数和多重精度对数函数的程序。
主要特点有：
1. 可计算到10^15或者更高。
2. 可分段计算。
3. 当n小于10亿时，使用自己写的多重精度大数乘以32bit整数的子程序,计算n的阶乘，fac(n),以10^9进制表示.
当n>=10亿且n是100万的倍数时，使用斯特林公式的普通形式(而不是指数形式)计算 fac(n)，以10^9进制表示,记作fac_0
当n>10亿且n不是100万的倍数时，令n=base+d(base是100万的整数倍，d>0 且 d<base)，使用自己写的多重精度大数乘以64bit整数的子程序,计算base+1到n的累乘积，记作fac_1, 则n的阶乘fac(n)=fac0*fac1
我的算法和无心人算法的差别，我的记作B，无心人的记作A。
1. B使用斯特林公式的普通形式，而A使用斯特林公式的对数形式。
2. B不计算n的阶乘的常用对数，而是直接计算n的阶乘的值，以10亿进制的表示形式。
3. A的核心部分为计算64bit整数的的常用对数，使用C语言函数库，A受制于double型数的精度，无法提高精度。 B的和核心部分计算一个多精度数与64bit整数的乘积，可通过增加多精度数数的长度来提高精度。
4. 性能的PK. 主要看 多精度数与64bit整数的乘积 VS log10(double),那个算得更快。