数值分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
1. (程序题)连续开方再连续平方
用C++语言的开方函数sqrt，连续K次开方，然后再连续K次平方，用双精度计算，结果肯定回不到原始值。
如果要求相差必须大于eps，请问最小的K是多少。注意，必须用C++，否则有可能过不了。
输入格式：
一行，一个整数N和一个小于0.1的实数（格式，比如1.3e-5）。
输出格式：
N-1行，每行两个整数，第一个数是2到N，从小到大，第二个整数为对应的K，比如第三行，一个整数是4，如果eps=1.0e-10，K=18，就是连续18次开方再连续18次平方，与4的误差才超过1.0e-10。
每行两个整数之间有3个空格。
输入样例：
20 1.e-10
输出样例：
2   19
3   20
4   18
5   18
6   18
7   17
8   18
9   17
10   19
11   20
12   19
13   16
14   18
15   16
16   18
17   19
18   17
19   19
20   18
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
#include<math.h>
#include<iostream>
using namespace std;

/*
	用C++语言的开方函数sqrt，连续K次开方，然后再连续K次平方，
	用双精度计算，结果肯定回不到原始值。如果要求相差必须大于eps，请问最小的K是多少
*/

int main(){

	int N=0;
	double n=0;
	cin>>N>>n;
	
	for(int i=2;i<=N;i++){	//进行2到N的最小k计算 
		
		int flag=0;	//开方次数标识 
		double a=i;	//双精度i用于开方 
		
		while(1){	//不断开方，直到找到最小k 
			
			a=sqrt(a);
			flag++;
			double temp=a;
			for(int k=flag;k>0;k--){
				
				temp*=temp;
				
			}
		
			if(fabs(temp-i)>n){
				break;
			}
			
		}
		cout<<i<<"   "<<flag<<endl;
	}
	
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
拉格朗日插值：
给定N个点（xi, yi），请用拉格朗日插值，计算y(t)，并计算基函数li(t)的绝对值之和
输入格式（数之间有一个空格）：
第一行两个整数N和M（1<N<20, 0<M<10），两个整数之间有一个空格
第二行2N个实数，x0, y0, x1, y1,,,,。用于插值的数据点
第三行M个实数，t0, t1, t2,,,,。需要计算的插值点
输出格式(精确到小数点后3位）：
M行，每行三个数，t, y(t)和基函数li(t)的绝对值之和（数之间有一个空格）
输入样例：
4 1
3.646 3.849 0.573 -0.298 1.774 -1.571 2.897 1.374
0.777
输出样例：
0.777 -1.008 1.597
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
#include<stdio.h>
#include<iostream>
#include <malloc.h>
#include<math.h>
using namespace std;
#define maxsize 21


double* radixfunctions(double t,int N,double x[]) {

	double* radix = (double*)malloc(maxsize * sizeof(double));
	for (int i = 0; i < N; i++) {
		radix[i] = 1;
		for (int j = 0; j < N; j++) {
			if (i != j) {
				radix[i] *= (x[i] - x[j]);
			}
		}
	}

	double* radixResult = (double*)malloc(maxsize * sizeof(double));

	for (int i = 0; i < N; i++) {
		radixResult[i] = 1;
		for (int j = 0; j < N; j++) {
			if (i!=j) {
				radixResult[i] *= (t - x[j]);
			}
		}
		radixResult[i] /= radix[i];

	}

	return radixResult;
}

int main()
{
	int N = 0;
	int M = 0;

	double x[maxsize] = { 0 };
	double y[maxsize] = { 0 };
	double t[maxsize] = { 0 };
	cin >> N >> M;
	int flag = 0;
	char d;
	for (int k = 0; k < N; k++) {
		cin >> x[flag] >> y[flag];
		flag++;
	}
	
	flag = 0;
	for (int j = 0; j < M; j++) {
		cin >> t[flag];
		flag++;
	}

	for (int i = 0; i < M; i++) {
		
		double* radixresult = radixfunctions(t[i], N, x);
		double Aresult = 0;
		double Bresult = 0;
		for (int j = 0; j < N; j++) {

			Aresult += radixresult[j] * y[j];
			Bresult += fabs(radixresult[j]);

		}
		printf("%.3lf %.3lf %.3lf\n", t[i], Aresult, Bresult);

	}
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
牛顿插值：

给定N个点（xi, yi），请用牛顿插值，输出各阶均差，并计算y(t)
输入格式：
第一行两个整数N和M（N<20, M<10），两个整数之间有一个空格
第二行2N个实数，x0, y0, x1, y1,,,,。用于插值的数据点
第三行M个实数，t0, t1, t2,,,,。需要计算的插值点

输出格式（精确到小数点后3位）：
第一行，N个实数，从0阶均差到N-1阶均差，数之间有一个空格
第二行，N个实数，最后整理成多项式的形式，系数a0,a1,,,an-1
第三行，M个实数，用牛顿插值公式计算的M个插值值，即用t0, t1, t2,,,,计算的结果。

输入样例：
3 1
2.643 0.53 2.532 0.237 1.594 -1.236
2.161

输出样例：
0.530 2.640 1.019
0.375 -2.635 1.019
-0.560
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
#include<stdio.h>
#include<iostream>
#include <malloc.h>
#include<math.h>
using namespace std;
#define maxsize 21


double** average(double y[], double x[], int N) {

	//牛顿插值法 求均差
	double** p = new double* [N];
	for (int i = 0; i < N; i++) {

		p[i] = new double[N];

	}

	for (int i = 0; i < N; i++) {

		p[0][i] = y[i];

	}

	for (int i = 1; i < N; i++) {
		for (int j = i; j < N; j++) {

			p[i][j] = (p[i - 1][j] - p[i - 1][j - 1]) / (x[j] - x[j - i]);

		}
	}
	return p;

}

double NewtonPoly(double y[],double x[], double t, int N) {
	//求Newton插值多项式
	double** Average = average(y, x, N);
	double result = Average[0][0];
	for (int i = 1; i < N; i++) {

		double temp = Average[i][i];

		for (int j = 0; j <i; j++) {

			temp = temp * (t - x[j]);
		}

		result += temp;
	}
	return result;
}


double* PolyFactors(double y[], double x[], int N) {
	//求Newton多项式的系数

	double** Average = average(y, x, N);
	double result[maxsize][maxsize] = { 0 };
	result[0][0] = 1;
	for (int j = 1; j < N; j++) {
		for (int k = 0; k <= 1; k++) {
			if(k==0)result[j][k] = 1;
			if (k == 1)result[j][k] = -x[j-1];
			
		}
	}

	int flag = 1;
	for (int i = 2; i < N; i++) {
		
		double* temp = (double*)malloc(sizeof(double) * flag);
		for (int k = 0; k < 2; k++) {
			temp[k] = result[i][k];
		}
		
		for (int m = 0; m <= 1; m++) {
			for (int j = 0; j <= flag; j++) {

				if (m == 0) {
					if (j == 0)result[i][j] *= 1;
					else {
						result[i][j] = 1 * result[i - 1][j];
					}
				}
				if (m == 1) {
					result[i][j + 1] += temp[1] * result[i - 1][j];
				}


			}

		}

		flag++;
	}

	double BeautifulResult[maxsize] = { 0 };


	for (int j = 0; j<N; j++) {

		int flag = 0;

		for (int i = j; i < N ; i++) {

			BeautifulResult[j] += Average[i][i] * result[i][flag];
			flag++;
		}

	}
	
	return BeautifulResult;
}

int main() {

	int N = 0;
	int M = 0;

	double x[maxsize] = { 0 };
	double y[maxsize] = { 0 };
	double t[maxsize] = { 0 };
	cin >> N >> M;
	int flag = 0;

	for (int k = 0; k < N; k++) {
		cin >> x[flag] >> y[flag];
		flag++;
	}

	flag = 0;
	for (int j = 0; j < M; j++) {
		cin >> t[flag];
		flag++;
	}
	
	double** Average = average(y, x, N);
	double* polyFator = PolyFactors(y,x,N);


	double result[maxsize] = { 0 };
	for (int i = 0; i < M; i++) {

		result[i] = NewtonPoly(y, x, t[i], N);

	}
	for (int i = 0; i < N; i++) {
		printf("%.3lf ", Average[i][i]);
		if (i == N - 1)printf("\n");
	}
	for (int i = 0; i < N; i++) {
		printf("%.3lf ", polyFator[i]);
		if (i == N - 1)printf("\n");
	}
	for (int i = 0; i < M; i++) {
		printf("%.3lf ", result[i]);
		if (i == N - 1)printf("\n");
	}

}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
1. (程序题)
三次样条插值

给定N个点（xi, yi），请用三次样条插值，计算y(t)，两端边界条件是二次导数为0.

输入格式：
第一行两个整数N和M（3<N<20, 0<M<10），两个整数之间有一个空格
第二行2N个实数，x0, y0, x1, y1,,,,。用于插值的数据点，数据之间有一个空格，xi从小到大输入。
第三行M个实数，t0, t1, t2,,,,。需要计算的插值点，数据之间有一个空格

输出格式：
M行，每行2个实数，ti及用三次样条计算的值y(ti)，数之间有一个空格，精确到小数点后三位，后面的0也要输出。输出顺序按原始顺序。



输入样例：
4 7
1.382 -2.085 1.722 -0.74 1.973 -0.655 3.026 0.552
1.706 1.613 2.06 2.756 2.365 1.885 2.941
输出样例：

1.706 -0.771
1.613 -1.036
2.060 -0.655
2.756 0.079
2.365 -0.466
1.885 -0.638
2.941 0.400
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
#include<iostream>
#include<malloc.h>
using namespace std;
#define maxsize 21

//求系数Mi
void spline(double x[],double y[],double M[],int N){

    // 先求出α[i]
    double *a = (double *)malloc(N * sizeof(double));
    a[0] = 1;
    a[N - 1] = 0;
    for (int i = 1; i < N-1;i++){
        a[i] = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
    }

    // 再求出β[i]
    double* bet = (double*)malloc(N * sizeof(double));
    bet[0] = 3 * (y[1] - y[0]) / (x[1] - x[0]);
    bet[N - 1] = 3 * (y[N - 1] - y[N - 2]) / (x[N - 1] - x[N - 2]);

    for (int i = 1; i < N - 1;i++){
        bet[i] = 3 * ((1 - a[i]) * ((y[i] - y[i - 1]) / (x[i] - x[i - 1])) +
                      a[i] * (y[i + 1] - y[i]) / (x[i + 1] - x[i]));
    }

    //接着求A[i],B[i]
    double* A = (double*)malloc((N-1) * sizeof(double));
    double* B = (double*)malloc((N-1)* sizeof(double));
    A[0] = -a[0] / 2;
    B[0] = bet[0] / 2;

    for (int i = 1; i < N - 1;i++){

        A[i] = -a[i] / (2 + (1 - a[i]) * A[i - 1]);
        B[i] = (bet[i] - (1 - a[i]) * B[i - 1]) / (2 + (1 - a[i]) * A[i - 1]);

    }

    //终于可以求M[i]啦
    M[N - 1] = (bet[N - 1] - (1 - a[N - 1]) * B[N - 2]) /
               (2 + (1 - a[N - 1]) * A[N - 2]);

    for (int i = N-2; i >=0;i--){

        M[i] = A[i] * M[i + 1] + B[i];

    }

}

int findindex(double x[],double t,int N){

    int index = 0;
    int i = 0;
    for ( i = 0; i < N;i++){
        if(t>x[i])continue;
        if(t<=x[i]){
            break;
        }
    }
    index = i-1;
    return index;
}


int main(){

    //处理输入,顺便排个序
    int N = 0;
    int M = 0;

    double x[maxsize] = {0};
    double y[maxsize] = {0};
    double t[maxsize] = {0};
    cin >> N >> M;
    int flag = 0;

    for (int k = 0; k < N; k++) {
        cin >> x[flag] >> y[flag];
        flag++;
    }

    flag = 0;
    for (int j = 0; j < M; j++) {
        cin >> t[flag];
        flag++;
    }

    double result[maxsize] = {0};
    double m[maxsize]={0};
    spline(x, y, m, N);
    for (int i = 0; i < M;i++){
        int dex = findindex(x, t[i], N);
        result[i] =
            (1 + 2 * (t[i] - x[dex]) / (x[dex + 1] - x[dex])) *
                ((t[i] - x[dex + 1]) / (x[dex] - x[dex + 1])) *
                ((t[i] - x[dex + 1]) / (x[dex] - x[dex + 1])) * y[dex] +
            (1 + 2 * (t[i] - x[dex + 1]) / (x[dex] - x[dex + 1])) *
                ((t[i] - x[dex]) / (x[dex + 1] - x[dex])) *
                ((t[i] - x[dex]) / (x[dex + 1] - x[dex])) * y[dex + 1] +
            (t[i] - x[dex]) * ((t[i] - x[dex + 1]) / (x[dex] - x[dex + 1])) *
                ((t[i] - x[dex + 1]) / (x[dex] - x[dex + 1])) * m[dex] +
            (t[i] - x[dex + 1]) * ((t[i] - x[dex]) / (x[dex + 1] - x[dex])) *
                ((t[i] - x[dex]) / (x[dex + 1] - x[dex])) * m[dex + 1];
    }

    for (int i = 0; i < M;i++){
        printf("%.3f %.3f\n", t[i], result[i]);
    }
    system("pause");
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
线性拟合

给定n个数据点(xi，yi)，用线性y=a+bx，进行回归，输出a、b、R2
设所有yi的平均值是ya, （yi-ya）的平方之和为tss，（a+bxi-ya）的平方之和为ess
R2=ess/tss。

输入格式：
第一行，整数n，n<100
第二行，n个实数，xi，数据之间有空格
第三行，n个实数，yi，数据之间有空格
输出格式：
一行，a、b、R2，保留小数点后三位


输入样例：
5
3 4 2 1 2
15 18 11 8 13


输出样例：
5.154 3.269 0.958
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
#include<iostream>
#include<stdio.h>
#include<malloc.h>
#include<math.h>
using namespace std;
#define maxsize 101

//线性拟合,此处求系数a,b,R2
void linefit(double x[],double y[] ,int N,double result[]){

    double temp1 = 0;
    double temp2 = 0;

    //求L11，和L1y
    for (int i = 0; i < N;i++){
        temp1 += pow(x[i],2);
        temp2 += x[i];
    }
    temp2 = pow(temp2, 2) / N;
    temp1 -= temp2;//这里放的是L11
    temp2 = 0;
    double temp3 = 0;
    double temp4 = 0;
    for (int i = 0; i < N;i++){
        temp2 += y[i] * x[i];
        temp3 += x[i];
        temp4 += y[i];
    }
    temp2 = temp2 - temp3 * temp4 / N;//这里放的是L1y
    temp3 /= N;//x1的平均值x1a
    temp4 /= N;//y的平均值ya

    result[1] = temp2 / temp1;//a1
    result[0] = temp4 - (temp2 / temp1) * temp3;//a0

    //接下来求R2
    double tss = 0;
    double ess = 0;
    for (int i = 0; i < N;i++){
        tss += pow(y[i] - temp4, 2);
        ess += pow(result[0] + result[1] * x[i] - temp4, 2);
    }
    result[2] = ess / tss;
}


int main(){
    int N = 0;
    double x[maxsize] = {0};
    double y[maxsize] = {0};
    double result[maxsize] = {0};
    cin >> N;
    for (int i = 0; i < N;i++){
        cin >> x[i];
    }
    for (int i = 0; i < N;i++){
        cin >> y[i];
    }
    linefit(x, y, N, result);
    for (int i = 0; i < 3;i++){
        printf("%.3f", result[i]);
        if(i<2)
            printf(" ");
    }
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
两棱柱公共部分体积
问题描述：
有一个侧棱与z轴平行的棱柱P1和一个侧棱与y轴平行的棱柱P2。它们都向两端无限延伸，底面分别是包含M个顶点和N个顶点的凸多边形，其中第i个顶点的坐标分别是(X1i，Y1i)和(X2i，Y2i)。请计算这两个棱柱公共部分的体积。
限制条件：3≤M, N≤100，-100≤X1i,Y1i,X2i,Y2i≤100

输入格式：
第一行，m
接着m个点坐标，一个点一行(都是实数）

X11 Y11
X12 Y12
……
X1m Y1m

接着一个整数，n
接着n个点坐标，一个点一行(都是实数）
X21 Y21
X22 Y22
……
X2n Y2n
这些顶点位置的顺序在xy平面或xz平面上按逆时针顺序给出，且多边形都是凸的。

输出格式：
输出两个棱镜P1和P2相交的体积，精确到小数点后3位。

输入样例：
4 3
7 2
3 3
0 2
3 1
4 2
0 1
8 1
输出样例：
4.708
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
逐次分半求积：

给定误差eps=1.e-6，被积函数（1+x*x*sin(px)*sin(px))/(1+q*x*x)，积分区间[a, b]，请用Romberg（龙贝格）求积法计算数值积分，输出数值积分区间数和积分结果。

Romberg（龙贝格）求积法请见下面算法

输入格式：

第一行四个实数，p, q, a, b，两个整数之间有一个空格

输出格式（实数精确到小数点后6位）：

第一行，一个整数和一个实数，Romberg法的积分小区间数（即积分点减1）和数值积分结果。

第二行，作为对比，输出逐次分半辛普森法的结果，一个整数和一个实数，积分小区间数和数值积分结果。


输入样例：

1.696 2.924 0.168 5.675



输出样例：

64 1.520629

128 1.520626
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
#include<iostream>
#include<cmath>
#include<vector>
using namespace std;
#define eps 1.e-6
 
//用麦克劳林公式,计算sin（px）的近似值狗屎
double maclaurinSin(double x,double p){
    double result = 0;
    double use = 0;
    for (int i = 1;true;i+=2){
        use = pow(x,i);
        for (int j = 1; j <= i; j++) {
            use /= j;
        }
        result += pow(-1.0, i / 2)*use;
        if(use<=1.e-10){
            return result;
            break;
        }
    }
}

//计算函数的值
double functionValue(double x,double p,double q){
    double result = 0;
    double temp = x * p;
    result =
        (1 + x * x * sin(temp) * sin(temp)) / (1 + q * x * x);
    return result;
}

void romberg(double a,double b,double p,double q){

    int k = 0;
    int n = 1;
    double h = (b - a) / 2;
    double Rom[60] = {0};
    int flag = 0;
    double temp = h * (functionValue(a, p, q) + functionValue(b, p, q));
    Rom[flag] = temp;
    flag++;
    for (k = 1; k <= 9;k++){
        int m;
        double F = 0;
        for (int i = 1; i <= n; i++) {
            F += functionValue(a + (2 * i - 1) * h, p, q);
        }
        for (m=0; m <= k;m++){
            double need = 0;
            if(m==0){
                need = Rom[flag-k]/2+h*F;
                Rom[flag] = need;
                flag++;
            }else{
                need =(pow(4,m)*Rom[flag-1]-Rom[flag-1-k])/(pow(4,m)-1);
                Rom[flag] = need;
                flag++;
            }

        }
        if(fabs(Rom[flag-1]-Rom[flag-2])<eps||k==9){
            printf("%.0f %.6f",pow(2,k), Rom[flag-1]);
            break;
        }
        h /= 2;
        n *= 2;
    }
}

void simpson(double a, double b, double p, double q){
    double F1 = 0;
    double F2 = 0;
    int n = 2;
    double h = (b - a) / 4;
    F1 = functionValue(a, p, q) + functionValue(b, p, q);
    F2 = functionValue((a + b) / 2, p, q);
    double S = (b - a) * (F1 + 4 * F2) / 6;
    for (; true; ){
        double F3 = 0;
        for (int i = 1; i <= n;i++){
            F3 += functionValue(a + (2 * i - 1) * h, p, q);
        }
        double s = h * (F1 + 2 * F2 + 4 * F3) / 3;
        if(fabs(s-S)<15*eps){
            printf("%d %f", 2*n, s);
            break;
        }else{
            h /= 2;
            n *= 2;
            F2 += F3;
            S = s;
        }
    }
}

    int main() {
    double p = 0;
    double q = 0;
    double a = 0;
    double b = 0;
    cin >> p >> q >> a >> b;
    romberg(a, b, p, q);
    cout << endl;
    simpson(a, b, p, q);
    system("pause");
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
1. (程序题)
矩阵条件数：
给定矩阵A。求A的条件数，范数采用无穷范数。

输入格式：

第一行一个整数n，矩阵的行数(n<20)

接着n行，每行是矩阵的一行数据，按行顺序。每行n个实数。



输出格式：

一个实数，矩阵A的条件数，精确到小数点后4位。



输入样例：
3
102 60 89
38 45 91
97 50 108


输出样例：
34.6863
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
#include<iostream>

#include<stdio.h>

using namespace std;

#define MaxSize 21

#define EPS 1.0e-9





//高斯约旦法求矩阵的逆，采用列主元消去法

void matrix_gauss(double A[][MaxSize],int N) {



 for (int pivot = 0; pivot < N; pivot++) {

 double B[MaxSize] = { 0 };

 B[pivot] = 1.0;

 if (pivot < N - 1) {

 int maxrow = pivot;

 double maxval = abs(A[pivot][pivot]);

 for (int row = pivot + 1; row > N; row--) {

 double val = abs(A[row][pivot]);

 if (val > maxval) {

 maxval = val;

 maxrow = row;

 }

 }

 if (maxrow != pivot) {

 for (int k = 0; k < N; k++) {

 swap(A[maxrow][k], A[pivot][k]);

 }

 }

 }



 double coef = 1.0 / A[pivot][pivot];

 if (abs(coef) > EPS) {

 for (int col = 0; col < N; col++) {

 A[pivot][col] *= coef;

 }

 B[pivot] *= coef;

 }



 for (int row = 0; row < N; row++) {

 if (row == pivot)continue;

 coef = 1.0 * A[row][pivot];

 if (abs(coef) > EPS) {

 for (int col = 0; col < N; col++) {

 A[row][col] -= coef * A[pivot][col];

 }

 B[row] -= coef * B[pivot];

 }

 }



 for (int row = 0; row < N; row++) {

 A[row][pivot] = B[row];

 }

 }



}



int main() {

 int N;

 cin >> N;

 double MaxValue1 = 0.0;

 double MaxValue2 = 0.0;

 double matr[MaxSize][MaxSize] = { 0 };



 for (int i = 0; i < N; i++) {

 for (int j = 0; j < N; j++) {

 cin >> matr[i][j];

 }

 }



 for (int i = 0; i < N; i++) {

 double temp = 0.0;

 for (int j = 0; j < N; j++) {

 temp += abs(matr[i][j]);

 }

 if (temp > MaxValue1) {

 MaxValue1 = temp;

 }

 }

 matrix_gauss(matr, N);

 for (int i = 0; i < N; i++) {

 double temp = 0.0;

 for (int j = 0; j < N; j++) {

 temp += abs(matr[i][j]);

 }

 if (temp > MaxValue2) {

 MaxValue2 = temp;

 }

 }

 double result = MaxValue1 * MaxValue2;

 printf("%.4f", result);

}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
#include<iostream>
#include<stdio.h>
#include<math.h>
using namespace std;
#define EPS 1.e-6

void Jacobi(double* a[], int N) {
	double* oldResult = (double*)malloc(N * sizeof(double));
	double* NewResult = (double*)malloc(N * sizeof(double));

	//用于判断收敛误差
	double BNRM = 0.0;
	for (int i = 0; i < N; i++) {
		BNRM = BNRM + a[i][N] * a[i][N];
	}
	BNRM = sqrt(BNRM);

	if (oldResult != NULL && NewResult != NULL) {
		//初始值设为0
		for (int i = 0; i < N; i++) {
			NewResult[i] = 0.0;
		}
		for (int iter = 0; iter < N * 100; iter++) {

			for (int i = 0; i < N; i++) {
				oldResult[i] = NewResult[i];
			}

			for (int i = 0; i < N; i++) {
				double RESID = a[i][N];
				for (int j = 0; j < N; j++) {
					if (j != i) {
						RESID = RESID - a[i][j] * oldResult[j];
					}
				}
				NewResult[i] = RESID / a[i][i];
			}

			//收束判定
			double VAL = 0.0;
			for (int i = 0; i < N; i++) {
				VAL = VAL + pow((NewResult[i] - oldResult[i]), 2);
			}
			VAL = sqrt(VAL);
			if (VAL / BNRM <= EPS) {
				break;
			}
		}
	}

	
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
1. (程序题) 构造迭代公式
给定一个线性方程组Ax=b，请输出迭代公式x=Bx+f的B和f



输入格式：

先输入一个整数n，A矩阵的阶（n<10）

按增广矩阵的格式输入，共n行，每行n+1个实数，分别是增广矩阵的一行数据，数据之间有一个空格。



输出格式：

B和f按增广矩阵的形式输出，即共n行，每行n+1个实数，分别是增广矩阵的一行数据，数据之间有一个空格。

共输出三种格式的结果：Jacobi、Gauss-Seidel及松弛影子0.7的超松弛迭代法。

每种迭代结果之后输出一个空行，精确到小数点后4位。

注意：每行后面没有空格



输入样例：

3

167 17 70 1742

62 245 88 2012

67 38 176 324



输出样例：

0.0000  -0.1018  -0.4192  10.4311

-0.2531  0.0000  -0.3592  8.2122

-0.3807  -0.2159  0.0000  1.8409



0.0000  -0.1018  -0.4192  10.4311

0.0000  0.0258  -0.2531  5.5725

0.0000  0.0332  0.2142  -3.3332



0.3000  -0.0713  -0.2934  7.3018

-0.0531  0.3126  -0.1995  4.4551

-0.0719  -0.0283  0.4083  -1.3305
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
#include<iostream>
#include<stdio.h>
#include<math.h>
#include <iomanip>
#include<vector>
using namespace std;
#define EPS 1.e-6
#define MaxSize 11

//简单地求一下Jacobi迭代法的系数AX=B,X=BX+F
void JacobiFactor(double a[][MaxSize], int N) {
	
	//定义一个N*N+1的容器来放置B,F
	vector<vector<double>>result(N);
	for (int i = 0; i < N; i++) {
		result[i].resize(N + 1);
	}

	//容器初始化
	for (int i = 0; i < N; i++) {
		for (int j = 0; j <= N; j++) {
			result[i][j] =a[i][j];
		}
	}

	//进行Jacobi的系数处理AX=B 变成X=BX+F
	for (int i = 0; i < N; i++) {
		for (int j = 0; j <= N; j++) {
			if (i == j) {
				result[i][j] = 0.0;
			}
			else if (j == N) {
				result[i][j] = a[i][j] / a[i][i];
			}
			else {
				result[i][j] = (-1.0) * a[i][j] / a[i][i];
			}
		}
	}

	for (int i = 0; i < N; i++) {
		for (int j = 0; j <= N; j++) {
			cout << setprecision(4) << result[i][j];
			if (j != N)cout << " ";
			if (j == N)cout << endl;
		}
	}

}




//求Gauss-Seidel迭代法的系数
void GaussSeidelFactor(double a[][MaxSize], int N) {
	//定义一个N*N+1的容器来放置系数
	vector<vector<double>>result(N);
	for (int i = 0; i < N; i++) {
		result[i].resize(N + 1);
	}

	//容器初始化
	for (int i = 0; i < N; i++) {
		for (int j = 0; j <= N; j++) {
			result[i][j] =0;
		}
	}

	//进行Gauss-Seidel的系数处理
	for (int i = 0; i < N; i++) {
		for (int j = 0; j <= N; j++) {
			if (i == 0) {
				if (i == j) {
					result[i][j] = 0.0;
				}
				else if(j==N){
					result[i][j] = a[i][j] / a[i][i];
				}
				else {
					result[i][j] = (-1.0) * a[i][j] / a[i][i];
				}
			}
			else {
				if (j < i) {
					for (int t = 0; t <= N; t++) {
						result[i][t] = result[i][t] - a[i][j] * result[j][t] / a[i][i];
					}
				}
				if (j > i &&j!=N) {
					result[i][j]= result[i][j]-a[i][j] / a[i][i];
				}
				if (j == N) {
					result[i][j] = result[i][j] + a[i][j] / a[i][i];
				}
			}
		}
	}

	for (int i = 0; i < N; i++) {
		for (int j = 0; j <= N; j++) {
			cout<< setprecision(4) << result[i][j];
			if (j != N)cout << " ";
			if (j == N)cout << endl;
		}
	}
}


//求松弛迭代法的系数
void SORFactor(double a[][MaxSize], int N,double W) {
	//定义一个N*N+1的容器来放置系数
	vector<vector<double>>result(N);
	for (int i = 0; i < N; i++) {
		result[i].resize(N + 1);
	}

	//容器初始化
	for (int i = 0; i < N; i++) {
		for (int j = 0; j <= N; j++) {
			result[i][j] = 0;
		}
	}

	//进行松弛法的系数处理
	for (int i = 0; i < N; i++) {
		for (int j = 0; j <= N; j++) {
			if (i == 0) {
				if (i == j) {
					result[i][j] = 1-W;
				}
				else if (j == N) {
					result[i][j] = W *a[i][j] / a[i][i];
				}
				else {
					result[i][j] = (-1.0) * W * a[i][j] / a[i][i];
				}
			}
			else {
				if (j < i) {
					for (int t = 0; t <= N; t++) {
						result[i][t] = result[i][t] - W*a[i][j] * result[j][t] / a[i][i];
					}
				}
				if (j > i && j != N) {
					result[i][j] = result[i][j]- W * a[i][j] / a[i][i];
				}
				if (j == N) {
					result[i][j] = result[i][j] +W * a[i][j] / a[i][i];
				}
				if (j == i) {
					result[i][j] = result[i][j] + (1 - W);
				}
			}
		}
	}

	for (int i = 0; i < N; i++) {
		for (int j = 0; j <= N; j++) {
			cout<< setprecision(4) << result[i][j];
			if (j != N)cout << " ";
			if (j == N)cout << endl;
		}
	}
}


int main() {
	int N = 0;
	double a[MaxSize][MaxSize] = { 0 };
	cin >> N;
	for (int i = 0; i < N; i++) {
		for (int j = 0; j <= N; j++) {
			cin >> a[i][j];
		}
	}
	cout << setiosflags(ios::fixed);
	JacobiFactor(a, N);
	cout << endl;
	GaussSeidelFactor(a, N);
	cout << endl;
	SORFactor(a, N, 0.7);
}


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
Gauss-Seidel迭代法求解线性方程组：

用Gauss-Seidel迭代法求解线性方程组Ax=b。输出迭代公式x=Bx+f的矩阵B和向量f以及10次迭代的结果。



输入格式：

第一行一个整数n，矩阵的行数(n<20)

接着n行，每行是矩阵的一行数据，按行顺序。每行n+1个实数，前n个实数是A的系数，最后一个是b的对应行的数。



输出格式：

前N行是矩阵B和向量f。

第i行N+1个实数，前N个实数是矩阵B的第i行的数据，最后一个数据是f的第i个数据。

第N+1行，是迭代10次后的解向量。 所有数精确到小数点后3位，数与数之间有空格，结束的位置没有空格



输入样例：

4
6 7 3 0 84
2 8 5 6 99
7 1 7 3 97
3 4 5 8 82



输出样例：

0.000 -1.167 -0.500 0.000 14.000
0.000 0.292 -0.500 -0.750 8.875
0.000 1.125 0.571 -0.321 -1.411
0.000 -0.411 0.080 0.576 1.444
2.885 5.203 10.058 0.273
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
#include<iostream>
#include<stdio.h>
#include<math.h>
#include <iomanip>
#include<vector>
using namespace std;
#define EPS 1.e-6
#define MaxSize 100


//Gauss-Seidel迭代法
void GaussSeidel(double a[][MaxSize], int N,int n) {
	//定义一个N*N+1的容器来放置系数
	vector<vector<double>>factors(N);
	for (int i = 0; i < N; i++) {
		factors[i].resize(N + 1);
	}

	//容器初始化
	for (int i = 0; i < N; i++) {
		for (int j = 0; j <= N; j++) {
			factors[i][j] =0;
		}
	}

	//进行Gauss-Seidel的系数处理
	for (int i = 0; i < N; i++) {
		for (int j = 0; j <= N; j++) {
			if (i == 0) {
				if (i == j) {
					factors[i][j] = 0.0;
				}
				else if(j==N){
					factors[i][j] = a[i][j] / a[i][i];
				}
				else {
					factors[i][j] = (-1.0) * a[i][j] / a[i][i];
				}
			}
			else {
				if (j < i) {
					for (int t = 0; t <= N; t++) {
						factors[i][t] = factors[i][t] - a[i][j] * factors[j][t] / a[i][i];
					}
				}
				if (j > i &&j!=N) {
					factors[i][j]= factors[i][j]-a[i][j] / a[i][i];
				}
				if (j == N) {
					factors[i][j] = factors[i][j] + a[i][j] / a[i][i];
				}
			}
		}
	}

	for (int i = 0; i < N; i++) {
		for (int j = 0; j <= N; j++) {
			if (abs(factors[i][j]) < EPS) {
				cout << setprecision(3) << 0.0;
			}
			else {
				cout << setprecision(3) << factors[i][j];
			}
			if (j != N)cout << " ";
			if (j == N)cout << endl;
		}
	}
	
	double* resultB = (double*)malloc(N * sizeof(double));
	if (resultB) {
		for (int i = 0; i < N; i++) {
			resultB[i] = 0.0;
		}
		for (int iter = 0; iter < 10; iter++) {

			for (int i = 0; i < N; i++) {
				double RESID = factors[i][N];
				for (int j = 0; j < N; j++) {
					RESID = RESID + factors[i][j] * resultB[j];
				}
				resultB[i] = RESID;
			}
		}
		for (int i = 0; i < N; i++) {
			cout << setprecision(3) << resultB[i];
			if (i != N - 1)cout << " ";
		}
		cout << endl;
	}

}




int main() {
	int N = 0;
	int n = 10;
	double a[MaxSize][MaxSize] = { 0 };
	cin >> N;
	for (int i = 0; i < N; i++) {
		for (int j = 0; j <= N; j++) {
			cin >> a[i][j];
		}
	}
	cout << setiosflags(ios::fixed);
	GaussSeidel(a, N,n);
}

二分法结果的真伪：针对一个非线性方程求根，输入二分法初始区间[a, b]。计算后，输出求得的结果及二分区间的次数，请判断这样的结果是否合理，合理意味存在某种非线性方程会产生这样的结果，不合理意味程序计算的结果不对或者就是伪造的结果。误差在1e-10范围内就是合理。 二分法计算过程如下：img 输入格式：第一行一个整数n，测试例子数接着n行，每行包括a, b, x, nit，[a, b]是初始区间，x是求得的解，nit是二分的次数。 输出格式：只有一行，每个例子输出“Yes”或“No”，合理是“Yes”，否则是“No”，结果之间有一个空格。 输入样例： 52.930 18.504 13.443210449218 122.601 15.841 9.707676611342 22.137 23.513 16.405287905720 182.579 21.821 18.547057947643 72.437 22.663 18.417603153228 18 输出样例：Yes No No No Yes

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
#include<iostream>
#include<stdio.h>
#include<vector>
#include <iomanip>
using namespace std;


int judge(long double a, long double b,long double goal) {

	long double num = goal;
	long double left = a;
	long double right = b;
	int counts = 0;
	long double mid = (a+b)/2.0;
	do {
		cout << "第" <<  counts << "次中分" << endl;
		cout << setiosflags(ios::fixed);
		cout << "left=" << setprecision(10) << left << endl;
		cout << "rigth=" << setprecision(10) << right << endl;
		cout << "goals=" << setprecision(10) << num << endl;
		cout << "mid=" << setprecision(10) << mid << endl;

		if (abs(mid-num)<=1.0e-10) {
			break;
		}
		if (num > mid) {
			cout << "num>mid" << endl;
			left = mid;
			right = right;
			mid = (left + right) / 2.0;
			counts += 1;
		}
		else {
			cout << "num<mid" << endl;
			right = mid;
			left = left;
			mid = (left + right) / 2.0;
			counts += 1;

		}

	} while (true);

	return counts;
}

int main() {
	int N = 0;
	cin >> N;
	//定义一个N*N+1的容器来放置系数
	vector<vector<long double>>data(N);
	for (int i = 0; i < N; i++) {
		data[i].resize(4);
	}

	//容器初始化
	for (int i = 0; i < N; i++) {
		for (int j = 0; j <= 3; j++) {
			cin >> data[i][j];
		}
	}

	//进行迭代次数的比对
	for (int i = 0; i < N; i++) {
		int flag = 0;
		flag = judge(data[i][0], data[i][1], data[i][2]);
		cout << setiosflags(ios::fixed);
		cout << flag << " ";
		if (abs(flag - data[i][3]) <= 1.0e-10) {
			cout << "Yes";
		}
		else {
			cout << "No";
		}
		if (i != N - 1)cout << " ";
		cout << endl;
	}
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
迭代加速法(Aitken法)求非线性方程的解

给定非线性方程exp(x)-a*x*x+b*x-c=0，其中a, b, c>0。请用牛顿法和Aitken法求其中一个根，初始值取为0，求得的前后两个近似解相差的绝对值小于1.e-8退出迭代，或者迭代次数达到50就退出迭代。Aitken的初始迭代公式就是牛顿迭代公式。



输入格式：

多行，每行一个测试用例，三个实数，a, b, c。直至“0 0 0”



输出格式：

多行，每行一个测试用例的结果。牛顿迭代法的迭代次数及根的结果，Aitken迭代法的迭代次数及根的结果，根精确到小数点后4位。数据之间有一个空格。



输入样例：

114 26 64

108 52 63

0 0 0



输出样例：

44 9.1440 22 9.1440

22 9.0377 11 9.0377
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
#include<iostream>
#include<vector>
#include<cmath>
#include <iomanip>
using namespace std;
#define eps 1.0e-8


//求函数的值
template<class _Ty>
 _Ty FunctionValue(_Ty x, vector<_Ty> factor) {
	_Ty result = exp(x)-factor[0] * x * x + factor[1] * x - factor[2];
	return result;
}

//求导数值
 template<class _Ty>
 _Ty FunctionCValue(_Ty x, vector<_Ty> factor) {
	 _Ty result = exp(x) - 2 * factor[0] * x + factor[1];
	 return result;
 }


 //牛顿法
 void RootNewton(double x, vector<double> factor, int js) {
	 int flag = 0;
	 int r = js;
	 double result = x;
	 double y = FunctionValue(x, factor);
	 double y0 = FunctionCValue(x, factor);
	 result = result - y / y0;
	 while (abs(y / y0) >= eps and flag < 50) {
		 y = FunctionValue(result, factor);
		 y0 = FunctionCValue(result, factor);
		 result = result - y / y0;
		 flag += 1;
	 }
	 cout << setprecision(4) << flag << " " << result << " ";
 }

 //牛顿法??
 void RootNewton2(double x, vector<double> factor, int js) {
	 int flag = 0;
	 int r = js;
	 double result = x;
	 double y = FunctionValue(x, factor);
	 double y0 = FunctionCValue(x, factor);
	 result = result - y / y0;
	 while (abs(y / y0) >= eps and flag < 50) {
		 y = FunctionValue(result, factor);
		 y0 = FunctionCValue(result, factor);
		 result = result - y / y0;
		 y = FunctionValue(result, factor);
		 y0 = FunctionCValue(result, factor);
		 result = result - y / y0;
		 flag += 1;
	 }
	 cout << setprecision(4) << flag << " " << result ;
 }


//埃特金法求
 void RootAitken(double x, vector<double> factor, int js) {
	 int r = 0;
	 double k = x;
	 int flag = 0;
	 double x0 = x;
	 while ((flag == 0) && (r != js)) {
		 r += 1;
		 double u = FunctionValue(x0, factor);
		 double v = FunctionValue(u, factor);
		 if (abs(u - v) < eps) {
			 x0 = v;
			 flag = 1;
		 }
		 else {
			 x0 = v - (v - u) * (v - u) / (v - 2.0 * u + x0);
		 }
		 k = x0;
		 r = js - r;
	 }
	 cout << setprecision(4) << r << " " << k;
 }

int main() {
	vector<double>v;
	double a = 1.0;
	vector<double >factors(3);
	int flag = 0;
	while (cin>>a) {
		v.push_back(a);
		if (abs(a) <= 1.0e-6)flag += 1;
		else flag = 0;
		if (flag == 3)break;
	}
	v.pop_back();
	v.pop_back();
	v.pop_back();
	for (int i = 0; i < v.size() / 3; i++) {
		for (int j = 0; j < 3; j++) {
			factors[j] = v[3 * i + j];
		}
		cout << setiosflags(ios::fixed);
		RootNewton(0.0,factors,50);
		RootNewton2(0.0, factors, 50);
		cout << endl;
	}
	
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
弦位法求解高阶多项式：

给定n阶多项式Pn(x)，先用弦位法求解一个根，两个初始值取为0和1，求得的前后两个近似解相差的绝对值小于1.e-8退出迭代，或者迭代次数达到50就退出迭代。求得一个根之后，比如5，用（x-5）除Pn(x)，得到一个新的n-1阶多项式，然后重复用这个过程，求出所有根。



输入格式：

第一行，一个整数n，多项式的阶（n<8）

第二行，（n+1）个实数，是多项式的系数，顺序为常数项、x的系数、x平方的系数等等，最后一个系数保证为1。数与数之间有一个空格。数据确保根为实数而且互异。



输出格式：

n个实数，依次求得的根，精确到小数点后4位，数与数之间有一个空格。

注意：用其它方法求得的根，可能顺序不同，而且有可能最后一两位有效数值不同。一定要用题目给出的方法求。这种方法也是通用的求高阶多项式的方法。



输入样例：

3

-5.7536 13.7022 -6.9327 1.0000



输出样例：

0.5716 2.9578 3.4033
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
#include<iostream>
#include<cmath>
#include<iomanip>
#include<vector>
using namespace std;
#define eps 1.0e-8

template<class _Ty>
inline _Ty FunctionValue(vector<_Ty>factors,vector<_Ty>results, _Ty x) {
	int n = factors.size();
	_Ty result = 0;
	for (int i = n - 1; i >= 0; i--) {
		result += pow(x, i) * factors[i];
	}
	for (int i = 0; i < results.size(); i++) {
		result = result / (x - results[i]);
	}
	return result;
}

template<class _Ty>
inline void SecantMethod(vector<_Ty>&factors) {

	int N = factors.size();
	vector<_Ty>result;

	for (int i = 0; i < N-1; i++) {
		_Ty x0 = 0.0;
		_Ty x1 = 1.0;
		_Ty temp = 0.0;
		size_t flag = 0;
		do {
			temp = FunctionValue(factors,result, x1) * (x1 - x0) /
				(FunctionValue(factors,result, x1) - FunctionValue(factors,result, x0));
			x0 = x1;
			x1 = x1 - temp;
			flag += 1;
		} while (abs(temp) > eps && flag <= 50);
		result.push_back(x1);
	}

	//vector 越界sad
	cout << setiosflags(ios::fixed);
	cout << endl;
	for (int i = 0; i < result.size(); i++) {
		cout <<setprecision(4)<< result[i]<<" ";
	}

}

int main() {

	int N = 0;
	cin >> N;
	vector<double>factors(N + 1);
	for (int i = 0; i <= N; i++) {
		cin >> factors[i];
	}
	SecantMethod(factors);
	
}
1
2
3