关于zju1007.奇怪！

[jwd_1_cool]

这道题目我以前已经通过了！
但今天找不到源码，又重新写了下，竟然怎么也AC不了！

我查了下，用mmmcd(超超) 给的那个公式
y(x)=(1-x)((2-x)((3-x)((4-x)((5-x)((6-x)((7-x)((8-x)R+1/40320)+1/35280)+1/4320)+1/600)+1/96)+1/18)+1/4)+1
其中R=∑ 1.0 / (k(k+1)(k+2)(k+3)(k+4)(k+5)(k+6)(k+7)(k+8)(k+x))

但还是答案错。 谁能帮忙给个程，看看！谢谢！特别是 mmmcd(超超) 兄弟的。

[jwd_1_cool]

通过了！但是 uva 上有道类似题目！ uva731 精度要求更高，我没通过！哪位仁兄有空看看！

[BlueSky2008]

帮你召唤一下 mmmcd(超超) ：）

[mmmcd]

那些系数的规律是：1/(n*n!)
4=2*2！
16=3*3！
18=4*4！
...
35280=7*7！
...

[mmmcd]

打错了：
18=3*3！
96=4*4！

[jwd_1_cool]

不知道为什么，我用你的公式，最后一位是错的！无论N取几。

但少用加速几个就可以通过了！

所以我还好是不知道如何解决UVA 那道！你看看！
http://online-judge.uva.es/p/v7/731.html

[jwd_1_cool]

#include<stdio.h>
int main()
{
long k; /* dummy summation variable */
long n; /* number of terms to sum */
double kd; /* double version of k */
double x; /* the variable */
double sum; /* used for computation summations */
double px, p1; /* psi(x) and psi(1) */
double dpx; /* value of the first accelerated series at x */
double dp2,dp3,dp4,dp5,dp6,dp7,dp8; /* value of the first accelerated series at 2 */
double ddpx,dddpx,ddddpx,dddddpx,ddddddpx,dddddddpx,ddddddddpx; /* value of the second accelerates series at x */
long i; /* counter for elements of table */

p1 = 1.0;
dp2 = 0.25;
dp3=1.0/18.0;
dp4=1.0/96.0;
dp5=1.0/600.0;
dp6=1.0/4320.0;
dp7=1.0/35280.0;
dp8=1.0/40320.0;
n = 12000;

/* loop for the 2001 table entries */
for (i=0; i<2001; i++ )
{
x = i * 0.001;
sum = 0.0;
for( k=1; k<n; k++ ){
kd = k;
sum += 1.0/(kd*(kd+1)*(kd+2)*(kd+x) );
}
//ddddddddpx=(8.0-x)*sum;
//dddddddpx=(7.0-x)*(ddddddddpx+dp8);
//dddddddpx=(7.0-x)*sum;
//ddddddpx=(6.0-x)*(dddddddpx+dp7);
//dddddpx = (5.0-x)*(dp6+ddddddpx);
//ddddpx = (4.0-x)*(dp5+dddddpx);
//dddpx = (3.0-x)*(ddddpx+dp4);
//dddpx= (3.0-x)*sum;
//ddpx = (2.0-x)*(dp3+dddpx);
ddpx=(2.0-x)*sum;
dpx = (1.0-x)*(dp2+ddpx);
px = p1+dpx;
printf("%5.3f %16.12f\n", x, px );
}
return 1;
}
不知道为什么再上去，数据就出错！火死了！

[boodweb]

呵呵，我一开始也没想到，后来请教了starfish
这是说明，其实还可以推下去的，可以更快
/*
Algorithm:
We know that phi(1) = 1.
Let f(x) = phi(x) - phi(1),
then we get f(x) = (1-x)* ∑[ 1/(k * (k+1) * (k + x)) ],
and f(1) = 0, f(2) = -0.25.
Set g(x) = f(x) - (x-1)*f(2),
and we can get g(x) = (1-x) * (2-x) * ∑[1/(k * (k + 1) * (k + 2) * (k + x) )
].
So phi(x) = f(x) + phi(1) = g(x) + (x-1)*f(2) + phi(1)

More over, we use the inequation given by the problem,
and we can see that ∑1/(k^r) < 1/ [(r-1)* ((n-1)^(r-1))], for k = n-1, ...,∞
.
Let ∑1/(k^r) < 1/ [(r-1)* ((n-1)^(r-1))] < ESP,
where ESP = 1e-12 given by the problem,
and we can calculate the max k we should count.

*/

[jwd_1_cool]

继续推下去是
f(3)=-(1/18)
g(x) = (1-x) * (2-x) * (3-x)*∑[1/(k * (k + 1) * (k + 2) * (k+3) * (k + x) )


phi(x) = f(x) + phi(1) = g(x) + (x-1)*f(2) + (x-2)*f(3) + phi(1)


是吗???

[jwd_1_cool]

不对,如何再推下去,加速啊!

[mmmcd]

里面说：
As an example, the sample output below shows 4 acceptable lines out of 3001, which might appear in the output file.

The values of x should start at 0.00 and increase by 0.1 until the line with x=300.00 is output.


我也没解决，
期待高手

[jwd_1_cool]

等待！

[boodweb]

应该不是很难的吧，自己多想想吧，我懒得重推了

[mmmcd]

可是ZJU 1007 AC的程序过不了UVA731

[jwd_1_cool]

讲一下推导公式啊，否则我们怎么想啊！
你给出的东西里面 又没讲清楚，我想也没想不出来啊！

[mmmcd]

1/(k(k+x)) - 1/(k(k+1)) = (1-x)/(k(k+1)(k+x))

sigma 1/(k(k+x)) - sigma 1/(k(k+1)) = (1-x) * sigma 1/(k(k+1)(k+x))

sigma 1/(k(k+x)) = (1-x) sigma 1/(k(k+1)(k+x)) + sigma 1/(k(k+1))
= (1-x) sigma 1/(k(k+1)(k+x)) + 1/1 - 1/2 + 1/2 - 1/3 + ...
= (1-x) sigma 1/(k(k+1)(k+x)) + 1/1

1/(k(k+1)(k+x)) - 1/(k(k+1)(k+2)) = (2-x)/(k(k+1)(k+2)(k+x))

sigma 1/(k(k+1)(k+x)) - sigma 1/(k(k+1)(k+2))
= (2-x) * sigma 1/(k(k+1)(k+2)(k+x))

sigma 1/(k(k+1)(k+x))
= (2-x) * sigma 1/(k(k+1)(k+2)(k+x)) + sigma 1/(k(k+1)(k+2))
= (2-x) * sigma 1/(k(k+1)(k+2)(k+x))+1/2/(1*2)-1/2/(2*3)+1/2/(2*3)-1/2/(3*4)+...
= (2-x) * sigma 1/(k(k+1)(k+2)(k+x))+1/2/(1*2)

sigma 1/(k(k+1)(k+2)(k+x)) - sigma 1/(k(k+1)(k+2)(k+3)) =
sigma (3-x)/(k(k+1)(k+2)(k+3)(k+x)) = (3-x)*sigma 1/(k(k+1)(k+2)(k+3)(k+x))

...

sigma 1/(k(k+1)(k+2)...(k+n-1)(k+x)) - sigma 1/(k(k+1)(k+2)...(k+n)) =
sigma (n-x)/(k(k+1)(k+2)...(k+n)(k+x))
=(n-x)*sigma 1/(k(k+1)(k+2)...(k+n)(k+x))

sigma 1/(k(k+1)(k+2)...(k+n-1)(k+x))=
(n-x)*sigma 1/(k(k+1)(k+2)...(k+n)(k+x)) + 1/n/n!

[boodweb]

的确，不管加不加速，x=300.0的时候一看就不对了，可能是积累误差太大了，看来要改进，让我想想
btw: uva真夸张,zju只到2.0阿...

[boodweb]

开始讨厌C++了，明明x=5.0,结果floor(x)=4.0，一点不人性化，害的我还要写自己的floor...还好终于过了

基本思想跟starfish得差不多：计算phi(x)的时候，用上面starfish的方法，先推导
f(x)=phi(x)-phi(m)，其中m=floor(x)
再计算
g(x)=f(x)-(x-m)f(m+1)=(m-x)(m+1-x)∑[1/k/(k+m)/(k+m+1)/(k+x)]
phi(x)=f(x)+phi(m)=g(x)+(x-m)f(m+1)+phi(m)
即可

#include <stdio.h>
#include <iostream>
#include <cmath>

using namespace std;

const double ESP=1e-13;
const double E=1e-3;
const int r=4;
double phi_m=1.0; //phi(m)
double myfloor(double x)
{
double t=floor(x);
if(x-t>E+0.9) return t+1.0;
return t;
}
double f(int m)
{
double ret=0.0;
int k;
for(k=1;k<=m-1;k++) ret+=1.0/k;
ret*=-1.0/(m-1)/m;
ret+=1.0/m/m;
return ret;
}
//phi(x)=Σ(1/k/(k+x)),k=1..infinity
double phi(double x,int n)//r=4
{
double ret=0.0;
int k;
int m;
double t=myfloor(x);
if(x<=2.0+E) m=1;
else if(x-t<=E)
m=t-1;
else
m=t;
for(k=1;k<=n;k++)
{
ret+=1.0/k/(k+m)/(k+m+1)/(k+x);
}
ret*=(m-x)*(m+1-x);
ret+=phi_m-(m-x)*f(m+1);
if(x>=2.0-E && x-t<E)
phi_m=ret;
return ret;
}
int main()
{
double x;
int n;
n=(int)pow(1.0/(r-1.0)/ESP,1.0/(r-1.0))+2;
for(x=0.0;x<=300.0+E;x+=0.1)
{
printf("%6.2f %16.12f\n", x, phi(x,n) );
}
return 0;
}