Skip to content
Snippets Groups Projects
Commit e5167aba authored by Oliver Sander's avatar Oliver Sander Committed by sander@PCPOOL.MI.FU-BERLIN.DE
Browse files

use stable formula to solve quadratic equation

[[Imported from SVN: r2416]]
parent 2b8a71c2
Branches
Tags
No related merge requests found
......@@ -16,17 +16,17 @@ class TruncatedCGSolver : public IterativeSolver<VectorType>
typedef typename VectorType::field_type field_type;
/** \brief Returns positive root of a quadratic equation
\todo Uses the standard formula which is known to be unstable
*/
static field_type positiveRoot(const field_type& a, const field_type& b, const field_type& c) {
field_type root1 = (b + std::sqrt(b*b - 4*a*c)) / (2*a);
field_type root2 = (b - std::sqrt(b*b - 4*a*c)) / (2*a);
field_type p = b/a;
field_type q = c/a;
std::cout << "roots: " << root1 << ", " << root2 << std::endl;
std::cout << "values: " << a*root1*root1 + b*root1 + c
<< ", " << a*root2*root2 + b*root2 + c << std::endl;
field_type root1 = - p/2 - sign(p) * std::sqrt(p*p/4 - q);
field_type root2 = q/root1;
// There must be one positive and one negative root
assert(root1*root2 <= 0);
return std::max(root1, root2);
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment