From 9a0f56a891ec127ffd1d0a407170be7e0f3794d6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 29 Dec 2021 09:55:54 -0500 Subject: [PATCH] Bug fix in quartic solver when d2 is close to zero This bug was resulted in lost particles for torus models with R ~ r. Bug was reported to author of quartic solver and a fix was made in the quarticpp project: https://github.com/cridemichel/quarticpp/issues/1 --- src/external/quartic_solver.cpp | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/external/quartic_solver.cpp b/src/external/quartic_solver.cpp index cfefa02ea3..17410654fd 100644 --- a/src/external/quartic_solver.cpp +++ b/src/external/quartic_solver.cpp @@ -11,11 +11,6 @@ constexpr double CUBIC_RESCAL_FACT = 3.488062113727083e+102; constexpr double QUART_RESCAL_FACT = 7.156344627944542e+76; constexpr double MACHEPS = std::numeric_limits::epsilon(); -double max3(double a, double b, double c) -{ - return std::max(std::max(a, b), c); -} - double solve_cubic_analytic_depressed_handle_inf(double b, double c) { /* find analytically the dominant root of a depressed cubic x^3+b*x+c @@ -479,7 +474,7 @@ void quartic_solver(double coeff[5], std::complex roots[4]) * better) */ if (realcase[0] == -1 || (std::abs(d2) <= - MACHEPS * oqs::max3(std::abs(2. * b / 3.), std::abs(phi0), l1 * l1))) { + MACHEPS * (std::abs(2. * b / 3.) + std::abs(phi0) + l1 * l1))) { double d3 = d - l3 * l3; double err0 = 0.0; if (realcase[0] == 1)