summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--libdimension/polynomial.c48
1 files changed, 24 insertions, 24 deletions
diff --git a/libdimension/polynomial.c b/libdimension/polynomial.c
index 136b0b5..68970ee 100644
--- a/libdimension/polynomial.c
+++ b/libdimension/polynomial.c
@@ -320,11 +320,10 @@ static inline size_t
dmnsn_solve_linear(const double poly[2], double x[1])
{
x[0] = -poly[0];
- if (x[0] >= dmnsn_epsilon) {
+ if (x[0] >= dmnsn_epsilon)
return 1;
- } else {
+ else
return 0;
- }
}
/** Solve a normalized quadratic polynomial algebraically. */
@@ -336,13 +335,13 @@ dmnsn_solve_quadratic(const double poly[3], double x[2])
double s = sqrt(disc);
x[0] = (-poly[1] + s)/2.0;
x[1] = (-poly[1] - s)/2.0;
- if (x[1] >= dmnsn_epsilon) {
+
+ if (x[1] >= dmnsn_epsilon)
return 2;
- } else if (x[0] >= dmnsn_epsilon) {
+ else if (x[0] >= dmnsn_epsilon)
return 1;
- } else {
+ else
return 0;
- }
} else {
return 0;
}
@@ -370,20 +369,12 @@ dmnsn_solve_cubic(double poly[4], double x[3])
x[0] = -2.0*msqrtp3*cos(4.0*atan(1.0)/3.0 - theta) - bdiv3;
x[1] = -(x[0] + x[2] + poly[2]);
- if (x[2] >= dmnsn_epsilon) {
+ if (x[2] >= dmnsn_epsilon)
return 3;
- } else if (x[1] >= dmnsn_epsilon) {
+ else if (x[1] >= dmnsn_epsilon)
return 2;
- } else if (x[0] >= dmnsn_epsilon) {
- return 1;
- } else {
- return 0;
- }
- } else {
+ } else if (disc >= dmnsn_epsilon) {
/* One real root */
- if (disc < 0.0)
- disc = 0.0;
-
double cbrtdiscq = cbrt(sqrt(disc/108.0) + fabs(q)/2.0);
double abst = cbrtdiscq - p/(3.0*cbrtdiscq);
@@ -392,13 +383,22 @@ dmnsn_solve_cubic(double poly[4], double x[3])
} else {
x[0] = abst - bdiv3;
}
-
- if (x[0] >= dmnsn_epsilon) {
- return 1;
- } else {
- return 0;
- }
+ } else if (fabs(p) < dmnsn_epsilon) {
+ /* Equation is a perfect cube */
+ x[0] = -bdiv3;
+ } else {
+ /* Two real roots; one duplicate */
+ double t1 = 3.0*q/p, t2 = -t1/2.0;
+ x[0] = dmnsn_max(t1, t2);
+ x[1] = dmnsn_min(t1, t2);
+ if (x[1] >= dmnsn_epsilon)
+ return 2;
}
+
+ if (x[0] >= dmnsn_epsilon)
+ return 1;
+ else
+ return 0;
}
/* Uspensky's algorithm */