/* poly/solve_quartic.c * * Copyright (C) 2003 CERN and K.S. K\"{o}lbig * * Converted to C and implemented into the GSL Library - Sept. 2003 * by Andrew W. Steiner and Andy Buckley * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or (at * your option) any later version. * * This program is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ /* solve_quartic.c - finds the real roots of * x^4 + a x^3 + b x^2 + c x + d = 0 */ #include #include #include #include /* Work needed to remove the complex numbers from this * algorithm: when done, compilation will work without these * includes */ #include #include #define SWAPD(a,b) do { double tmp = b ; b = a ; a = tmp ; } while(0) int gsl_poly_solve_quartic (double a, double b, double c, double d, double *x0, double *x1, double *x2, double *x3) { gsl_complex i, zarr[4], w1, w2, w3; double r4 = 1.0 / 4.0; double q2 = 1.0 / 2.0, q4 = 1.0 / 4.0, q8 = 1.0 / 8.0; double q1 = 3.0 / 8.0, q3 = 3.0 / 16.0; double u[3], v[3], v1, v2, disc; double aa, pp, qq, rr, rc, sc, tc, q, h; int k1 = 0, k2 = 0, mt; gsl_complex tmp[3]; GSL_SET_COMPLEX (&i, 0.0, 1.0); GSL_SET_COMPLEX (&zarr[0], 0.0, 0.0); GSL_SET_COMPLEX (&zarr[1], 0.0, 0.0); GSL_SET_COMPLEX (&zarr[2], 0.0, 0.0); GSL_SET_COMPLEX (&zarr[3], 0.0, 0.0); GSL_SET_COMPLEX (&w1, 0.0, 0.0); GSL_SET_COMPLEX (&w2, 0.0, 0.0); GSL_SET_COMPLEX (&w3, 0.0, 0.0); /* Deal easily with the cases where the quartic is degenerate. The * ordering of solutions is done explicitly. */ if (b == 0 && c == 0) { if (d == 0) { if (a > 0) { *x0 = -a; *x1 = 0.0; *x2 = 0.0; *x3 = 0.0; } else { *x0 = 0.0; *x1 = 0.0; *x2 = 0.0; *x3 = -a; } return 4; } else if (a == 0) { if (d > 0) { return 0; } else { *x0 = -sqrt (sqrt (d)); *x1 = -(*x0); return 2; } } } /* For non-degenerate solutions, proceed by constructing and * solving the resolvent cubic */ aa = a * a; pp = b - q1 * aa; qq = c - q2 * a * (b - q4 * aa); rr = d - q4 * (a * c - q4 * aa * (b - q3 * aa)); rc = q2 * pp; sc = q4 * (q4 * pp * pp - rr); tc = -(q8 * qq * q8 * qq); /* This code solves the resolvent cubic in a convenient fashion * for this implementation of the quartic. If there are three real * roots, then they are placed directly into u[]. If two are * complex, then the real root is put into u[0] and the real * and imaginary part of the complex roots are placed into * u[1] and u[2], respectively. Additionally, this * calculates the discriminant of the cubic and puts it into the * variable disc. */ { double qcub = (rc * rc - 3 * sc); double rcub = (2 * rc * rc * rc - 9 * rc * sc + 27 * tc); double Q = qcub / 9; double R = rcub / 54; double Q3 = Q * Q * Q; double R2 = R * R; double CR2 = 729 * rcub * rcub; double CQ3 = 2916 * qcub * qcub * qcub; disc = (CR2 - CQ3) / 2125764.0; if (R == 0 && Q == 0) { u[0] = -rc / 3; u[1] = -rc / 3; u[2] = -rc / 3; } else if (CR2 == CQ3) { double sqrtQ = sqrt (Q); if (R > 0) { u[0] = -2 * sqrtQ - rc / 3; u[1] = sqrtQ - rc / 3; u[2] = sqrtQ - rc / 3; } else { u[0] = -sqrtQ - rc / 3; u[1] = -sqrtQ - rc / 3; u[2] = 2 * sqrtQ - rc / 3; } } else if (CR2 < CQ3) { double sqrtQ = sqrt (Q); double sqrtQ3 = sqrtQ * sqrtQ * sqrtQ; double theta = acos (R / sqrtQ3); if (R / sqrtQ3 >= 1.0) theta = 0.0; double norm = -2 * sqrtQ; u[0] = norm * cos (theta / 3) - rc / 3; u[1] = norm * cos ((theta + 2.0 * M_PI) / 3) - rc / 3; u[2] = norm * cos ((theta - 2.0 * M_PI) / 3) - rc / 3; } else { double sgnR = (R >= 0 ? 1 : -1); double A = -sgnR * pow (fabs (R) + sqrt (R2 - Q3), 1.0 / 3.0); double B = Q / A; u[0] = A + B - rc / 3; u[1] = -0.5 * (A + B) - rc / 3; u[2] = -(sqrt (3.0) / 2.0) * fabs (A - B); } } /* End of solution to resolvent cubic */ /* Combine the square roots of the roots of the cubic * resolvent appropriately. Also, calculate 'mt' which * designates the nature of the roots: * mt=1 : 4 real roots (disc == 0) * mt=2 : 0 real roots (disc < 0) * mt=3 : 2 real roots (disc > 0) */ if (disc == 0) { u[2] = u[1]; } if (disc <= 0) { mt = 2; v[0] = fabs (u[0]); v[1] = fabs (u[1]); v[2] = fabs (u[2]); v1 = GSL_MAX (GSL_MAX (v[0], v[1]), v[2]); if (v1 == v[0]) { k1 = 0; v2 = GSL_MAX (v[1], v[2]); } else if (v1 == v[1]) { k1 = 1; v2 = GSL_MAX (v[0], v[2]); } else { k1 = 2; v2 = GSL_MAX (v[0], v[1]); } if (v2 == v[0]) { k2 = 0; } else if (v2 == v[1]) { k2 = 1; } else { k2 = 2; } w1 = gsl_complex_sqrt_real (u[k1]); w2 = gsl_complex_sqrt_real (u[k2]); } else { mt = 3; GSL_SET_COMPLEX (&w1, u[1], u[2]); GSL_SET_COMPLEX (&w2, u[1], -u[2]); w1 = gsl_complex_sqrt (w1); w2 = gsl_complex_sqrt (w2); } /* Solve the quadratic in order to obtain the roots * to the quartic */ q = qq; if (gsl_complex_abs (gsl_complex_mul (w1, w2)) != 0.0) { w3 = gsl_complex_mul_real (gsl_complex_inverse (gsl_complex_mul (w1, w2)), -q / 8.0); } h = r4 * a; zarr[0] = gsl_complex_add_real (gsl_complex_add (gsl_complex_add (w1, w2), w3), -h); zarr[1] = gsl_complex_add_real (gsl_complex_add (gsl_complex_negative (gsl_complex_add (w1, w2)), w3), -h); zarr[2] = gsl_complex_add_real (gsl_complex_sub (gsl_complex_sub (w2, w1), w3), -h); zarr[3] = gsl_complex_add_real (gsl_complex_sub (gsl_complex_sub (w1, w2), w3), -h); /* Arrange the roots into the variables z0, z1, z2, z3 */ if (mt == 2) { /* No real roots. Return 0; */ return 0; } else if (mt == 3) { *x0 = GSL_REAL (zarr[0]); *x1 = GSL_REAL (zarr[1]); } /* * Sort the roots as usual: main sorting by ascending real part, secondary * sorting by ascending imaginary part */ if (mt == 1) { /* Roots are all real, sort them by the real part */ if (*x0 > *x1) SWAPD (*x0, *x1); if (*x0 > *x2) SWAPD (*x0, *x2); if (*x0 > *x3) SWAPD (*x0, *x3); if (*x1 > *x2) SWAPD (*x1, *x2); if (*x2 > *x3) { SWAPD (*x2, *x3); if (*x1 > *x2) SWAPD (*x1, *x2); } return 4; } else { /* 2 real roots */ if (*x0 > *x1) SWAPD (*x0, *x1); } return 2; }