收藏 不显示删除回复显示所有回复显示星级回复显示得分回复 高效递归算法求最接近π的分数[问题点数:100分，结帖人:northwolves]

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

答案:85910/27346

这是怎么做的？

怎么实现无重复出现历遍数字？

猜测这可能用到动态规划或是贪心算法,   先从分子和分母的最高位开始,逐位确定.因为数字所在的位置越高,对其最终分数值离PIE的偏差的影响就越大.   具体的还没考虑过

其实枚举就好了，数据不大。

PIE是常量，数据又不大，遍吧～
2903040

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

n值是什么???

这个就是Pi的连分数逼近得到的值355/113
下一个数103993/33102越界了

贴个递归遍历的

#include <stdio.h>
#include <math.h>

inline void swap(int &e1, int &e2) { int t = e1; e1 = e2; e2 = t; }

void  permutation(int *figs, int start, int &denominator, int &numerator)
{
if(start == 9) {
int d=0, n=0;
for(int i = 9, pown=1; i >0; i--, pown *=10)
{
d += figs[i]*pown; i--;
n += figs[i]*pown;
}
double p1 = fabs(double(d)/n -3.1415926);
double p2 = fabs(double(denominator)/numerator-3.1415926);
if(p1 < p2)
{
denominator = d;
numerator  = n;
}

}
else  {
int n = (start <2)?  9 : 10;
for(int i = start; i < n; i++)
{
swap(figs[start], figs[i]);
permutation(figs, start+1, denominator, numerator);
swap(figs[start], figs[i]);
}
}
}

void main()
{
int figs[10] = {1,2,3,4,5,6,7,8,9,0};
int denominator=12345, numerator=67890;

permutation(figs, 0, denominator, numerator);
printf("%d/%d\n",denominator, numerator);
}

狼的参考答案是错误的，
比如58908/18751就比85910/27346好

穷举结果显示是99733/31746

#include <stdio.h>
#include <math.h>

#define PI 3.1415926535898
int main(){
    int i;
    double min_error=1.0;
    int min_u,min_d;
    for(i=10000;i <=31831;i++){
        double v=floor(i*PI+0.5);
        int k=(int)v;
        if(k>=100000)k=99999;
        double err=fabs((double)k/i-PI);
        if(err <min_error){
            min_error=err;
            min_u=k;
            min_d=i;
        }
    }
    printf("Best Estimation:%d/%d\n",min_u,min_d);
}

To mathe:

狼的参考答案是错误的，
比如58908/18751就比85910/27346好
---------------------
mathe兄,分子分母是0-9的不重复全排列

不好意思，没有仔细看。结果是正确的，挺有意思，正好等于355/113
#include <stdio.h>
#include <math.h>
#include <algorithm>
using namespace std;

char digits[10]={0,1,2,3,4,5,6,7,8,9};
#define PI 3.1415926535898
#define NUM(a,b,c,d,e)  ((((((a)*10)+(b))*10+(c))*10+(d))*10+(e))
#define UP(d)          NUM(d[0],d[1],d[2],d[3],d[4])
#define DOWN(d)        NUM(d[5],d[6],d[7],d[8],d[9])
int main(){
    double min_error=1.0;
    int min_u,min_d;
    do{
        int d,u;
        double x;
        if(digits[0]==0||digits[5]==0)continue;
        if(digits[0] <3*digits[5])continue;
        u=UP(digits);
        d=DOWN(digits);
        x=fabs((double)u/d-PI);
        if(x <min_error){
            min_error=x;
            min_u=u;
            min_d=d;
        }
    }while(next_permutation(digits,digits+10));
    if(min_error <1.0){
        printf("Found %d/%d\n",min_u,min_d);
    }
}

贴个效率较高的递归：

/*
*  回溯+最大-最小剪枝
*/
#include <cstdlib>
#include <iostream>
#include <windows.h>

using namespace std;

const double  pi=3.1415926;
double    best_pi=3.14;
pair <int,int>  best;
int  index[10];
int  pos;

double _abs( double n ){
      return  n>0? n : -n;
};

int  power[10][5]={
            0, 0, 0, 0, 0,
      1,10,100,1000,10000,
      2,20,200,2000,20000,
      3,30,300,3000,30000,
      4,40,400,4000,40000,
      5,50,500,5000,50000,
      6,60,600,6000,60000,
      7,70,700,7000,70000,
      8,80,800,8000,80000,
      9,90,900,9000,90000
      };
     
void  f( int N, int D, int pos ){
      if( pos==-1 ){
          if( _abs(double(N)/D-pi) <_abs(best_pi-pi) ){
              best=pair <int,int>(N,D);           
              best_pi=double(N)/D;
          }
          return ;
      }
      int  i, j, temp1, temp2;
      for( i=0; i <10; ++i )
          if( index[i]==0 ){
              index[i]=1;
              for( j=0; j <10; ++j )
                    if( index[j]==0 ){
                        index[j]=1;
                        temp1=N+power[i][pos];
                        temp2=D+power[j][pos];
                        if( temp1+power[1][pos]>best_pi*temp2
                        && temp1 <best_pi*(temp2+power[1][pos]) )
                            f( temp1, temp2, pos-1 );
                        index[j]=0;
                    }
              index[i]=0;
          }
};
     
int main(int argc, char *argv[]){
    int time=GetTickCount();
    f( 0, 0, 4 );
    cout < <best.first < <"/" < <best.second < <"==" < <best_pi < <endl;
    cout < <"Use time : " < <GetTickCount()-time < <" ms\n";
    system("PAUSE");
    return EXIT_SUCCESS;
}

哥哥们，麻烦你们代码加点注释啊，看不懂啊。。。

我有一个提议，不知是否妥当。请大家看看……
我觉得可以利用这个问题考验一下大家的计算能力。
设这样的一个分数为Pn，它满足：
1、它的分子和分母所包含的所有数字，恰好由n套0到9的数字构成。（这里去掉了首位数字不能为0的条件。所谓“n套0到9的数字”是指n个0、n个1、……和n个9的一共10n个数字。）
2、它是满足上面1的条件中的所有数中最接近于圆周率Pai的那个。
上面大家已经计算出了P1=85910/27346，那么P2呢？其它的呢？大家也可以预测一下计算力的极限在哪里。
呵呵。

这么推广，势必要用大数运算库了。

对2套0到9的数字构成的分子与分母，算了一下，8416355741/2679009238=3.1415926535897969904648757317947。
由于pi我用的是double类型，可能精度不够，因此p2未必就是8416355741/2679009238。

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

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的精度已经挺了不起了。

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

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

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

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

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

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

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

学习一下牛人的算法。

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

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

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

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

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

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

#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);
}

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

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

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

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

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

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

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

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

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