nerdcruft/libcruft-util
1
0
Fork 0
Code Issues Releases Activity

polynomial: set NaNs as we go

Browse Source
This commit is contained in:
Danny Robson 2015-01-29 15:38:07 +11:00
parent e98bd6338e
commit 57d5538717
1 changed files with 8 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
polynomial.cpp
View File

@ -80,11 +80,7 @@ namespace util { namespace polynomial {
return {s[0], s[1], std::numeric_limits<float>::quiet_NaN () }; return {s[0], s[1], std::numeric_limits<float>::quiet_NaN () };
} }
std::array<float,3> s = { std::array<float,3> s;
std::numeric_limits<float>::quiet_NaN (),
std::numeric_limits<float>::quiet_NaN (),
std::numeric_limits<float>::quiet_NaN ()
};
// Normalise to x^3 + ax^2 + bx + c = 0 // Normalise to x^3 + ax^2 + bx + c = 0
const float a = _b / _a; const float a = _b / _a;
@ -103,10 +99,13 @@ namespace util { namespace polynomial {
{ {
if (almost_zero (q)) { if (almost_zero (q)) {
s[0] = 0.f; s[0] = 0.f;
s[1] = std::numeric_limits<float>::quiet_NaN ();
s[2] = std::numeric_limits<float>::quiet_NaN ();
} else { } else {
const float u = std::cbrt (-q); const float u = std::cbrt (-q);
s[0] = 2 * u; s[0] = 2 * u;
s[1] = -u; s[1] = -u;
s[2] = std::numeric_limits<float>::quiet_NaN ();
} }
} else if (D < 0) { } else if (D < 0) {
const float phi = std::acos (-q / std::sqrt (-p * p * p)) / 3.f; const float phi = std::acos (-q / std::sqrt (-p * p * p)) / 3.f;
@ -116,11 +115,14 @@ namespace util { namespace polynomial {
s[1] = -t * std::cos (phi + PI_f / 3.f); s[1] = -t * std::cos (phi + PI_f / 3.f);
s[2] = -t * std::cos (phi - PI_f / 3.f); s[2] = -t * std::cos (phi - PI_f / 3.f);
} else { } else {
float u = std::cbrt (std::sqrt (D) + std::fabs (q)); float u = std::cbrt (std::sqrt (D) + abs (q));
if (q > 0.f) if (q > 0.f)
s[0] = -u + p / u; s[0] = -u + p / u;
else else
s[0] = u - p / u; s[0] = u - p / u;
s[1] = std::numeric_limits<float>::quiet_NaN ();
s[2] = std::numeric_limits<float>::quiet_NaN ();
} }
// Resubstitute a / 3 from above // Resubstitute a / 3 from above