Produce a table of the values of the series
or the 3001 values of x, x = 0.0, 0.1, 0.2, …, 300.00. All entries of the table must have an absolute error less than 10−10
. This problem is based on a problem from Hamming (1962), when mainframes were very slow by today’s microcomputer standards.
Format of function:
void Series_Sum( double sum[] );
where double sum[] is an array of 3001 entries, each contains a value of ϕ(x) for x = 0.0, 0.1, 0.2, …, 300.00.
Sample program of judge:
#include <stdio.h>
void Series_Sum( double sum[] );
int main()
{
int i;
double x, sum[3001];
Series_Sum( sum );
x = 0.0;
for (i=0; i<3001; i++)
printf("%6.2f %16.12f\n", x + (double)i * 0.10, sum[i]);
return 0;
}
/* Your function will be put here */
Sample Output:
0.00 1.644934066848
0.10 1.534607244904
...
1.00 1.000000000000
...
2.00 0.750000000000
...
300.00 0.020942212934
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
which implies ϕ(1)=1.0. One can then produce a series for ϕ(x)−ϕ(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 again on the second series to produce a third sequence, which converges even more rapidly than the second.
The following equation is helpful in determining how may items are required in summing the series above.
解法:
- 注意到当x是整数时, ϕ ( x ) = ∑ i = 1 x 1 i \phi(x)=\sum\limits_{i=1}^x\frac{1}i ϕ(x)=i=1∑xi1
- 需要用公式将x不是整数的情况收敛到整数。
ϕ ( x ) − ϕ ( i ) = ∑ i − x k ( k + x ) ( k + i ) \phi(x)-\phi(i)=\sum\frac{i-x}{k(k+x)(k+i)} ϕ(x)−ϕ(i)=∑k(k+x)(k+i)i−x
∑ 1 k ( k + x ) ( k + i ) = ∑ k + x + i + 1 − x k ( k + x ) ( k + i ) ( k + i + 1 ) = ∑ 1 k ( k + i ) − ∑ i + 1 − x k ( k + x ) ( k + i ) ( k + i + 1 ) = ϕ ( i ) − ( i + 1 − x ) ∑ 1 k ( k + x ) ( k + i ) ( k + i + 1 ) \sum\frac{1}{k(k+x)(k+i)}=\sum\frac{k+x+i+1-x}{k(k+x)(k+i)(k+i+1)}=\sum\frac{1}{k(k+i)}-\sum \frac{i+1-x}{k(k+x)(k+i)(k+i+1)}=\phi(i)-(i+1-x)\sum \frac{1}{k(k+x)(k+i)(k+i+1)} ∑k(k+x)(k+i)1=∑k(k+x)(k+i)(k+i+1)k+x+i+1−x=∑k(k+i)1−∑k(k+x)(k+i)(k+i+1)i+1−x=ϕ(i)−(i+1−x)∑k(k+x)(k+i)(k+i+1)1
∑
1
k
(
k
+
x
)
(
k
+
i
)
(
k
+
i
+
1
)
=
∑
k
+
x
+
i
+
2
−
x
k
(
k
+
x
)
(
k
+
i
)
(
k
+
i
+
1
)
(
k
+
i
+
2
)
=
∑
i
+
2
−
x
k
(
k
+
x
)
(
k
+
i
)
(
k
+
i
+
1
)
(
k
+
i
+
2
)
+
∑
1
k
(
k
+
i
)
(
k
+
i
+
1
)
(
k
+
i
+
2
)
\sum \frac{1}{k(k+x)(k+i)(k+i+1)}=\sum\frac{k+x+i+2-x}{k(k+x)(k+i)(k+i+1)(k+i+2)}=\sum\frac{i+2-x}{k(k+x)(k+i)(k+i+1)(k+i+2)}+\sum\frac{1}{k(k+i)(k+i+1)(k+i+2)}
∑k(k+x)(k+i)(k+i+1)1=∑k(k+x)(k+i)(k+i+1)(k+i+2)k+x+i+2−x=∑k(k+x)(k+i)(k+i+1)(k+i+2)i+2−x+∑k(k+i)(k+i+1)(k+i+2)1
其中
∑ 1 k ( k + i ) ( k + i + 1 ) ( k + i + 2 ) = 1 2 ( ∑ 1 k ( k + i ) ( k + i + 1 ) − ∑ 1 k ( k + i + 1 ) ( k + i + 2 ) ) = 1 2 ( ϕ ( i ) − ϕ ( i + 1 ) − ( ϕ ( i + 1 ) − ϕ ( i + 2 ) ) ) \sum\frac{1}{k(k+i)(k+i+1)(k+i+2)}=\frac12(\sum\frac{1}{k(k+i)(k+i+1)}-\sum\frac{1}{k(k+i+1)(k+i+2)})=\frac12(\phi(i)-\phi(i+1)-(\phi(i+1)-\phi(i+2))) ∑k(k+i)(k+i+1)(k+i+2)1=21(∑k(k+i)(k+i+1)1−∑k(k+i+1)(k+i+2)1)=21(ϕ(i)−ϕ(i+1)−(ϕ(i+1)−ϕ(i+2)))
可得收敛式:
ϕ
(
x
)
=
ϕ
(
i
)
+
(
i
−
x
)
(
ϕ
(
i
)
−
(
i
+
1
−
x
)
(
∑
i
+
2
−
x
k
(
k
+
x
)
(
k
+
i
)
(
k
+
i
+
1
)
(
k
+
i
+
2
)
+
1
2
ϕ
(
i
)
−
ϕ
(
i
+
1
)
+
1
2
ϕ
(
i
+
2
)
)
)
)
\phi(x)=\phi(i)+(i-x)(\phi(i)-(i+1-x)(\sum\frac{i+2-x}{k(k+x)(k+i)(k+i+1)(k+i+2)}+\frac12\phi(i)-\phi(i+1)+\frac12\phi(i+2))))
ϕ(x)=ϕ(i)+(i−x)(ϕ(i)−(i+1−x)(∑k(k+x)(k+i)(k+i+1)(k+i+2)i+2−x+21ϕ(i)−ϕ(i+1)+21ϕ(i+2))))
提供一种具体算法的思路:
- 外循环计算 ϕ ( i ) \phi(i) ϕ(i)
- 内循环计算 x x x在 i i i到 i + 1 i+1 i+1时的 ϕ ( x ) \phi(x) ϕ(x)
- 计算1000项精度即可满足条件
- 积分式是用来估算项数的,不是用来算n到无穷的值的,n以后的值可忽略。
上代码
double phi(int x) {
double ret = 0.0;
for (int i = 1; i <= x; i++) {
ret += 1.0 / i;
}
return ret / (double)x;
}
void Series_Sum(double sum[]) {
int n = 1000;
sum[0] = 1.644934066848;
for (int i = 0; i < 300; i++) {
double phii = phi(i + 1);
double phii1 = phi(i + 2);
double phii2 = phi(i + 3);
for (int j = 1; j <= 10; j++) {
int x = 10*i + j;
double dx = x * 0.1;
sum[x] = 0;// ((double)i + 3.0 - dx) / (4.0 * (double)n * (double)n * (double)n * (double)n);
for (int k = 1; k < n; k++) {
sum[x] += 1.0 / ((double)k * ((double)k + dx) * ((double)k + (double)i + 1.0) * ((double)k + (double)i + 2.0) * ((double)k + (double)i + 3.0));
}
sum[x] = phii + ((double)i + 1.0 - dx) * (((double)i + 2.0 - dx) * (((double)i + 3.0 - dx) * sum[x] + 0.5 * phii - phii1 + 0.5 * phii2) + phii - phii1);
}
}
}