From 7e0f021a9aec35fd8e6725e87e3313b101d26f5e Mon Sep 17 00:00:00 2001 From: Tobias Klauser Date: Sun, 27 Jan 2008 11:37:44 +0100 Subject: Initial import (2.0.2-6) --- reference/C/CONTRIB/SNIP/cubic.c | 61 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) create mode 100755 reference/C/CONTRIB/SNIP/cubic.c (limited to 'reference/C/CONTRIB/SNIP/cubic.c') diff --git a/reference/C/CONTRIB/SNIP/cubic.c b/reference/C/CONTRIB/SNIP/cubic.c new file mode 100755 index 0000000..cced533 --- /dev/null +++ b/reference/C/CONTRIB/SNIP/cubic.c @@ -0,0 +1,61 @@ +/* +** CUBIC.C - Solve a cubic polynomial +** public domain by Ross Cottrell +*/ + +#include +#include + +void SolveCubic(double a, + double b, + double c, + double d, + int *solutions, + double *x) +{ + long double a1 = b/a, a2 = c/a, a3 = d/a; + long double Q = (a1*a1 - 3.0*a2)/9.0; + long double R = (2.0*a1*a1*a1 - 9.0*a1*a2 + 27.0*a3)/54.0; + double R2_Q3 = R*R - Q*Q*Q; + + double theta; + + if (R2_Q3 <= 0) + { + *solutions = 3; + theta = acos(R/sqrt(Q*Q*Q)); + x[0] = -2.0*sqrt(Q)*cos(theta/3.0) - a1/3.0; + x[1] = -2.0*sqrt(Q)*cos((theta+2.0*M_PI)/3.0) - a1/3.0; + x[2] = -2.0*sqrt(Q)*cos((theta+4.0*M_PI)/3.0) - a1/3.0; + } + else + { + *solutions = 1; + x[0] = pow(sqrt(R2_Q3)+fabs(R), 1/3.0); + x[0] += Q/x[0]; + x[0] *= (R < 0.0) ? 1 : -1; + x[0] -= a1/3.0; + } +} + +#ifdef TEST + +int main(void) +{ + double a1 = 1.0, b1 = -10.5, c1 = 32.0, d1 = -30.0; + double a2 = 1.0, b2 = -4.5, c2 = 17.0, d2 = -30.0; + double x[3]; + int solutions; + + SolveCubic(a1, b1, c1, d1, &solutions, x); + + /* should get 3 solutions: 2, 6 & 2.5 */ + + SolveCubic(a2, b2, c2, d2, &solutions, x); + + /* should get 1 solution: 2.5 */ + + return 0; +} + +#endif /* TEST */ -- cgit v1.2.3-54-g00ecf