diff options
Diffstat (limited to 'experimental/Intersection/QuarticRoot.cpp')
| -rw-r--r-- | experimental/Intersection/QuarticRoot.cpp | 24 |
1 files changed, 22 insertions, 2 deletions
diff --git a/experimental/Intersection/QuarticRoot.cpp b/experimental/Intersection/QuarticRoot.cpp index 839c3b973b..b2e73cba50 100644 --- a/experimental/Intersection/QuarticRoot.cpp +++ b/experimental/Intersection/QuarticRoot.cpp @@ -120,7 +120,7 @@ static int cubicRootsX(double A, double B, double C, double D, double s[3]) { if (approximately_zero(A)) { // we're just a quadratic return quadraticRootsX(B, C, D, s); } - if (approximately_zero(D)) { + if (approximately_zero(D)) { // 0 is one root int num = quadraticRootsX(A, B, C, s); for (int i = 0; i < num; ++i) { if (approximately_zero(s[i])) { @@ -130,6 +130,16 @@ static int cubicRootsX(double A, double B, double C, double D, double s[3]) { s[num++] = 0; return num; } + if (approximately_zero(A + B + C + D)) { // 1 is one root + int num = quadraticRootsX(A, A + B, -D, s); + for (int i = 0; i < num; ++i) { + if (approximately_equal(s[i], 1)) { + return num; + } + } + s[num++] = 1; + return num; + } double a, b, c; { double invA = 1 / A; @@ -197,7 +207,7 @@ int quarticRoots(const double A, const double B, const double C, const double D, } int num; int i; - if (approximately_zero(E)) { + if (approximately_zero(E)) { // 0 is one root num = cubicRootsX(A, B, C, D, s); for (i = 0; i < num; ++i) { if (approximately_zero(s[i])) { @@ -207,6 +217,16 @@ int quarticRoots(const double A, const double B, const double C, const double D, s[num++] = 0; return num; } + if (approximately_zero(A + B + C + D + E)) { // 1 is one root + num = cubicRootsX(A, A + B, -(D + E), -E, s); // note that -C==A+B+D+E + for (i = 0; i < num; ++i) { + if (approximately_equal(s[i], 1)) { + return num; + } + } + s[num++] = 1; + return num; + } double u, v; /* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */ const double invA = 1 / A; |