wendt@bocklin.UUCP (12/10/85)
Hi Folks: Here's a version of Bennet Todd's cubic equation solver that works on everything I tried. Part of the problem was that pow doesn't work if the first argument is negative, even on simple cases (ie cube root of -1). The complex math has been much cut back and could even be completely eliminated. Some of the changes you may be less fond of; I'm using struct functions and struct args. The command line syntax is simply "cubic c3 c2 c1 c0", where ci is the coefficient on x^i. The source for the algebraic changes was Chrystal: Textbook of Algebra. This is a very old book that I got out of the library, probably any algebra book of a similar era would have a similar treatment. The older ones are better because they're less abstract. People actually solved cubics by hand in those days. I had fun working on cubic. Even Macsyma has trouble with the cubic I was trying to solve. It gives me back three roots but it won't express them numerically as a+bi (gets floating exeptions), and all three of them seem to have imaginary parts. So at least now we have a working (albeit numeric) cubic solver. Now, who's going to do quartic? 3 2 560 z - 1232 z - 475 z + 3 (d6) - ---------------------------- = 0 560 >From "An Analysis of Alpha-Beta Pruning" by Knuth & Moore. Alan Wendt U of A CS arizona!wendt ---- don't process this with csh. or sh. --- ---- cut here and don't run though anything --- #include <stdio.h> #include <math.h> typedef struct complex { double real, imag; } COMPLEX; COMPLEX Complex(), /* returns a complex made from the real args */ ctimes(), /* product -- dyadic */ cplus(); /* sum -- dyadic */ #define Pi 3.14159265358979323846 /* * cubic -- apply closed-form solution to an arbitrary cubic * in real coefficients * based on the exposition of the algorithm in the CRC Standard * Mathematical Tables, 27'th edition, CRC Press, Boca Raton FL, * 1984, on page 9, and Chrystal: Textbook of Algebra * * syntax: cubic c3 c2 c1 c0 * Where the eqn is c3 x^3 + c2 x^2 + c1 x + c0 = 0 * ci are arbitrary constants in format suitable * for conversion using "%f" format in scanf(). */ /* * Written by Bennett Todd * Date: 10/15/85 -- initial release * This code is released to the public domain; do with it as you wish. * If you use it, credit would be appreciated. * Please send any bug fixes, enhancements, or general comments to * Bennett Todd * Duke Computation Center * Durham, NC 27706-7756 * +1 919 684 3695 * UUCP: ...{decvax,seismo,philabs,ihnp4,akgua}!mcnc!ecsvax!duccpc!bet * BITNET: dbtodd@tucc * * Changed by: Alan Wendt, U of AZ CS, arizona!wendt. */ main(argc, argv) int argc; char **argv; { /* Sorry about the variable names; they are chosen to match CRC */ double o, p, q, r, a, b, t, l, m; COMPLEX root_1, root_2, root_3, omega, omegasquared; omega = Complex(-0.5, sqrt(3.0) / 2.0); /* third roots of unity */ omegasquared = Complex(-0.5, -sqrt(3.0) / 2.0); /* (attempt to) parse command line */ if (argc != 5) { fprintf(stderr, "syntax: cubic c3 c2 c1 c0\n"); fprintf(stderr, "Indicate missing terms with zeroes.\n"); exit(1); } if (sscanf(argv[1], "%lf", &o) != 1) { fprintf(stderr, "cubic: cannot parse %s\n", argv[1]); exit(1); } if (sscanf(argv[2], "%lf", &p) != 1) { fprintf(stderr, "cubic: cannot parse %s\n", argv[2]); exit(1); } if (sscanf(argv[3], "%lf", &q) != 1) { fprintf(stderr, "cubic: cannot parse %s\n", argv[3]); exit(1); } if (sscanf(argv[4], "%lf", &r) != 1) { fprintf(stderr, "cubic: cannot parse %s\n", argv[4]); exit(1); } /* sanity check */ if (o == 0) { fprintf(stderr, "cubic: sorry buddy, that's a quadratic.\n"); exit(1); } p /= o; q /= o; r /= o; /* Change of variable from "x^3 + px^2 + qx + r = 0" to * "x'^3 + ax' + b = 0" by propagating the substitution * "x = x' - p/3". */ a = q - p * p / 3.0; b = (2.0 * p * p * p - 9.0 * p * q + 27.0 * r) / 27.0; t = b * b / 4.0 + a * a * a / 27.0; if (t < 0.0) { /* as b^2 / 4 + a^3/27 < 0, a < 0, -a and rho > 0 */ double rho, theta; rho = pow(-a, 1.5) / sqrt(27.0); theta = acos(b / rho / 2.0); rho = -2.0 * pow(rho, 1.0 / 3.0); root_1 = Complex(cos(theta / 3.0) * rho, 0.0); root_2 = Complex(cos((2.0 * Pi + theta) / 3.0) * rho, 0.0); root_3 = Complex(cos((4.0 * Pi + theta) / 3.0) * rho, 0.0); } else { t = sqrt(t); l = b / 2 + t; m = b / 2 - t; /* pow doesn't work for negative bases */ if (l >= 0) l = - pow(l, 1.0 / 3.0); else l = pow(-l, 1.0 / 3.0); if (m >= 0) m = -pow(m, 1.0 / 3.0); else m = pow(-m, 1.0 / 3.0); root_1 = Complex(l + m, 0.0); root_2 = cplus( ctimes(omega, Complex(l, 0.0)), ctimes(omegasquared, Complex(m, 0.0)) ); root_3 = cplus( ctimes(omegasquared, Complex(l, 0.0)), ctimes(omega, Complex(m, 0.0)) ); } p /= 3.0; printf("%.16g", root_1.real - p); if (root_1.imag != 0) printf("%s%.16g i", root_1.imag < 0 ? " - " : " + ", fabs(root_1.imag)); printf("\n%.16g", root_2.real - p); if (root_2.imag != 0) printf("%s%.16g i", root_2.imag < 0 ? " - " : " + ", fabs(root_2.imag)); printf("\n%.16g", root_3.real - p); if (root_3.imag != 0) printf("%s%.16g i", root_3.imag < 0 ? " - " : " + ", fabs(root_3.imag)); printf("\n"); } /* complex math routines. */ COMPLEX Complex(x, y) double x, y; { COMPLEX result; result.real = x; result.imag = y; return result; } COMPLEX ctimes(a, b) COMPLEX a, b; { COMPLEX result; result.real = a.real*b.real - a.imag*b.imag; result.imag = a.real*b.imag + a.imag*b.real; return result; } COMPLEX cplus(a, b) COMPLEX a, b; { COMPLEX result; result.real = a.real + b.real; result.imag = a.imag + b.imag; return result; } /* I'll stop breaking the law when they repeal it. */