htswanalysis/MACS/lib/gsl/gsl-1.11/poly/zsolve_cubic.c

   1 /* poly/zsolve_cubic.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* zsolve_cubic.c - finds the complex roots of x^3 + a x^2 + b x + c = 0 */
  21
  22 #include <config.h>
  23 #include <math.h>
  24 #include <gsl/gsl_math.h>
  25 #include <gsl/gsl_complex.h>
  26 #include <gsl/gsl_poly.h>
  27
  28 #define SWAP(a,b) do { double tmp = b ; b = a ; a = tmp ; } while(0)
  29
  30 int
  31 gsl_poly_complex_solve_cubic (double a, double b, double c,
  32                               gsl_complex *z0, gsl_complex *z1,
  33                               gsl_complex *z2)
  34 {
  35   double q = (a * a - 3 * b);
  36   double r = (2 * a * a * a - 9 * a * b + 27 * c);
  37
  38   double Q = q / 9;
  39   double R = r / 54;
  40
  41   double Q3 = Q * Q * Q;
  42   double R2 = R * R;
  43
  44   double CR2 = 729 * r * r;
  45   double CQ3 = 2916 * q * q * q;
  46
  47   if (R == 0 && Q == 0)
  48     {
  49       GSL_REAL (*z0) = -a / 3;
  50       GSL_IMAG (*z0) = 0;
  51       GSL_REAL (*z1) = -a / 3;
  52       GSL_IMAG (*z1) = 0;
  53       GSL_REAL (*z2) = -a / 3;
  54       GSL_IMAG (*z2) = 0;
  55       return 3;
  56     }
  57   else if (CR2 == CQ3)
  58     {
  59       /* this test is actually R2 == Q3, written in a form suitable
  60          for exact computation with integers */
  61
  62       /* Due to finite precision some double roots may be missed, and
  63          will be considered to be a pair of complex roots z = x +/-
  64          epsilon i close to the real axis. */
  65
  66       double sqrtQ = sqrt (Q);
  67
  68       if (R > 0)
  69         {
  70           GSL_REAL (*z0) = -2 * sqrtQ - a / 3;
  71           GSL_IMAG (*z0) = 0;
  72           GSL_REAL (*z1) = sqrtQ - a / 3;
  73           GSL_IMAG (*z1) = 0;
  74           GSL_REAL (*z2) = sqrtQ - a / 3;
  75           GSL_IMAG (*z2) = 0;
  76         }
  77       else
  78         {
  79           GSL_REAL (*z0) = -sqrtQ - a / 3;
  80           GSL_IMAG (*z0) = 0;
  81           GSL_REAL (*z1) = -sqrtQ - a / 3;
  82           GSL_IMAG (*z1) = 0;
  83           GSL_REAL (*z2) = 2 * sqrtQ - a / 3;
  84           GSL_IMAG (*z2) = 0;
  85         }
  86       return 3;
  87     }
  88   else if (CR2 < CQ3)  /* equivalent to R2 < Q3 */
  89     {
  90       double sqrtQ = sqrt (Q);
  91       double sqrtQ3 = sqrtQ * sqrtQ * sqrtQ;
  92       double theta = acos (R / sqrtQ3);
  93       double norm = -2 * sqrtQ;
  94       double r0 = norm * cos (theta / 3) - a / 3;
  95       double r1 = norm * cos ((theta + 2.0 * M_PI) / 3) - a / 3;
  96       double r2 = norm * cos ((theta - 2.0 * M_PI) / 3) - a / 3;
  97
  98       /* Sort r0, r1, r2 into increasing order */
  99
 100       if (r0 > r1)
 101         SWAP (r0, r1);
 102
 103       if (r1 > r2)
 104         {
 105           SWAP (r1, r2);
 106
 107           if (r0 > r1)
 108             SWAP (r0, r1);
 109         }
 110
 111       GSL_REAL (*z0) = r0;
 112       GSL_IMAG (*z0) = 0;
 113
 114       GSL_REAL (*z1) = r1;
 115       GSL_IMAG (*z1) = 0;
 116
 117       GSL_REAL (*z2) = r2;
 118       GSL_IMAG (*z2) = 0;
 119
 120       return 3;
 121     }
 122   else
 123     {
 124       double sgnR = (R >= 0 ? 1 : -1);
 125       double A = -sgnR * pow (fabs (R) + sqrt (R2 - Q3), 1.0 / 3.0);
 126       double B = Q / A;
 127
 128       if (A + B < 0)
 129         {
 130           GSL_REAL (*z0) = A + B - a / 3;
 131           GSL_IMAG (*z0) = 0;
 132
 133           GSL_REAL (*z1) = -0.5 * (A + B) - a / 3;
 134           GSL_IMAG (*z1) = -(sqrt (3.0) / 2.0) * fabs(A - B);
 135
 136           GSL_REAL (*z2) = -0.5 * (A + B) - a / 3;
 137           GSL_IMAG (*z2) = (sqrt (3.0) / 2.0) * fabs(A - B);
 138         }
 139       else
 140         {
 141           GSL_REAL (*z0) = -0.5 * (A + B) - a / 3;
 142           GSL_IMAG (*z0) = -(sqrt (3.0) / 2.0) * fabs(A - B);
 143
 144           GSL_REAL (*z1) = -0.5 * (A + B) - a / 3;
 145           GSL_IMAG (*z1) = (sqrt (3.0) / 2.0) * fabs(A - B);
 146
 147           GSL_REAL (*z2) = A + B - a / 3;
 148           GSL_IMAG (*z2) = 0;
 149         }
 150
 151       return 3;
 152     }
 153 }