~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

# Linux Cross ReferenceTina4/src/math/numeric/eqn.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

```  1 /**@(#)Quadratic & cubic roots
2 **/
3
4 #include <math.h>
5 #include <tina/sys.h>
6 #include <tina/sysfuncs.h>
7 #include <tina/math.h>
8
9 Bool    quadratic_roots(double a, double b, double c, double *x1, double *x2)
10 {
11     double  disc = b * b - 4.0 * a * c;
12
13     /** assumes non-degenerate quadratic **/
14
15     if (b < 0)
16     {
17         a *= -1.0;
18         b *= -1.0;
19         c *= -1.0;
20     }
21     if (disc >= 0)
22     {                           /**real roots **/
23         double  q = -0.5 * (b + sqrt(disc));
24
25         *x1 = q / a;
26         *x2 = c / q;
27         return (true);
28     } else
29     {
30         *x1 = -b / (2.0 * a);   /** real part **/
31         *x2 = sqrt(-disc) / (2.0 * a);  /** imag part **/
32         return (false);
33     }
34 }
35
36 Bool    cubic_roots(double a, double b, double c, double d, double *x1, double *x2, double *x3)
37 {
38     double  disc, q, r;
39
40     /** assumes non-degenerate cubic **/
41     b /= a;
42     c /= a;
43     d /= a;
44
45     q = (b * b - 3.0 * c) / 9.0;
46     r = (2.0 * b * b * b - 9.0 * b * c + 27.0 * d) / 54.0;
47     disc = q * q * q - r * r;
48
49     if (disc >= 0)
50     {
51         double  theta;
52         double  sq = sqrt(q), b3 = b / 3.0, sq3;
53
54         sq3 = sq * sq * sq;
55         if (sq3 == 0.0)
56             theta = 0.0;        /**arbitrary**/
57         else
58             theta = acos(r / sq3);
59
60         *x1 = -2.0 * sq * cos(theta / 3.0) - b3;
61         *x2 = -2.0 * sq * cos((theta + TWOPI) / 3.0) - b3;
62         *x3 = -2.0 * sq * cos((theta + 2.0 * TWOPI) / 3.0) - b3;
63         return (true);
64     } else
65     {
66         double  s = pow(fabs(r) + sqrt(-disc), 1.0 / 3.0);
67         double  b3 = b / 3.0;
68         double  x;
69
70         if (r >= 0.0)
71             x = -(s + q / s);
72         else
73             x = (s + q / s);
74
75         x -= b3;
76
77         b += x;
78         c += x * b;
79
80         *x1 = x;                /** real root **/
81         return (quadratic_roots(1.0, b, c, x2, x3));    /** complex roots **/
82     }
83 }
84 ```

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~