NA(Numerical Analysis) Lab总结

数值分析的8个LAB做完了，这八道题，可以说是自己做过的课程里边最难的作业了，最难但并不玄学！值得记下来，来纪念自己在这八道题上花费的时间和精力！

这八道题似乎很多年没有变过了，所以，说不定对之后的学弟学妹估计有些帮助吧~

Q1. Numerical Summation of a Series

题目

6-1 Numerical Summation of a Series （30 分）

Produce a table of the values of the series

$\phi(x) = \sum_{k = 1}^{\infty} \frac{1}{k(k+x)} ......(1)$

for 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:

1	void Series_Sum( double sum[] );

where double sum[] is an array of 3001 entries, each contains a value of $\phi(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:

$\frac{1}{k(k+1)} = \frac{1}{k} - \frac{1}{k + 1} ......(2)$

which implies $\phi(1)=1.0$. One can then produce a series for $\phi(x) - \phi(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.

$\sum_{k = n}^{\infty}\frac{1}{k^r} < \int_{n - 1}^{\infty} \frac{1}{x^r}dx (k > 1, r > 1) ......(3)$

思路

先求1，2，…，300的值，整数的值公式比较直观，可以直接求得

按照提示上的方法，推导出递推公式，即用另外几个数的值来表示当前位置的值，比较复杂，也比较有趣

对于在[i, i + 1]区间内的小数，用递推公式运用后几个整数的值求得，我用的是i+1, i+2, i+3, i+4四个数

注意

边界条件，如果用当前数的后四个数，到了297的时候，就不能继续向后走了，我选择的是297，298，299，300四个数来计算这个区间的数
注意C语言中乘法的int型和double型，尽管很注意，但是自己还是犯了错误

点评

应该是最难也是最有意思的一道题了，这样的题让人做的很有感觉！

推导公式还是挺痛苦的，思路理不清，大概用了三个时间段，总共五六个小时吧，所以一开始做的时候不要怀疑自己而去网上找答案，自己推，推出来之后还是很有成就感的！

Q2. Root of a Polynomial

题目

6-2 Root of a Polynomial （40 分）

A polynomial of degree n has the common form as p(x)=c, $p(x) = cnx^n + c{n - 1}x^{n - 1} + \dots + c_1x + c_0$ . Your task is to write a function to find a root of a given polynomial in a given interval.

Format of function:

1	double Polynomial_Root(int n, double c[], double a, double b, double EPS);

where int n is the degree of the polynomial; double c[] is an array of n+1 coefficients $cn, c{n - 1}, …, c_1$ and $c_0$ of the given polynomial; double a and b are the two end-points of the given interval; and double EPS is the accuracy of the root.

The function must return the root.

Note: It is guaranteed that a unique real number r exists in the given interval such that $p(r)=0$.

Sample program of judge:

#include <stdio.h>
#include <math.h>
#define ZERO 1e-13 /* X is considered to be 0 if |X|<ZERO */
#define MAXN 11    /* Max Polynomial Degree + 1 */
double Polynomial_Root(int n, double c[], double a, double b, double EPS);
int main()
{
    int n;
    double c[MAXN], a, b;
    double EPS = 0.00005;
    int i;
    scanf("%d", &n);
    for (i=n; i>=0; i--) 
        scanf("%lf", &c[i]);
    scanf("%lf %lf", &a, &b);
    printf("%.4f\n", Polynomial_Root(n, c, a, b, EPS));
    return 0;
}
/* Your function will be put here */

Sample Input:

1 2	2 1 1.5 -1 0 2

Sample Output:

0.5000

思路

算法用Modified-Newton是可以的

需要求1阶和2阶导数，不过多项式求导很简单

初始点自己定，为了保证收敛，可以将给定的区间分成N等份，逐个试验是不是满足条件，我分的是100份

注意

这道题需要注意的地方还挺多的

考虑分母为0的情况，如果分母为0，就可以认为不符合，直接取下一个初值重新来
考虑-0的情况，这应该是OJ的问题了，不过，自己稍微加一下没啥
初始点的选择，不能随便取一个或者几个，可以取几十个几百个，我取100个迭代也没有超时

点评

坑的地方挺多的，自己也花了挺多的时间

还有，modified-newton method，真的很不明显，一开始在PPT和课本上都没有找到，复习的时候再课本的某段文字里发现了，不过可能modified-newton本来就不是一个专有名词吧

Q3. There is No Free Lunch

题目

6-3 There is No Free Lunch （40 分）

One day, CYLL found an interesting piece of commercial from newspaper: the Cyber-restaurant was offering a kind of “Lunch Special” which was said that one could “buy one get two for free”. That is, if you buy one of the dishes on their menu, denoted by $di$ with price $p_i$, you may get the two neighboring dishes $d{i - 1}$ and $d{i+1}$ for free! If you pick up $d_1$, then you may get $d_2$ and the last one $d_n$ for free, and if you choose the last one $d_n$, you may get $d{n - 1}$ and $d_1$ for free.

However, after investigation CYLL realized that there was no free lunch at all. The price $p_i$ of the i-th dish was actually calculated by adding up twice the cost $c_i$ of the dish and half of the costs of the two “free” dishes. Now given all the prices on the menu, you are asked to help CYLL find the cost of each of the dishes.

Format of function:

1	void Price( int n, double p[] );

where int n satisfies that 2< n ≤10000; double p[] is passed into the function as an array of prices $p_i$ and then will be returned storing the costs of the dishes.

Sample program of judge:

#include <stdio.h>
#define Max_size 10000 /* max number of dishes */
void Price( int n, double p[] );
int main()
{
    int n, i;
    double p[Max_size];
    scanf("%d", &n);
    for (i=0; i<n; i++) 
        scanf("%lf", &p[i]);
    Price(n, p);
    for (i=0; i<n; i++)
        printf("%.2f ", p[i]);
    printf("\n");
    return 0;
}
/* Your function will be put here */

Sample Input:

1	12 23.64 17.39 12.77 16.62 10.67 14.85 12.68 26.90 28.30 15.59 37.99 23.18

Sample Output:

1	9.20 5.58 3.24 7.00 1.99 6.36 2.25 10.01 11.52 0.50 17.65 4.88

思路

由题意不难得出一个方程组，其系数矩阵很有特征

自己知道有两种思路：

一种是对最后一行和第一行做处理，其它的都是一个完美的三对角矩阵，这样就可以比较轻松的解出来

另一种是直接用Jacobi迭代或者Gauss-Siedel迭代解方程

题目的本意我想肯定是前者，但是我用的后者，因为简单粗暴

注意

一个特别需要注意的地方，如果用迭代法做，自己想开一个二维数组存系数矩阵，并不可以，这道题的Max_size给的很大，会直接段错误，我估计case里边也有很大的系数矩阵，就算动态分配估计也不可以吧..
不过由于这个矩阵比较特殊，主对角线上的值就是最大的，而且每行的和相同，所以可以轻松的推出Jacobi迭代中自己需要的信息，而且写起来更简单

点评

好吧，用了简单粗暴的方式解决了本来还需要推导一番的问题，如果先做了第四题，那么这道题应该很容易了hh

Q4. Compare Methods of Jacobi with Gauss-Seidel

题目

6-4 Compare Methods of Jacobi with Gauss-Seidel （30 分）

Use Jacobi and Gauss-Seidel methods to solve a given $n\times n$ linear system $A\vec{x} = \vec{b}$ with an initial approximation $\vec{x^0}$
Note: When checking each $a{ii}$, first scan downward for the entry with maximum absolute value ($a{ii}$ included). If that entry is non-zero, swap it to the diagonal. Otherwise if that entry is zero, scan upward for the entry with maximum absolute value. If that entry is non-zero, then add that row to the i-th row.

Format of functions:

1
2
3

int Jacobi( int n, double a[][MAX_SIZE], double b[], double x[], double TOL, int MAXN );
int Gauss_Seidel( int n, double a[][MAX_SIZE], double b[], double x[], double TOL, int MAXN );

where int n is the size of the matrix of which the entries are in the array double a[][MAX_SIZE] and MAX_SIZE is a constant defined by the judge; double b[] is for $\vec{b}$, double x[] passes in the initial approximation $\vec{x^0}$ and return the solution $\vec{x}$; double TOL is the tolerance for $\lVert \vec{x^{k+1}} - \vec{x^k} \rVert$; and finally int MAXN is the maximum number of iterations.

The function must return:

k if there is a solution found after k iterations;
0 if maximum number of iterations exceeded;
1 if the matrix has a zero column and hence no unique solution exists;
2 if there is no convergence, that is, there is an entry of $\vec{x^K}$that is out of the range [-bound, bound] where bound is a constant defined by the judge.

Sample program of judge:

#include <stdio.h>
#include <math.h>
#define MAX_SIZE 10
#define bound pow(2, 127)
#define ZERO 1e-9 /* X is considered to be 0 if |X|<ZERO */
enum bool { false = 0, true = 1 };
#define bool enum bool
int Jacobi( int n, double a[][MAX_SIZE], double b[], double x[], double TOL, int MAXN );
int Gauss_Seidel( int n, double a[][MAX_SIZE], double b[], double x[], double TOL, int MAXN );
int main()
{
    int n, MAXN, i, j, k;
    double a[MAX_SIZE][MAX_SIZE], b[MAX_SIZE], x[MAX_SIZE];
    double TOL;
    scanf("%d", &n);
    for (i=0; i<n; i++) {
        for (j=0; j<n; j++)
            scanf("%lf", &a[i][j]);
        scanf("%lf", &b[i]);
    }
    scanf("%lf %d", &TOL, &MAXN);
    printf("Result of Jacobi method:\n");
    for ( i=0; i<n; i++ )
        x[i] = 0.0;
    k = Jacobi( n, a, b, x, TOL, MAXN );
    switch ( k ) {
        case -2:
            printf("No convergence.\n");
            break;
        case -1: 
            printf("Matrix has a zero column.  No unique solution exists.\n");
            break;
        case 0:
            printf("Maximum number of iterations exceeded.\n");
            break;
        default:
            printf("no_iteration = %d\n", k);
            for ( j=0; j<n; j++ )
                printf("%.8f\n", x[j]);
            break;
    }
    printf("Result of Gauss-Seidel method:\n");
    for ( i=0; i<n; i++ )
        x[i] = 0.0;
    k = Gauss_Seidel( n, a, b, x, TOL, MAXN );
    switch ( k ) {
        case -2:
            printf("No convergence.\n");
            break;
        case -1: 
            printf("Matrix has a zero column.  No unique solution exists.\n");
            break;
        case 0:
            printf("Maximum number of iterations exceeded.\n");
            break;
        default:
            printf("no_iteration = %d\n", k);
            for ( j=0; j<n; j++ )
                printf("%.8f\n", x[j]);
            break;
    }
    return 0;
}
/* Your function will be put here */

Sample Input 1:

3
2	-1	1	-1
2	2	2	4
-1	-1	2	-5
0.000000001 1000

Sample Output 1:

Result of Jacobi method:
No convergence.
Result of Gauss-Seidel method:
no_iteration = 37
1.00000000
2.00000000
-1.00000000

Sample Input 2:

3
1	2	-2	7
1	1	1	2
2	2	1	5
0.000000001 1000

Sample Output 2:

Result of Jacobi method:
no_iteration = 232
1.00000000
2.00000000
-1.00000000
Result of Gauss-Seidel method:
No convergence.

Sample Input 3:

5
2 1 0 0 0 1
1 2 1 0 0 1
0 1 2 1 0 1
0 0 1 2 1 1
0 0 0 1 2 1
0.000000001 100

Sample Output 3:

Result of Jacobi method:
Maximum number of iterations exceeded.
Result of Gauss-Seidel method:
no_iteration = 65
0.50000000
0.00000000
0.50000000
0.00000000
0.50000000

思路

思路很清晰，两个算法也都给出了，就是把它实现一遍

课本上原模原样的算法抄上去就可以了

注意需要在之前按照题意对矩阵进行调整

注意

由于自己一遍过了，所以觉得没啥需要注意的，调整矩阵，题目中也提示到了

点评

第一道做出的题目，抄算法，所以觉得这道题没啥..

Q5. Approximating Eigenvalues

题目

6-5 Approximating Eigenvalues （40 分）

Approximate an eigenvalue and an associated eigenvector of $a$ given $n \times n$ matrix $A$ near $a$ given value $p$ and a nonzero vector $
\vec{x} = (x_1, \dots, x_n)^T$

Format of function:

1	int EigenV(int n, double a[][MAX_SIZE], double *lambda, double v[], double TOL, int MAXN);

where int is the size of the matrix of which the entries are in the array double a[][MAX_SIZE] and MAX_SIZE is a constant defined by the judge; double *lambda is passed into the function as an initial approximation $p$ of the eigenvalue and is supposed to be returned as a more accurate eigenvalue; double v[] is passed into the function as the initial vector $\vec{x}$ and is supposed to be returned as the associated eigenvector with unit $\lVert \cdot \rVert$ norm; double TOL is the tolerance for the eigenvalue; and finally int MAXN is the maximum number of iterations.

The function must return:

1 if there is a solution;
0 if maximum number of iterations exceeded;
−1 if p is just the accurate eigenvalue.

Sample program of judge:

#include <stdio.h>
#define MAX_SIZE 10
int EigenV(int n, double a[][MAX_SIZE], double *lambda, double v[], double TOL, int MAXN);
int main()
{
    int n, MAXN, m, i, j, k;
    double a[MAX_SIZE][MAX_SIZE], v[MAX_SIZE];
    double lambda, TOL;
    scanf("%d", &n);
    for (i=0; i<n; i++) 
        for (j=0; j<n; j++) 
            scanf("%lf", &a[i][j]);
    scanf("%lf %d", &TOL, &MAXN);
    scanf("%d", &m);
    for (i=0; i<m; i++) {
        scanf("%lf", &lambda);
        for (j=0; j<n; j++)
            scanf("%lf", &v[j]);
        switch (EigenV(n, a, &lambda, v, TOL, MAXN)) {
            case -1: 
                printf("%12.8f is an eigenvalue.\n", lambda );
                break;
            case 0:
                printf("Maximum number of iterations exceeded.\n");
                break;
            case 1:
                printf("%12.8f\n", lambda );
                for (k=0; k<n; k++)
                    printf("%12.8f ", v[k]);
                printf("\n");
                break;
        }
    }
    return 0;
}
/* Your function will be put here */

Sample Input 1:

3
1 2 3
2 3 4
3 4 5
0.0000000001 1000
1
-0.6 1 1 1

Sample Output 1:

1 2	-0.62347538 1.00000000 0.17206558 -0.65586885

Sample Input 2:

2
0 1
1 0
0.0000000001 10
2
1.0 1 1
100.0 1 0

Sample Output 2:

1 2	1.00000000 is an eigenvalue. Maximum number of iterations exceeded.

思路

近似特征值的算法用反幂法是可以的

我按照书上的算法走的

书上的算法Step1是给p赋值，本题目中给定了初值，直接用它就可以

需要解线性方程组，直接解似乎会超时，需要做一个LU Factorization然后之后的每次计算用L和U矩阵而不是原始的A矩阵，如果LU factorization做不下去，那就是题目中说的给定值就是初始值了

注意

有两个需要注意的地方

不要改掉原来的系数矩阵，因为每次迭代解方程的系数矩阵需要在原矩阵对角线上减去lambda，要在原来的矩阵上减，而不是一直减下去
课本上的算法有坑！至少英文版第七版上是不对的，Step8：$\mu = y_p$，而新的$p$是在Step9算的，所以Step8应该在Step9之后！

点评

虽然这两个坑好像很多人踩，但我都没有踩到，也算是比较幸运的吧

代码量还是挺大的，比较烦，需要LU分解和解方程，不过都有算法，自己想想也不复杂

Q6. Cubic Spline

题目

6-6 Cubic Spline （40 分）

Construct the cubic spline interpolant $S$ for the function $f$, defined at points $x_0 < x_1 < \dots < x_n$, satisfying some given boundary conditions. Partition a given interval into $m$ equal-length subintervals, and approximate the function values at the endpoints of these subintervals.

Format of functions:

1
2
3

void Cubic_Spline(int n, double x[], double f[], int Type, double s0, double sn, double a[], double b[], double c[], double d[]);
double S( double t, double Fmax, int n, double x[], double a[], double b[], double c[], double d[] );

The function Cubic_Spline is for calculating the set of coefficients of
$S(x)$ where $S(x) = S^{[j]}(x) = aj + b_j(x - x{j - 1}) + cj(x - x{j - 1})^2 + dj(x - x{j - 2})^3$. The function S is for obtaining an approximation of $f(t)$ using the spline function $S(t)$. The parameters are defined as the following:

int n is the number of interpolating points−1;
double x[] contains n+1 distinct real numbers $x_0 < x_1 < \dots < x_n$;
double f[] contains n+1 real numbers $f(x_0), f(x_1), \dots ,f(x_n)$;
int Type is the type of boundary conditions, where 1 corresponds to the clamped boundary condition $S’(x_0) = s_0, S’(x_n) = s_n$ and 2 corresponds to the natural boundary condition $S’’(x_0) = s_0, S’’(x_n) = s_n$;
double s0 and double sn are $s_0$ and $s_n$ respectively;
double a[], b[], c[], and d[] contain the set of coefficients of $S(x)$;
double t is an endpoint of a subinterval;

and finally double Fmax is the default value that $S(t)$ will assume if t is out of the range $[x_0, x_n]$

Sample program of judge:

#include <stdio.h>
#define MAX_N 10
void Cubic_Spline(int n, double x[], double f[], int Type, double s0, double sn, double a[], double b[], double c[], double d[]);
double S( double t, double Fmax, int n, double x[], double a[], double b[], double c[], double d[] );
int main()
{
    int n, Type, m, i;
    double x[MAX_N], f[MAX_N], a[MAX_N], b[MAX_N], c[MAX_N], d[MAX_N];
    double s0, sn, Fmax, t0, tm, h, t;
    scanf("%d", &n);
    for (i=0; i<=n; i++) 
        scanf("%lf", &x[i]);
    for (i=0; i<=n; i++) 
        scanf("%lf", &f[i]);
    scanf("%d %lf %lf %lf", &Type, &s0, &sn, &Fmax);
    Cubic_Spline(n, x, f, Type, s0, sn, a, b, c, d);
    for (i=1; i<=n; i++)
        printf("%12.8e %12.8e %12.8e %12.8e \n", a[i], b[i], c[i], d[i]);
    scanf("%lf %lf %d", &t0, &tm, &m);
    h = (tm-t0)/(double)m;
    for (i=0; i<=m; i++) {
        t = t0+h*(double)i;
        printf("f(%12.8e) = %12.8e\n", t, S(t, Fmax, n, x, a, b, c, d));
    }
    return 0;
}
/* Your functions will be put here */

Sample Input:

2
0.0 1.0 2.0
0.0 1.0 2.0
1 1.0 1.0 0.0
0.0 3.0 2

Sample Output:

0.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 
1.00000000e+00 1.00000000e+00 0.00000000e+00 0.00000000e+00 
f(0.00000000e+00) = 0.00000000e+00
f(1.50000000e+00) = 1.50000000e+00
f(3.00000000e+00) = 0.00000000e+00

思路

算法名称已经告诉了，书上也有现成的算法

有一个区别在于，本题目的Natural Cubic Spline并不是两端的二阶导为0，而是等于s0和sn. 这个需要对原算法做出一定的调整，其实就是推导出来的系数矩阵中的系数和右侧的首尾的值改一下就可以了..

注意

很坑的一点，很容易忽视的一点！！

如果发现自己总是2、3点过不了，有时候不一定是Spline算法的问题，而是自己写的S函数有问题！！我当时2、3点过不了，就是S函数的问题，注意区间、边界值之类的，反正多试几次应该就有了~

点评

唉，最后两个点，死活过不了，提交次数最多的一道题了！！

没想到啊没想到，竟然是S的问题而不是算法的问题，气死我了！！

Q7. Orthogonal Polynomials Approximation

题目

6-7 Orthogonal Polynomials Approximation （40 分）

Given a function $f$ and a set of $m>0$ distinct points $x1 < x_2 < \dots < x_m$. You are supposed to write a function to approximate $f$ by an orthogonal polynomial using the exact function values at the given m points with a weight $w(x_i)$ assigned to each point $x_i$. The total error $err =
sum{i=1}^{m}w(x_i)|f(x_i) - P_n(x_i)|^2$ must be no larger than a given tolerance.

Format of function:

1	int OPA( double (f)(double t), int m, double x[], double w[], double c[], double eps );

where the function pointer double (*f)(double t) defines the function $f$; int m is the number of points; double x[] contains points $x_1 < x_2 < \dots < x_m$; double w[] contains the values of a weight function at the given points x[]; double c[] contains the coefficients $c_1 < c_2 < \dots < c_m$ of the approximation polynomial $P_n(x) = c_0 + c_1x + \dots + c_nx^n$; double *eps is passed into the function as the tolerance for the error, and is supposed to be returned as the value of error. The function OPA is supposed to return the degree of the approximation polynomial.

Note: a constant Max_n is defined so that if the total error is still not small enough when n = Max_n, simply output the result obtained when n = Max_n.

Sample program of judge:

#include <stdio.h>
#include <math.h>
#define MAX_m 200
#define MAX_n 5
double f1(double x)
{
    return sin(x);
}
double f2(double x)
{
    return exp(x);
}
int OPA( double (*f)(double t), int m, double x[], double w[], double c[], double *eps );
void print_results( int n, double c[], double eps)
{	
    int i;
    printf("%d\n", n);
    for (i=0; i<=n; i++)
        printf("%12.4e ", c[i]);
    printf("\n");
    printf("error = %9.2e\n", eps);
    printf("\n");
}
int main()
{
    int m, i, n;
    double x[MAX_m], w[MAX_m], c[MAX_n+1], eps;
    m = 90;
    for (i=0; i<m; i++) {
        x[i] = 3.1415926535897932 * (double)(i+1) / 180.0;
        w[i] = 1.0;
    }
    eps = 0.001;
    n = OPA(f1, m, x, w, c, &eps);
    print_results(n, c, eps);
    m = 200;
    for (i=0; i<m; i++) {
        x[i] = 0.01*(double)i;
        w[i] = 1.0;
    }
    eps = 0.001;
    n = OPA(f2, m, x, w, c, &eps);
    print_results(n, c, eps);
    return 0;
}
/* Your function will be put here */

Sample Output:

3
 -2.5301e-03   1.0287e+00  -7.2279e-02  -1.1287e-01 
error =  6.33e-05
4
  1.0025e+00   9.6180e-01   6.2900e-01   7.0907e-03   1.1792e-01 
error =  1.62e-04

思路

算法的话，按照PPT上的算法即可，不需要做太大的改动，注意需要求内积，自己写个内积函数

注意

很需要注意的一点，就是题目要求的是标准系数，而不是算法求出来每个项的系数，这是很需要注意的！！

当你算出来error是一样的，但是系数全部不一样的时候，绝大多数都是这个错误！

点评

课本上没有这个算法，自己搞懂原理之后，写出来一个，样例是可以过的，但是总是有些点过不了，改用PPT的算法就过了，差不多都是一样的，除了error的计算，它是一步步减的，我是最后算的，可是本质上是一样的，唉，不知道为什么，不过AC了，就不去想它了。

Q8. Shape Roof

题目

6-8 Shape Roof （40 分）

The kind of roof shown in Figure 1 is shaped from plain flat rectangular plastic board in Figure 2.

Figure 1

Figure 2

The transection of the roof is a sine curve which fits the function $y(x)=lsin(lx)$ in centimeters. Given the length of the roof, your task is to calculate the length of the flat board needed to shape the roof.

Format of function:

1	double Integral(double a, double b, double (*f)(double x, double y, double z), double eps, double l, double t);

where the function Integral returns the value of a definit integral of a given function pointed by f over a given interval from double a to double b. The integrant function f will be defined in the sample program of judge. Other parameters are: double eps, the accuracy of the resulting value of the integral; double l and double t, the altitude and the frequency of the sine curve, respectively.

Note: the length of the flat board which Integral returns must be in meters.

Sample program of judge:

#include <stdio.h>
#include <math.h>
double f0( double x, double l, double t )
{
    return sqrt(1.0+l*l*t*t*cos(t*x)*cos(t*x));
}
double Integral(double a, double b, double (*f)(double x, double y, double z), 
double eps, double l, double t);
int main()
{
    double a=0.0, b, eps=0.005, l, t;
    scanf("%lf %lf %lf", &l, &b, &t);
    printf("%.2f\n", Integral(a, b, f0, eps, l, t));
    return 0;
}
/* Your function will be put here */

Sample Input:

2 100 1

Sample Output:

1.68

思路

就是一个积分，甚至积分公式的函数都给出来了，直接代入值就可以了

分成足够多的小段，矩形累加起来是可以过的

注意

没啥注意的，试几次总会过得

点评

笑死我了，题目本意肯定不是这样，不过，反正，四行代码，过了，就那样了，可能是今年的case要求不是很高吧，分了4000000段，都没有超时~

总结

八道题，除了第一道觉得很有意思之外，其它的基本上都有现成的算法，所以做过之后觉得还好..

感觉这种小数点后十几位的要求，有一定的玄学因素在里边，可能同样的算法，步骤顺序不同，结果就有差异..

不过这八道题做起来还是很有价值的，不是像PAT、LeetCode那样的一个算法知识点，需要的是对数学算法的实现，还是比较有趣的~