- 注册时间
- 2007-12-27
- 最后登录
- 1970-1-1
- 威望
- 星
- 金币
- 枚
- 贡献
- 分
- 经验
- 点
- 鲜花
- 朵
- 魅力
- 点
- 上传
- 次
- 下载
- 次
- 积分
- 40348
- 在线时间
- 小时
|
楼主 |
发表于 2009-6-10 17:40:52
|
显示全部楼层
现在将这个模型的期望数目重新推导了一下,得出一个比较准确的公式,然后用计算机计算
的确在N=32时无法超过1000了(只有824),而N=33时期望数目可以达到1262.
所以这个方法33已经是最好的结果了.-
- #define N 32
- #define MAXR 800
- double maxr;
- double f(double r[4],double x,int n)
- {
- double rt=0.0;
- int i;
- for(i=0;i<n;i++){
- rt*=x;
- rt+=r[n-1-i];
- }
- return rt;
- }
-
- double round(double x)
- {
- return (double)(int)(x+0.5);
- }
-
- #define ERROR 0.0001
- void solve3(double r[4],int m[])
- {
- double dr[3];
- double r0=100.0;
- double err;
- int c;
- int gc=0;
- dr[0]=r[1];dr[1]=2*r[2];dr[2]=3*r[3];
- start:
- c=0;
- do{
- double t=f(r,r0,4);
- err=t/f(dr,r0,3);
- r0-=err;
- c++;
- }while(fabs(err)>ERROR&&c<100);
- if(c>=100){
- r0=(double)rand();
- gc++;
- if(gc>=100){
- m[0]=m[1]=m[2]=-1;
- return;
- }
- goto start;
- }
- m[0]=(int)round(r0);
- dr[2]=r[3];
- dr[1]=r[2]-r0*dr[2];
- dr[0]=r[1]-r0*dr[1];
- double dd=dr[1]*dr[1]-4*dr[0]*dr[2];
- if(dd<0.0){
- m[1]=m[2]=-1;
- }else{
- dd=sqrt(dd);
- m[1]=(int)round((-dr[1]+dd)/(2*dr[2]));
- m[2]=(int)round((-dr[1]-dd)/(2*dr[2]));
- }
- }
-
- void get_m(int k,double q, int m[3])
- {
- int L=N-3*k;
- double r[4];
- double K=k*pow(3.0,(double)(k-2));
- double H=1-2*q*q+2*q*q*q;
- double Q=1-2*q*q+2*q*q*q*q;
- double HL=pow(H,L);
- double QL=pow(Q,L);
- double S=pow(6.0,k-2);
- double R=k*6*S+k*(k-1)/2.0*S*9.0;
- r[0]=-K;
- r[0]+=HL*K;
- r[1]=-3*HL*K;
- r[2]=3.0/2.0*HL*K;
- r[0]-=3.0/2.0*QL*K;
- r[1]+=11.0/2.0*QL*K;
- r[2]-=9.0/2.0*QL*K;
- r[3]=QL*K;
- r[1]+=0.5*QL*K*(K-1);
- r[2]-=3.0/2.0*QL*K*(K-1);
- r[3]+=QL*K*(K-1);
- r[3]+=QL*R*4.0;
- solve3(r,m);
- }
-
- void test(int k,double q, int m)
- {
- int L=N-3*k;
- double K=k*pow(3.0,(double)(k-2));
- double H=1-2*q*q+2*q*q*q;
- double Q=1-2*q*q+2*q*q*q*q;
- double HL=pow(H,L);
- double QL=pow(Q,L);
- double dm=m;
- double S=pow(6.0,k-2);
- double R=k*6*S+k*(k-1)/2.0*S*9.0;
- double r=K*dm-K*dm*(dm-1)*(dm-2)/2.0*HL-dm*(dm-1)*(dm-2)*(dm-3)/4.0*K*QL
- -K*(K-1)*dm*dm*(dm-1)*(dm-1)/4.0*QL-R*dm*dm*dm*dm*QL;
- if(r>=MAXR){
- printf("k=%d,q=%f,m=%d,r=%f\n",k,q,m,r);
- if(r>maxr)maxr=r;
- }
- }
-
- int _tmain(int argc, _TCHAR* argv[])
- {
- int k,s;
- double q;
- for(k=3;k<=8;k++){
- for(q=0.6667;q<0.707;q+=0.0001){
- int m[3];
- get_m(k,q,m);
- for(s=0;s<3;s++){
- if(m[s]>0&&m[s]<100000){
- test(k,q,m[s]);
- }
- }
- }
- }
- printf("maxr=%f\n",maxr);
- return 0;
- }
复制代码 |
|