ACM Numerical Summation of a Series (ZOJ-1007) 的N阶求解算法

把 zju 1007 Numerical Summation of a Series 运行时间缩到 0ms 级别

主要问题是如何减少求和需要的级数项

直接求解记为0-阶算法的话，至少要 2e12 项求和才能达到 0.5e-12 的精度

题目给出了用y(x)-y(1)的算法，记为1-阶算法吧，使用y(1) 导出 y(x).

余项为 o(1/(2*n^2)) 仍然要至少 1e6 项求和。

把 1007 在 ZOJ 上的运行时间刷成 0ms 的思路

完全可以顺着这个思路，使用M-阶公式进行插值。

使用y(1), y(2), y(3), ..., y(M) 来导出 y(x).

推导如下：

/*  1-order Indirect Sum:  R = o[1/(2*n^2)]
 *    (y(x)-y(m))/(m-x) = sum_{k=1, inf}{u(m; k,x)}; u(m; k,x) = 1/(k*(k+m)*(k+x))
 *  2-order Indirect Sum: [ using y(1), y(2); y(m) = (1 + 1/2 +...+ 1/m)/m ]
 *    t(k,x) = 1/(k(k+1)(k+2)(k+x));
 *    v(k,x) = (1-x)(2-x)*[ t(k,x) ] + (2-x) * v(k,1) - (1-x) * v(k,2);
 *    R = (1-x)(2-x)*sum_{k=n+1, inf}{ t(k,x) } < (1-x)(2-x)/(3*n^3)
 *    When 0 <= x <= 2, N > power(1*2/(3*0.5e-12), 1/3).

 *  N-order Indirect Sum:
 *    prod_{m=1,M}{1/(m-x)} - prod_{m=1,M}{1/(k+m)}
 *     = sum_{m=1,M}{ K(m)*[1/(m-x) - 1/(k+m)] }; where K(m) = prod_{i=1,M; i!=m}{1/(i-m)}
 *  ==> 1/(k(k+x)) = prod_{m=1,M}{ 1/(m-x) } * ( 1 / [prod_{m=1,M}{1/(k+m)} * k(k+x)] ) ...
 *     + sum_{m=1,M}{ 1/(k(k+m)) * K(m) * L(m,x)}; where L(m,x) = prod_{i=1,M; i!=m}{(i-x)}
 *    Obviously, K(m) * L(m,m) = 1;
 *    R = (1-x)...(M-x)*sum_{k=n+1, inf}{1/(k(k+1)...(k+M))} < (1-x)...(M-x)/((M+1)*n^(M+1))
 *    When 0 <= x <= 2, N > power(1*2*...*M/((M+1)*0.5e-12), 1/(M+1)).
 */

算法阶数和迭代次数的关系怎样呢？

相加使用的项数用这个公式计算：N = int(1.2 * pow(prod / (ORDER + 1) * 2e12, 1.0 / (ORDER + 1)));

 *  Choose an ORDER with a smaller M*N. [M = 10 is prefered] (slack factor = 1.2)
 *    ORDER(M)  :   1        2          3       4       5      6     7    8    9  (10)  11  12
 *    slacked N : 1.2e6  13,207  1,579  473  221  133  92  71  58  (50)   44  40

选用 10 阶算法看起来已经具有相当小的计算开销，再增加阶数意义也不大。

求前 50 项和就能达到 0.5e-12 的精度要求，每次计算花费的迭代次数是很少的。

ORDER=10 的时候在 ZOJ 上跑出来的结果是 0ms 和 164KB

贴下C代码：

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

#define ORDER 10
double summation(int N, double x, double* ref) {
    double kk, yx, prod, sum = 0;
    int i, m;
    for (kk = N; kk > 0.5; kk--) {			// (!) 'int k=N; k>0; k--' will overflow.
        prod = kk * (kk + x);
        for (i = 1; i <= ORDER; i++) prod *= (kk + i);
        sum += 1 / prod;					// sum of 'prod_{m=1,M}{1/(k+m)} * k(k+x)'
    }
    for (i = 1; i <= ORDER; i++) {
        for (m = 1 + (i == 1); m <= ORDER; m += (i == ++m))
            ref[m - 1] *= ((i - x) / (i - m));		// prod_{i=1,M; i!=m}{(i-x)/(i-m)}
        sum = sum * (i - x);
    }
    for (yx = sum, i = 1; i <= ORDER; i++)
        yx += ref[i - 1];
    return yx;
}
int main(void) {
    double tmp, prod, ym[ORDER], ref[ORDER];		// ref: y(m); m = 1, 2, ..., ORDER
    int i, N;
    tmp = prod = 1;
    for (i = 1; i <= ORDER; tmp += 1.0 / ++i)		// (!) '1/++i' will always be 0
        prod *= i, ym[i - 1] = tmp / i;
    N = int(1.2 * pow(prod / (ORDER + 1) * 2e12, 1.0 / (ORDER + 1)));	// slack by 1.2 for round-off
    for (i = 0; i < 2001; i++) {
        tmp = 0.001 * i;
        memcpy(ref, ym, ORDER * sizeof(double));
        printf("%5.3f %16.12f\n", tmp, summation(N, tmp, ref));
    }
    return 0;
}

原题大致是如下的：

-------------------------------------------------------------------------------------------------------------------------------------

 Numerical Summation of a Series     Time Limit: 10 Seconds   Memory Limit: 32768 KB   Special Judge

-------------------------------------------------------------------------------------------------------------------------------------
Produce a table of the values of the series: y = sum_{k=1,inf}{ 1 / [k*(k+x)]}, for the 2001 values of
x= 0.000, 0.001, 0.002, ..., 2.000. All table entries must be under absolute error of 0.5e-12. This
problem is from Hamming (1962), when mainframes were very slow by today's microcomputer standards.
  Input: This problem has no input.
  Output: Formatted as two columns with x and y(x) printed as in the C printf or the Pascal writeln.
e.g. printf("%5.3f %16.12f\n", x, psix ) or writeln(x:5:3, psix:16:12)
Here are 4 acceptable lines out of 2001. (x starts at 0.000 and increases by 0.001 until 2.000)
0.000   1.644934066848
...
0.500   1.227411277760
...
1.000   1.000000000000
...
2.000   0.750000000000
  Hint
  Roundoff would render any typical double precision computation useless for the desired precision.
To improve the convergence, note that 1/(k*(k+1)) = 1/k - 1/(k+1). The following inequality is helpful
in determining how may items are required in summing the series above:
  sum_{k=n,inf}{1/k^r} < int_{n-1,inf}{1/x^r}dx = {-1/[(r-1)*x^(r-1)]}_{n-1,inf} (k>1; r>=1)