找了一下Pari/Gp里面有intformal函数可以对多项式积分，不过里面要求所有符号有个默认顺序，而积分必须对默认符号积分，其实起来不是很方便，所以下面的代码使用前必须保证gw,gc变量没有定义，或者gw定义在gc之前，也就是调用函数reorder()可以看到gc没有定义或者gw,gc都定义了但是gw在前面
symbolp(n)=
{
local(p,q,p1,p2,p3);
global(gw,gc,gt1,gt2);
p=1;
for(u=1,n,
q=intformal(p,gw);
q=subst(q,gc,gt1);
q=subst(q,gw,gt2);
p1=subst(q,gt1,gc);
p1=subst(p1,gt2,(gc+gw)/2.0);
p2=subst(q,gt1,1.0-gw);
p3=subst(p2,gt2,1.0-gw);
p2=subst(p2,gt2,1-(gc+gw)/2.0);
p3=p3-p2;
p2=subst(q,gt1,gw);
p2=subst(p2,gt2,gw);
p=p1+p2+p3;
);
p
}
getdist(n)=
{
local(p,d1,d2,d3,d4);
global(gw,gc);
p=symbolp(n);
d1=p*(-gc^2/2.0+gw^2+gc-gw);
d2=p*(-gc^2+2*gw*gc-gw^2+gc-gw);
d3=intformal(d1);
d4=intformal(d2);
d1=subst(d3,gw,gc/2.0);
d2=subst(d4,gw,gc);
d4=subst(d4,gw,gc/2.0);
d2=d2-d4;
d1=d1+d2;
d1=intformal(d1);
2.0*subst(d1,gc,1.0)
}
然后计算结果如下:
(16:06) gp > getdist(2)
%17 = 0.2114583333333333333333333333
(16:07) gp > getdist(10)
%18 = 0.2203686807254583318975768873
(16:07) gp > getdist(20)
%19 = 0.2204074196271448826430741941
(16:07) gp > getdist(30)
%20 = 0.2204077318388700882261725891
(16:07) gp > getdist(40)
%21 = 0.2204077356077135199957228968
(16:07) gp > getdist(50)
%22 = 0.2204077356585942289207798388
(16:07) gp > getdist(60)
%23 = 0.2204077356593274980124179385
(16:07) gp > getdist(70)
%25 = 0.2204077356593385397246878697
(16:33) gp > getdist(99)
%1 = 0.2204077356593387142224517526
(16:10) gp > getdist(100)
%26 = 0.2204077356593387141426361209
也就是我们可以得出g(a,b)=(a+b)/2时的高精度结果，也就是每秒4.537次,同Fans的模拟结果匹配
14#里面zgg定义的gp(x1,x2)=(x1+x2+q/2)/(q+2)好像有问题，比如在极端情况x1~=x2时，gp(x1,x2)应该也取x1或x2,但是他这里定义的函数不满足，所以结果应该不对。不知道是不是他定义的gp函数不是我定义的g函数？
而对于Fans前面定义的情况，由于对应的g函数太复杂，符号计算太难。最好我们改成线性的，至少弄成折线，比如(a+b)<1时，g(a+b)=c*a+(1-c)*b,不然g(a+b)=(1-c)*a+c*b,其中c是一个大于1/2的数，这样定义的g就应该偏向于选择离1/2远的点。反之如果c小于1/2，那么偏向于选择离1/2近的点。
当然，由于g表达式更加复杂了，对应的符号计算代码也要修改了
发现只要g是分段，那么符号计算的结果就会越来越复杂，不是很合适
由于对于一般的g,直接积分计算很困难，我们可以考虑对p进行多项式拟合。也就是假设p是关于w,c的n次多项式(O(n^2)个系数），然后对p采样（采样数目大于待定系数），代入等式得到线性方程组，当然方程数目大于变量数目，Ax=b,然后解方程A'Ax=A'b应该可以得出一个比较好的p(w,c)
上面多项式可能不如改成三角函数，不然计算可能会比较变态。唯一比较麻烦的是在g比较复杂时，最后一步平均移动距离的计算公式可能比较复杂，实在不行，这一步可以采用数值积分。
而34#中得出g(a,b)=(a+b)/2时最优应该是错误的方法。理由只能是p和g是相互约束的，给定g即唯一确定p,同样给定p就唯一确定g,所以我们不能使用假设p是给定的条件求出g的最优值。
而34#方法得出错误的结论说明32#中迭代的方法也不能使用，因为它也必然会收敛到34#中的结果。
又弄了一个10^8的豆的，感觉上需要验证随机数的可靠性和考虑double的误差积累了。所以请Fans也有空算一下权重p的极小值吧。呵呵。
源码如下，我用vc2005：#include "stdafx.h"
#include
#include
#include
t ...
zgg___ 发表于 2013-2-19 15:45 http://bbs.emath.ac.cn/images/common/back.gif
19#中经常会出现越过一个豆子吃另外一个豆子，去除这种情况，可以得出一秒可以超过5.9次
这些天抽空对25层和30层的方法又检查了一下，觉得有的结论是对的，呵呵。
先说下我在30层为什么认为存在问题？
主要有二：Q1、25层中对平均吃豆时长的计算结果和试验结果相差大，比如说，对p=0时，即w总去吃最近的豆的情况，此时按照30层公式8的结果得0.2185013458，而试验结果是0.220+（和试验结果差距这么大，一定是有问题的呢，这个不光是Fans的试验，我自己的程序也是这个结果，呵呵。）Q2、是用25层公式7算出的曲线和30层代码第3行的结果是不一致的，不一致性在30层最后输出的图里可以看到，因为都是M的计算，所以和试验都搭不上边，也一定有问题的呢，呵呵。综合这两点，我也对“模型”产生过怀疑，呵呵，但是，我并不觉得“模型”有什么问题，可能在计算中会有计算错误，但是“模型”还是比较要命的，呵呵。
对于Q2的问题，发现了问题所在（正如mathe在42层说的），是因为在计算中没有考虑到g有可能落在x1和x2之外，比如说，当出现1/2<g<w<x1<x2的时候，w按照g判决的规则应该是去吃x2，而实际上w去吃的是x1，因此就不一致了。而这个反过来说明g判决不完全合理。而当g判决中1/2的权重为0，即p=0时，这个问题是不存在的。
所以我们可以着重计算几个较为简单的问题：
T1、当w吃豆的策略为总去吃最近的豆的时候，平均吃豆时长是多少？（这样就是一个数了，呵呵。）
T2、当w按照T1中策略时，当w吃豆稳定的时候，w吃豆后的分布fw是什么？w没有去吃的豆的分布f是什么？可以具体的去求当w吃豆后，位置位于0到1/4的几率是多少？当w吃豆后，剩下的豆，位置位于0到1/4的几率是多少？（这样就又变成求具体数了，呵呵。）
T3、假设w吃豆的策略是尽可能的去吃位置靠近0的豆子（这是一种笨策略，但是，这个解题上比较简单，或许有一点点启发性，呵呵），那么平均吃豆时长是多少？
对于T1，我的结果是25层中的公式8，我找到了一个更简单的表述，如下图。
对于T2，我的结果是：w吃的豆为1/4，w剩的豆为1/2-17^(-2/3)=0.348748。通式见图。
对于T3，相当于只有一个豆的情况，我的结果是1/3。
最后，我对于Q1，我也倾向于认为对不上的原因是迭代的次数不够，因为在T2和T3的求解数据试验中，迭代次数都达不到，而且我不认为抛弃掉前部数据可以解决。
总之，我觉得至少在p=0时计算或许还是可靠的，而g的判决式的选择存在着一些盲目性，或许g的选择应该满足：1、是x1、x2的函数；2、总是在x1、x2之间。
最后，按照mathe提出的说法，w总去吃最近的豆，这就是最优策略。我不能确定这种说法真的对，我有时候也模糊的有一点这个感觉，但是，数据试验的证据实在太强大了，我不敢相信这种感觉，呵呵。
听说$19#$的策略会越过一个豆子吃另外一个豆子，立马打个补丁重新运行。
目前已运行$48$小时，结果为：
$p=0.478\pm0.005$
$t=0.2191031\pm0.0000005$
其中$p$是$1/2$的权重，$t$是吃一个豆的平均时长。
49# KeyTo9_Fans
这是不是 意味着zgg的就是终极答案了？