高效递归算法求最接近π的分数

[northwolves]

求一个分数，分子5位数（第1位不是0），分母也是5位数（第1位不是0），分子和分母这10个数正好由0到9这10个数字组成（不缺也不重复）。求最接近π值的那个分数。

答案:85910/27346

[medie2005]

已经找到一组更好了：6020794258/1916478335=3.1415926535897939070623514249119

[mathe]

P2我算出来结果是
6020794258/1916478335
我是通过gmp来运算的，不过速度非常慢，已经花了1分半钟了。

[mathe]

medie2005用什么方法算的，看来速度也不快:)

[medie2005]

我还是用回溯+剪枝的方法做的。

其实，原来pi的标准值我是用double类型的，已经搜索得到6020794258/1916478335，用时是8秒，优化后是5秒。
但是，考虑到double类型的精度问题，不敢肯定p2=6020794258/1916478335。
今天上午，照着gmp手册，将类型统一为gmp的浮点类型，再算了一遍，运行结果是6020794258/1916478335。由于只是验算，时间我就没记下来。

[mathe]

还有一个问题，对于n比较大的时候，Pi的精确值也就成为一个问题了。其实还不如将目标数设成另外一个数字，比如
3.123456789012345678900123456789000123456789000012...

[zgg___]

mathe的主意不错，但是我想目标还是先定为Pi吧，毕竟人家的名气大嘛。
下面是Pi的1000位，我想计算到P5也该够用了，供大家拷贝粘贴。
3.1415926535897932384626433832795028841
971693993751058209749445923078164062862
089986280348253421170679821480865132823
066470938446095505822317253594081284811
174502841027019385211055596446229489549
303819644288109756659334461284756482337
867831652712019091456485669234603486104
543266482133936072602491412737245870066
063155881748815209209628292540917153643
678925903600113305305488204665213841469
519415116094330572703657595919530921861
173819326117931051185480744623799627495
673518857527248912279381830119491298336
733624406566430860213949463952247371907
021798609437027705392171762931767523846
748184676694051320005681271452635608277
857713427577896091736371787214684409012
249534301465495853710507922796892589235
420199561121290219608640344181598136297
747713099605187072113499999983729780499
510597317328160963185950244594553469083
026425223082533446850352619311881710100
031378387528865875332083814206171776691
473035982534904287554687311595628638823
537875937519577818577805321712268066130
01927876611195909216420199
可以看到medie2005的P2的精度已经挺了不起了。

[zgg___]

我算P2也是6020794258/1916478335，刚刚算出的，用了半天的时间，呵呵。

[mathe]

我估算了一下，我的程序计算P3，大概需要1个月。（算法应该筒medie没有本质区别）
从P1到P2到P3,每进一步都要花费前面数万倍的时间。
看来如果算法没有本质的突破（好像很难），那么如果通过互联网调动数万机器，也许能够算出P4,但是P5基本上不大可能算出了。

[oo]

学习一下牛人的算法。

不过好象大家没看清lz的限制条件

[medie2005]

如mathe所言，用硬搜的办法是不行了。

是否可以先不管特殊的数字构成这个条件，由连分数求得最优渐近分数，然后再利用数学知识，逐步求得pn？

目前正在考虑用farey序列来做。

大家有什么好想法，不妨说说啊。

[mathe]

我写了个专门针对P3的程序，
使用定点计算，估算了一下，必用GMP的快一倍左右，但是用一台机器还是花费时间太多了。
所以我写了个程序，将计算拆分成22部分。
这样就可以大家一起合作，至少现将P3解决了。
比如编译后的程序时xxx
那么运行的方法是
xxx n
其中n是10~31的一个数字。

[mathe]

#include <stdio.h>
#include <time.h>
#include <math.h>
#include <assert.h>
typedef struct tag_int6{
unsigned x[6];
}int6;

#define NUMS 3
long long d, u;
long long max_u;
int delta;
time_t start_time;
char strd[5*NUMS+1];
char stru[5*NUMS+2];
double err=1.0;
double twotimes64;
long long best_d, best_u;
int6 PI;
int6 N_PI[5*NUMS+1];
int6 NA1_PI[5*NUMS+1];

//dst=10*src
void time10(int6 *dst, int6 *src)
{
long long x0=src->x[0];
long long x1=src->x[1];
long long x2=src->x[2];
long long x3=src->x[3];
long long x4=src->x[4];
long long x5=src->x[5];
long long r;
long long c;
r=(x0<<1)+(x0<<3);
c=r>>32;
dst->x[0]=(unsigned)r;
r=(x1<<1)+(x1<<3)+c;
c=r>>32;
dst->x[1]=(unsigned)r;
r=(x2<<1)+(x2<<3)+c;
c=r>>32;
dst->x[2]=(unsigned)r;
r=(x3<<1)+(x3<<3)+c;
c=r>>32;
dst->x[3]=(unsigned)r;
r=(x4<<1)+(x4<<3)+c;
c=r>>32;
dst->x[4]=(unsigned)r;
r=(x5<<1)+(x5<<3)+c;
c=r>>32;
dst->x[5]=(unsigned)r;
}

void add_by_pi(int6 *dst)
{
unsigned x0=dst->x[0];
unsigned x1=dst->x[1];
unsigned x2=dst->x[2];
unsigned x3=dst->x[3];
unsigned x4=dst->x[4];
unsigned x5=dst->x[5];
long long r;
int c;
r=(long long)x0+PI.x[0];
dst->x[0]=(unsigned)r;
c=(r&0x100000000LL)!=0;
r=(long long)x1+(long long)PI.x[1]+c;
dst->x[1]=(unsigned)r;
c=(r&0x100000000LL)!=0;
r=(long long)x2+(long long)PI.x[2]+c;
dst->x[2]=(unsigned)r;
c=(r&0x100000000LL)!=0;
r=(long long)x3+(long long)PI.x[3]+c;
dst->x[3]=(unsigned)r;
c=(r&0x100000000LL)!=0;
r=(long long)x4+(long long)PI.x[4]+c;
dst->x[4]=(unsigned)r;
c=(r&0x100000000LL)!=0;
r=(long long)x5+(long long)PI.x[5]+c;
dst->x[5]=(unsigned)r;
}

void init_PI(){//Init PI = 3.14159....*(2**128);
PI.x[0]=0x03707344;
PI.x[1]=0x13198a2e;
PI.x[2]=0x85a308d3;
PI.x[3]=0x243f6a88;
PI.x[4]=0x3;
PI.x[5]=0;
}

int test(int L){
long long lu,uu;
int count[10];
int i;
lu = *(long long *)&N_PI[L-1].x[4];//Get the lower bound of u
uu = *(long long *)&NA1_PI[L-1].x[4];//Get the upper bound of u
if(uu>=max_u)uu=max_u-1;
delta=0;
while(lu<=uu){
sprintf(stru,"%*lld",L,lu);
memset(count,0,sizeof(count));
for(i=0;i<L;i++){
count[stru[i]-'0']++;
count[strd[i]-'0']++;
}
for(i=0;i<10&&count[i]<=NUMS;i++);
if(i>=10)return 1;
++lu;
++delta;
}
return 0;
}

void calc(){
double the_err;
unsigned x;
unsigned long long tail;
tail = *(unsigned long long *)&N_PI[5*NUMS-1].x[2];
// if(x==0&&delta==0||
// x==0xFFFFFFFF&&delta==1){
d/=10;
the_err=tail/twotimes64;
the_err-=delta;
the_err/=d;
the_err=fabs(the_err);
if(the_err<err){
err=the_err;
sscanf(stru,"%lld",&best_u);
best_d=d;
fprintf(stderr,"%lld/%lld,err=%g\n",best_u,best_d,err);
}
d*=10;
// }
}

void enum_num(int L){
int i,start,end;
if(L>=5*NUMS){
calc();
return;
}
if(L==3){
time_t now_time = time(NULL);
fprintf(stderr,"%c%c%c:\t%d seconds\n",strd[0],strd[1],strd[2], now_time-start_time);
}
start=(L==0);
end = (L==0)?3:9;
if(L>0){
time10(&N_PI[L],&N_PI[L-1]);//Initialize to d*Pi
}else{
memcpy(&N_PI[0],&PI,sizeof(N_PI[0]));//Initialize = 1*Pi since start from 1.
d=1;
max_u=10;
}
memcpy(&NA1_PI[L],&N_PI[L],sizeof(N_PI[L]));
add_by_pi(&NA1_PI[L]);//Initialize to (d+1)*Pi
for(i=start;i<=end;i++){
strd[L]=i+'0';
strd[L+1]='\0';
if(test(L+1)){
d*=10LL;
max_u*=10LL;
enum_num(L+1);
d/=10LL;
max_u/=10LL;
}
d++;
memcpy(&N_PI[L],&NA1_PI[L],sizeof(N_PI[L]));
add_by_pi(&NA1_PI[L]);
}
d-=(end-start+1);
}

int main(int argc, char *argv[])
{
int i;
start_time = time(NULL);
init_PI();
twotimes64=1.0;
for(i=0;i<16;i++)twotimes64*=16;
if(argc<=1){
enum_num(0);
}else{
int v=atoi(argv[1]);
if(v<10||v>31){
fprintf(stderr,"The parameter should be between 10 and 31\n");
return -1;
}
sprintf(strd,"%02d",v);
d=v*10;
max_u=1000;
{///Preparing N_PI[1];that's v*PI
memcpy(&N_PI[1],&PI,sizeof(PI));
for(i=2;i<=v/10;i++)
add_by_pi(&N_PI[1]);
time10(&N_PI[1],&N_PI[1]);
for(i=0;i<v%10;i++)
add_by_pi(&N_PI[1]);///Insert more call to add_by_pi when strd[0] is larger than 2
}
enum_num(2);
}
printf("Result:%d/%d\n",best_u,best_d);
}

[mathe]

有兴趣的朋友可以选择一两个数据运行一下，
这样当所有部分都运行完了，我们将所有输出中最后的结果收集起来就可以了。
现在我一个程序从头开始（也就是10参数已经在运行中），另外一个程序用了参数31。
要运行的朋友请随便选一个参数，然后把参数贴在这里，避免大家重复相同参数。

[mathe]

不错，VC2005可以识别long long,那么我不用修改就可以直接在VC2005中编译了。
我再占用一下15和25

[mathe]

还有一点需要说明，这个计算方法得到结果要求最终结果同PI的差值不小于10^(-24),不然的化，就说明计算精度不够高，得到结果不一定可靠。

[zgg___]

to medie2005:
“目前正在考虑用farey序列来做。 ”
为啥想起用farey序列来做？是不是上个月IBM题的后遗症，呵呵。

[medie2005]

“上个月IBM题”？呵呵，我没看过。

考虑用farey序列来做，是因为farey序列与Stern-Brocot树有关，找最接近的分数比较方便，不用穷举验证是否是最优。

[mathe]

花了个周末，现在对于P3搜索出来的最优结果是（还没有搜索完毕）：
965148832009552/307216414867379
对于前面提到的22个部分，现在已经搜索了10-17,25-27,30,31
18,19,28,29还在运行中，此外还缺20-24.
运行速度比原先想象的还要慢一些，发现后面部分数据比前面的要运行慢一些。

[mathe]

P3找到一个更好的结果了：
709507129685401/225843133696748

[mathe]

P3我确定是709507129685401/225843133696748了，搜索完毕了。