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

Linux Cross Reference
Tina4/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 ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.