polynomial: use newtons method after cubic solve

This commit is contained in:
Danny Robson 2015-01-22 00:27:46 +11:00
parent 13188ca83c
commit bd12519d94
3 changed files with 40 additions and 4 deletions

View File

@ -126,6 +126,44 @@ namespace util { namespace polynomial {
float sub = a / 3.f; float sub = a / 3.f;
for (auto &i: s) for (auto &i: s)
i -= sub; i -= sub;
// Run an iteration of Newtons method to make the results slightly
// more accurate, they're a little loose straight out of the bat.
float da = 3;
float db = 2 * a;
float dc = b;
for (auto &i: s) {
float deriv = da * i * i + db * i + dc;
if (almost_zero (deriv))
continue;
i = i - eval (coeffs, i) / deriv;
}
return s; return s;
} }
} } } }
//-----------------------------------------------------------------------------
template <size_t S>
float
util::polynomial::eval (const std::array<float,S> coeffs, const float x)
{
float x_ = 1.f;
float sum = 0.f;
for (size_t i = 0; i < S; ++i) {
sum += coeffs[S-i-1] * x_;
x_ *= x;
}
return sum;
}
//-----------------------------------------------------------------------------
template float util::polynomial::eval (std::array<float,1>, float);
template float util::polynomial::eval (std::array<float,2>, float);
template float util::polynomial::eval (std::array<float,3>, float);
template float util::polynomial::eval (std::array<float,4>, float);

View File

@ -31,7 +31,7 @@ namespace util {
template <size_t S> template <size_t S>
float float
eval (std::array<float,S>); eval (std::array<float,S>, float x);
} }
} }

View File

@ -33,9 +33,7 @@ main (int, char**)
std::sort (s.begin (), s.end ()); std::sort (s.begin (), s.end ());
for (size_t j = 0; j < 3; ++j) { for (size_t j = 0; j < 3; ++j) {
std::cerr << i.solutions[j] << "==" << s[j] << '\n'; bool equal = ieee_single::almost_equal (i.solutions[j], s[j]);
bool equal = ieee_single::almost_equal (i.solutions[j], s[j], 864'026'624);
bool invalid = std::isnan (i.solutions[j]) && std::isnan (s[j]); bool invalid = std::isnan (i.solutions[j]) && std::isnan (s[j]);
CHECK (equal || invalid); CHECK (equal || invalid);