阶乘和圆周率

无心人 · 发表于 2010-8-23 20:56:07

我重读这个帖子，期待对另外一个问题程序的改进

wayne · 发表于 2010-8-24 09:25:04

我期待这个题有一个新的突破

云梦 · 发表于 2012-10-31 08:09:30

如果允许小数的话，问题就变得很简单了。比如 n!=3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446 E1000000000000 那么n=94851898540.1396977006676482198750221184274755718919148371815155857816611145105279808242712122064124653944534478504904850688162913

clockman · 发表于 2012-12-1 20:22:20

1# medie2005 做新手任务，凑字数

郭先抢 · 发表于 2012-12-3 18:16:10

这个题目真的很有意思!

clockman · 发表于 2012-12-4 16:16:05

1# medie2005 难度确实挺大

gua64 · 发表于 2012-12-6 18:08:52

我曾发过帖子说2的幂能以任何数字开始这个问题和那帖子上的问题类似但更难做了无心人发表于 2008-12-4 11:38

Matlab?

liangbch · 发表于 2012-12-12 14:56:54

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

复制代码

liangbch · 发表于 2012-12-12 14:59:02

#include "stdlib.h"
#include "stdio.h"
#include "string.h"
#include "windows.h"
typedef __int64 INT64;
typedef unsigned long DWORD;
#define DEC_RADIX 1000000000
#define PRECISION_DIV9 4 //10进制数表示的数字个数除以9
const DWORD PI_VALUE[][2]={
{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,[esi+ecx*4-4]
xor ecx,ecx //carry
test esi,3 //the a must be 4byte align
jnz thisExit
jmp cmp000
loop_start:
mov eax,[esi]
mul dword ptr num
add eax,ecx
adc edx,0
div dword ptr _s_TEN9
mov [esi],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[2];
tArr[0]=PI_VALUE[shift][0];
tArr[1]=PI_VALUE[shift][1];
t=tArr[1]-num[1];
if (t<0)
{
carry=1;
t+=DEC_RADIX;
}
else
carry=0;
tArr[1]=t;
t=tArr[0]-num[0]-carry;
tArr[0]=t;
if (t<0)
{
return 999;
}
else
{
double diff;
diff =tArr[0] + tArr[1] * 0.000000001;
diff /= scales[shift];
return diff;
}
}
void print_number(DWORD n,DWORD res[],int len,double exp)
{
char buff[PRECISION_DIV9*9+9];
char *p;
int i,s;
INT64 e;
p=buff;
p+=sprintf(p,"%d",res[0]);
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[PRECISION_DIV9+1]; //fac(i)的结果的数组
DWORD time;
time=GetTickCount();
min_diff=1;
res_exp=0;
fac[0]=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[j]=fac[j-2];
fac[0]=carry / DEC_RADIX;
fac[1]=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[j]=fac[j-1];
fac[0]=carry;
if (fac_len<PRECISION_DIV9)
fac_len++;
else
res_exp+=9;
}
bc=Uint32BitCount(fac[0]);
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;
}

复制代码

liangbch · 发表于 2012-12-12 21:34:01

刚才仔细看了一下无心人的代码，明白了其所用的算法。其主要算法如下，从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

账号		自动登录	找回密码
密码			欢迎注册

[讨论] 阶乘和圆周率

评分

评分