ACM_1007_Numerical Summation of a Series

Numerical Summation of a Series

---

Time limit: 10 Seconds   Memory limit: 32768K   Special Judge
Total Submit: 37   Accepted Submit: 25  

---

Produce a table of the values of the series

Equation 1

for the 2001 values of x, x= 0.000, 0.001, 0.002, ..., 2.000. All entries of the table must have an absolute error less than 0.5e-12 (12 digits of precision). This problem is based on a problem from Hamming (1962), when mainframes were very slow by today's microcomputer standards.

Input

This problem has no input.

Output

The output is to be formatted as two columns with the values of x and y(x) printed as in the C printf or the Pascal writeln.

printf("%5.3f %16.12f/n", x, psix )		writeln(x:5:3, psix:16:12)

As an example, here are 4 acceptable lines out of 2001.

0.000   1.644934066848
...
0.500   1.227411277760
...
1.000   1.000000000000
...
2.000   0.750000000000

The values of x should start at 0.000 and increase by 0.001 until the line with x=2.000 is output.

Hint

The problem with summing the sequence in equation 1 is that too many terms may be required to complete the summation in the given time. Additionally, if enough terms were to be summed, roundoff would render any typical double precision computation useless for the desired precision.

To improve the convergence of the summation process note that

Equation 2

which implies y(1)=1.0. One can then produce a series for y(x) - y(1) which converges faster than the original series. This series not only converges much faster, it also reduces roundoff loss.

This process of finding a faster converging series may be repeated to produce sequences which converge more and more rapidly than the previous ones.

The following inequality is helpful in determining how may items are required in summing the series above.

Equation 3

这道题目如果不考虑效率方面的问题的话，应该是到很简单的题目，计算sum(1/(k*(k+x)))只要一个存换就可以了，加上对x的循环，只要2重循环。

但是题目要求的精度是小数后12位,运行时间是10秒。

根据Equation 3 的公式，如果直接计算的话，需要到n>10^13才可以达到这个精度，这个是不现实的。

所以Equation 2给出了一个方法，就是计算f(x)-f(1)这样可以推算出公式为

g(x)=f(x)-f(1)=sum((1-x)/(k*(k+1)*(k+x)))

这样程序只需要n>1.4*10^6就可以了。

根据提示完成程序后，一测试，发现需要30多秒。

只好想别的办法。

由于dis(x)=1/(k*(k+1)*(k+x)))-1/(k*(k+1)*k)非常接近，我们可以计算出当k从10^4到无穷大，对dis(x)按k求和的结构远小于10^(-15),所以10000以后的尾数不是计算的关键，我们可以用1/(k*(k+1)*k)来替代。重新编写后，提交，运行只用了0.7s，真正验证了算法的改进是数量级的这个概念。

下面是最后提交的源代码：

                  

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

int main()
{
 double dX = 0.000;
 double dSum = 0.0;

 double dStror1 = 0.0;
 int i, j;

 for(i=1400000; i>10000; i--)//这里需要倒序，主要是double精度的问题。
 {
  dStror1 += 1/((double)(i)*(double)(i)*(double)(i+1));//预先计算10000以后的和
  }
 
 for(i=0; i<2001; i++)
 {
  dSum = dStror1;
  
  for(j=10000; j>=1; j--)
   {
   dSum +=  1/(((double)(j)*(double)(j+1))*double(j+dX));//计算10000以内的
  }

  dSum *= (1-dX);
  printf("%5.3f %16.12f/n", dX, 1+dSum);

  dX += 0.001;
 }

 return 0;
}