1 #include <math.h>
2 /*************************************************
3 Function: solve_quadratic_equation
4 Description: 求一元二次方程(a*x^2 + b*x + c = 0)的所有实数根
5 Input: 方程的系数 p = {c, b, a}
6 Output: 方程的所有实数根x
7 Return: 实数根的个数
8 Author: 枫箫
9 Version: 1.0
10 Date: 2017.7.8
11 Others: 如有疑问、意见或建议, 请多多交流, 多多指教!邮箱:gxb31415926@163.com
12 *************************************************/
13 int solve_quadratic_equation(float p[], float x[])
14 {
15 #define EPS 1.0e-30
16 #define ZERO 1.0e-30
17 float a, b, c, delta, sqrtDelta;
18
19 a = p[2];
20 b = p[1];
21 c = p[0];
22
23 if (fabs(a - 0.0) < EPS)
24 {
25 if (fabs(b - 0.0) < EPS)
26 {
27 return 0;
28 }
29 else
30 {
31 x[0] = -c / b;
32 return 1;
33 }
34 }
35 else
36 {
37 delta = b * b - 4.0 * a * c;
38 if (delta > ZERO)
39 {
40 if (fabs(c - 0.0) < EPS) //若c = 0,由于计算误差,sqrt(b*b - 4*a*c）不等于b
41 {
42 x[0] = 0.0;
43 x[1] = -b / a;
44 }
45 else
46 {
47 sqrtDelta = sqrt(delta);
48 if (b > 0.0)
49 {
50 x[0] = (-2.0 * c) / (b + sqrtDelta); //避免两个很接近的数相减,导致精度丢失
51 x[1] = (-b - sqrtDelta) / (2.0 * a);
52 }
53 else
54 {
55 x[0] = (-b + sqrtDelta) / (2.0 * a);
56 x[1] = (-2.0 * c) / (b - sqrtDelta); //避免两个很接近的数相减,导致精度丢失
57 }
58 }
59 return 2;
60 }
61 else if (fabs(delta - 0.0) < EPS)
62 {
63 x[0] = x[1] = -b / (2.0 * a);
64 }
65 else
66 {
67 return 0;
68 }
69 }
70 #undef EPS
71 #undef ZERO
72 }
73
74
75 /*************************************************
76 Function: solve_cubic_equation
77 Description: 盛金公式求一元三次方程(a*x^3 + b*x^2 + c*x + d = 0)的所有实数根
78 A = b * b - 3.0 * a * c;
79 B = b * c - 9.0 * a * d;
80 C = c * c - 3.0 * b * d;
81 (1)当A = B = 0时，方程有一个三重实根
82 (2)当Δ = B^2－4AC>0时，方程有一个实根和一对共轭虚根
83 (3)当Δ = B^2－4AC = 0时，方程有三个实根，其中有一个两重根
84 (4)当Δ = B^2－4AC<0时，方程有三个不相等的实根
85 Input: 方程的系数 p = {d, c, b, a}
86 Output: 方程的所有实数根x
87 Return: 实数根的个数
88 Author: 枫箫
89 Version: 1.0
90 Date: 2017.7.8
91 Others: 如有疑问、意见或建议, 请多多交流, 多多指教!邮箱:gxb31415926@163.com
92 *************************************************/
93 int solve_cubic_equation(float p[], float x[])
94 {
95 #define EPS 1.0e-30
96 #define ZERO 1.0e-30
97 float a, b, c, d, A, B, C, delta;
98 float Y1, Y2, Z1, Z2, K, parm[3], roots[2], theta, T;
99
100 a = p[3] ;
101 b = p[2] ;
102 c = p[1] ;
103 d = p[0] ;
104
105 if (fabs(a - 0.0) < EPS)
106 {
107 parm[2] = b;
108 parm[1] = c;
109 parm[0] = d;
110
111 return solve_quadratic_equation(parm, x);
112 }
113 else
114 {
115 A = b * b - 3.0 * a * c;
116 B = b * c - 9.0 * a * d;
117 C = c * c - 3.0 * b * d;
118
119 delta = B * B - 4.0 * A * C;
120
121 if (fabs(A - 0.0) < EPS && fabs(B - 0.0) < EPS)
122 {
123 x[0] = x[1] = x[2] = -b / (3.0 * a);
124 return 3;
125 }
126
127 if (delta > ZERO)
128 {
129 parm[2] = 1.0;
130 parm[1] = B;
131 parm[0] = A * C;
132
133 solve_quadratic_equation(parm, roots);
134 Z1 = roots[0];
135 Z2 = roots[1];
136
137 Y1 = A * b + 3.0 * a * Z1;
138 Y2 = A * b + 3.0 * a * Z2;
139
140 if (Y1 < 0.0 && Y2 < 0.0) //pow函数的底数必须为非负数,必须分类讨论
141 {
142 x[0] = (-b + pow(-Y1, 1.0 / 3.0) + pow(-Y2, 1.0 / 3.0)) / (3.0*a);
143 }
144 else if (Y1 < 0.0 && Y2 > 0.0)
145 {
146 x[0] = (-b + pow(-Y1, 1.0 / 3.0) - pow(Y2, 1.0 / 3.0)) / (3.0*a);
147 }
148 else if (Y1 > 0.0 && Y2 < 0.0)
149 {
150 x[0] = (-b - pow(Y1, 1.0 / 3.0) + pow(-Y2, 1.0 / 3.0)) / (3.0*a);
151 }
152 else
153 {
154 x[0] = (-b - pow(Y1, 1.0 / 3.0) - pow(Y2, 1.0 / 3.0)) / (3.0*a);
155 }
156 return 1;
157 }
158 else if (fabs(delta - 0.0) < EPS)
159 {
160 if (fabs(A - 0.0) > EPS)
161 {
162 K = B / A;
163 x[0] = -b / a + K;
164 x[1] = x[2] = -0.5 * K;
165 return 3;
166 }
167 else
168 {
169 return 0;
170 }
171 }
172 else
173 {
174 if (A > 0.0)
175 {
176 T = (2.0 * A * b - 3.0 * a * B) / (2.0 * pow(A, 3.0 / 2.0));
177 if (T > 1.0) //由于计算误差,T的值可能略大于1(如1.0000001)
178 {
179 T = 1.0;
180 }
181 if (T < -1.0)
182 {
183 T = -1.0;
184 }
185 theta = acos(T);
186 x[0] = (-b - 2.0 * sqrt(A) * cos(theta / 3.0)) / (3.0 * a);
187 x[1] = (-b + sqrt(A) * (cos(theta / 3.0) + sqrt(3.0) * sin(theta / 3.0))) / (3.0 * a);
188 x[2] = (-b + sqrt(A) * (cos(theta / 3.0) - sqrt(3.0) * sin(theta / 3.0))) / (3.0 * a);
189 return 3;
190 }
191 else
192 {
193 return 0;
194 }
195 }
196 }
197 #undef EPS
198 #undef ZERO
199 }
200
201
202 /*************************************************
203 Function: solve_quartic_equation
204 Description: 费拉里法求一元四次方程(a*x^4 + b*x^3 + c*x^2 + d*x + e = 0)的所有实数根
205 Input: 方程的系数 p = {e, d, c, b, a}
206 Output: 方程的所有实数根x
207 Return: 实数根的个数
208 Author: 枫箫
209 Version: 1.0
210 Date: 2017.7.8
211 Others: 如有疑问、意见或建议, 请多多交流, 多多指教!邮箱:gxb31415926@163.com
212 *************************************************/
213 int solve_quartic_equation(float p[], float x[])
214 {
215 #define EPS 1.0e-30
216
217 float a, b, c, d, e;
218 float parm[4], roots[3];
219 float y, M, N;
220 float x1[2], x2[2];
221 int rootCount1, rootCount2, rootCount, i;
222 float MSquareTemp, MSquare, yTemp;
223
224 a = p[4];
225 b = p[3];
226 c = p[2];
227 d = p[1];
228 e = p[0];
229
230 if (fabs(a - 0.0) < EPS)
231 {
232 if (fabs(b - 0.0) < EPS)
233 {
234 parm[2] = c;
235 parm[1] = d;
236 parm[0] = e;
237 return solve_quadratic_equation(parm, x);
238 }
239 else
240 {
241 parm[3] = b;
242 parm[2] = c;
243 parm[1] = d;
244 parm[0] = e;
245 return solve_cubic_equation(parm, x);
246 }
247 }
248 else
249 {
250 b = b / a;
251 c = c / a;
252 d = d / a;
253 e = e / a;
254
255 parm[3] = 8.0;
256 parm[2] = -4.0 * c;
257 parm[1] = 2.0 * (b * d - 4.0 * e);
258 parm[0] = -e * (b * b - 4.0 * c) - d * d;
259
260 if (rootCount = solve_cubic_equation(parm, roots))
261 {
262 y = roots[0];
263 MSquare = 8.0 * y + b * b - 4.0 * c;
264 for (i = 1; i < rootCount; i++)
265 {
266 yTemp = roots[i];
267 MSquareTemp = 8.0 * yTemp + b * b - 4.0 * c;
268 if (MSquareTemp > MSquare)
269 {
270 MSquare = MSquareTemp;
271 y = yTemp;
272 }
273 }
274
275 if (MSquare > 0.0)
276 {
277 M = sqrt(MSquare);
278 N = b * y - d;
279 parm[2] = 2.0;
280 parm[1] = b + M;
281 parm[0] = 2.0 * (y + N / M);
282 rootCount1 = solve_quadratic_equation(parm, x1);
283
284 parm[2] = 2.0;
285 parm[1] = b - M;
286 parm[0] = 2.0 * (y - N / M);
287 rootCount2 = solve_quadratic_equation(parm, x2);
288
289 if (rootCount1 == 2)
290 {
291 x[0] = x1[0];
292 x[1] = x1[1];
293 x[2] = x2[0];
294 x[3] = x2[1];
295 }
296 else
297 {
298 x[0] = x2[0];
299 x[1] = x2[1];
300 x[2] = x1[0];
301 x[3] = x1[1];
302 }
303 return rootCount1 + rootCount2;
304 }
305 else
306 {
307 return 0;
308 }
309 }
310 else
311 {
312 return 0;
313 }
314 }
315 #undef EPS
316 }
317
318
319 void main()
320 {
321 float x[2], xx[3], xxx[4], p[5];
322 int rootCount;
323 float breakPointHere, x1, x2, a, b, c, d, e;
324
325 //(1)一元二次方程测试
326 //x^2 - 1000000.000001*x + 1 = 0
327 //0*x^2 - 10*x + 1 = 0
328 //0 * x ^ 2 - 10 * x + 1 = 0
329 //x^2 - 10000000*x + 0.01 = 0
330 //1.0e-20*x^2 - 2.0e-20*x + 1.0e-20 = 0
331
332 a = 1;
333 b = -1000000.000001;
334 c = 1;
335 p[0] = c;
336 p[1] = b;
337 p[2] = a;
338 rootCount = solve_quadratic_equation(p, x);
339 x1 = (-b + sqrt(b * b - 4.0 * a * c)) / (2.0 * a);
340 x2 = (-b - sqrt(b * b - 4.0 * a * c)) / (2.0 * a);
341 breakPointHere = 1.0;
342
343 //(2)一元三次方程测试
344 //(x-1)*(x^2+1)=0 (x^3 - x^2 + x - 1 = 0)
345 //(x-1)^3 = 0 (x^3 - 3*x^2 + 3*x - 1 = 0)
346 //(x-1)^2*(x-2)=0 (x^3 - 4*x^2 + 5*x - 2 = 0)
347 //(x-1)*(x-2)*(x-3) = 0 (x^3 - 6*x^2 + 11*x - 6 = 0)
348 //0*x^3 + x^2 - 2*x + 1 = 0
349 //0*x^3 + 0*x^2 - 2*x + 1 = 0
350 //0*x^3 + 0*x^2 + 0*x + 1 = 0
351
352 a = 1;
353 b = -6;
354 c = 11;
355 d = -6;
356 p[0] = d;
357 p[1] = c;
358 p[2] = b;
359 p[3] = a;
360 rootCount = solve_cubic_equation(p, xx);
361 breakPointHere = 1.0;
362
363 //(3)一元四次方程测试
364 //(x-1)*(x-2)*(x^2 + 1)=0 (x^4 - 3*x^3 + 3*x^2 - 3*x + 2 = 0)
365 //(x-1)^2*(x^2 + 1)=0 (x^4 - 2*x^3 + 2*x^2 - 2*x + 1 = 0)
366 //(x-1)*(x-2)*(x-3)*(x-4)=0 (x^4 - 10*x^3 + 35*x^2 - 50*x + 24 = 0)
367 //(x-1)^2*(x-2)^2 = 0 (x^4 - 6*x^3 + 13*x^2 - 12*x + 4 = 0)
368 //0*x^4 + x^3 - 3*x^2 + 3*x - 1 = 0
369 //0*x^4 + 0*x^3 + x^2 - 2*x + 1 = 0
370 //0*x^4 + 0*x^3 + 0*x^2 - 2*x + 1 = 0
371 a = 1;
372 b = -10;
373 c = 35;
374 d = -50;
375 e = 24;
376 p[0] = e;
377 p[1] = d;
378 p[2] = c;
379 p[3] = b;
380 p[4] = a;
381 rootCount = solve_quartic_equation(p, xxx);
382 breakPointHere = 1.0;
383 }
浙公网安备 33010602011771号