首页 > 代码库 > 一元二次、三次、四次方程的实数根(c语言)

一元二次、三次、四次方程的实数根(c语言)

  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 }

 

一元二次、三次、四次方程的实数根(c语言)