From 60d54e565f081f39be8fc77087f46eb27c55d9d0 Mon Sep 17 00:00:00 2001 From: Nick Date: Sun, 2 Jan 2022 17:58:09 -0800 Subject: [PATCH] Quartic roots. (#718) --- doc/roots/quartic_roots.qbk | 40 +++++ doc/roots/roots_overview.qbk | 1 + include/boost/math/tools/quartic_roots.hpp | 150 ++++++++++++++++++ .../performance/quartic_roots_performance.cpp | 41 +++++ test/Jamfile.v2 | 1 + test/engel_expansion_test.cpp | 1 + test/math_unit_test.hpp | 36 +++-- test/quartic_roots_test.cpp | 137 ++++++++++++++++ test/simple_continued_fraction_test.cpp | 1 + test/test_z_test.cpp | 3 +- 10 files changed, 400 insertions(+), 11 deletions(-) create mode 100644 doc/roots/quartic_roots.qbk create mode 100644 include/boost/math/tools/quartic_roots.hpp create mode 100644 reporting/performance/quartic_roots_performance.cpp create mode 100644 test/quartic_roots_test.cpp diff --git a/doc/roots/quartic_roots.qbk b/doc/roots/quartic_roots.qbk new file mode 100644 index 000000000..fe1665704 --- /dev/null +++ b/doc/roots/quartic_roots.qbk @@ -0,0 +1,40 @@ +[/ +Copyright (c) 2021 Nick Thompson +Use, modification and distribution are subject to the +Boost Software License, Version 1.0. (See accompanying file +LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +] + +[section:quartic_roots Roots of Quartic Polynomials] + +[heading Synopsis] + +``` +#include + +namespace boost::math::tools { + +// Solves ax⁴ + bx³ + cx² + dx + e = 0. +std::array quartic_roots(Real a, Real b, Real c, Real d, Real e); + +} +``` + +[heading Background] + +The `quartic_roots` function extracts all real roots of a quartic polynomial ax⁴+ bx³ + cx² + dx + e. +The result is a `std::array`, which has length four, irrespective of the number of real roots the polynomial possesses. +(This is to prevent the performance overhead of allocating a vector, which often exceeds the time to extract the roots.) +The roots are returned in nondecreasing order. If a root is complex, then it is placed at the back of the array and set to a nan. + +The algorithm uses the classical method of Ferrari, and follows [@https://github.com/erich666/GraphicsGems/blob/master/gems/Roots3And4.c Graphics Gems V], +with an additional Halley iterate for root polishing. +A typical use of a quartic real-root solver is to raytrace a torus. + +[heading Performance and Accuracy] + +On a consumer laptop, we observe extraction of the roots taking ~90ns. +The file `reporting/performance/quartic_roots_performance.cpp` allows determination of the speed on your system. + +[endsect] +[/section:quartic_roots] diff --git a/doc/roots/roots_overview.qbk b/doc/roots/roots_overview.qbk index 770432af5..30bd2e385 100644 --- a/doc/roots/roots_overview.qbk +++ b/doc/roots/roots_overview.qbk @@ -18,6 +18,7 @@ There are several fully-worked __root_finding_examples, including: [include roots_without_derivatives.qbk] [include roots.qbk] [include cubic_roots.qbk] +[include quartic_roots.qbk] [include root_finding_examples.qbk] [include minima.qbk] [include root_comparison.qbk] diff --git a/include/boost/math/tools/quartic_roots.hpp b/include/boost/math/tools/quartic_roots.hpp new file mode 100644 index 000000000..1394e7fae --- /dev/null +++ b/include/boost/math/tools/quartic_roots.hpp @@ -0,0 +1,150 @@ +// (C) Copyright Nick Thompson 2021. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +#ifndef BOOST_MATH_TOOLS_QUARTIC_ROOTS_HPP +#define BOOST_MATH_TOOLS_QUARTIC_ROOTS_HPP +#include +#include +#include + +namespace boost::math::tools { + +namespace detail { + +// Make sure the nans are always at the back of the array: +template +bool comparator(Real r1, Real r2) { + using std::isnan; + if (isnan(r2)) { return true; } + return r1 < r2; +} + +template +std::array polish_and_sort(Real a, Real b, Real c, Real d, Real e, std::array& roots) { + // Polish the roots with a Halley iterate. + using std::fma; + using std::abs; + for (auto &r : roots) { + Real df = fma(4*a, r, 3*b); + df = fma(df, r, 2*c); + df = fma(df, r, d); + Real d2f = fma(12*a, r, 6*b); + d2f = fma(d2f, r, 2*c); + Real f = fma(a, r, b); + f = fma(f,r,c); + f = fma(f,r,d); + f = fma(f,r,e); + Real denom = 2*df*df - f*d2f; + if (abs(denom) > (std::numeric_limits::min)()) + { + r -= 2*f*df/denom; + } + } + std::sort(roots.begin(), roots.end(), detail::comparator); + return roots; +} + +} +// Solves ax^4 + bx^3 + cx^2 + dx + e = 0. +// Only returns the real roots, as these are the only roots of interest in ray intersection problems. +// Follows Graphics Gems V: https://github.com/erich666/GraphicsGems/blob/master/gems/Roots3And4.c +template +std::array quartic_roots(Real a, Real b, Real c, Real d, Real e) { + using std::abs; + using std::sqrt; + auto nan = std::numeric_limits::quiet_NaN(); + std::array roots{nan, nan, nan, nan}; + if (abs(a) <= (std::numeric_limits::min)()) { + auto cbrts = cubic_roots(b, c, d, e); + roots[0] = cbrts[0]; + roots[1] = cbrts[1]; + roots[2] = cbrts[2]; + if (b == 0 && c == 0 && d == 0 && e == 0) { + roots[3] = 0; + } + return detail::polish_and_sort(a, b, c, d, e, roots); + } + if (abs(e) <= (std::numeric_limits::min)()) { + auto v = cubic_roots(a, b, c, d); + roots[0] = v[0]; + roots[1] = v[1]; + roots[2] = v[2]; + roots[3] = 0; + return detail::polish_and_sort(a, b, c, d, e, roots); + } + // Now solve x^4 + Ax^3 + Bx^2 + Cx + D = 0. + Real A = b/a; + Real B = c/a; + Real C = d/a; + Real D = e/a; + Real Asq = A*A; + // Let x = y - A/4: + // Mathematica: Expand[(y - A/4)^4 + A*(y - A/4)^3 + B*(y - A/4)^2 + C*(y - A/4) + D] + // We now solve the depressed quartic y^4 + py^2 + qy + r = 0. + Real p = B - 3*Asq/8; + Real q = C - A*B/2 + Asq*A/8; + Real r = D - A*C/4 + Asq*B/16 - 3*Asq*Asq/256; + if (abs(r) <= (std::numeric_limits::min)()) { + auto [r1, r2, r3] = cubic_roots(Real(1), Real(0), p, q); + r1 -= A/4; + r2 -= A/4; + r3 -= A/4; + roots[0] = r1; + roots[1] = r2; + roots[2] = r3; + roots[3] = -A/4; + return detail::polish_and_sort(a, b, c, d, e, roots); + } + // Biquadratic case: + if (abs(q) <= (std::numeric_limits::min)()) { + auto [r1, r2] = quadratic_roots(Real(1), p, r); + if (r1 >= 0) { + Real rtr = sqrt(r1); + roots[0] = rtr - A/4; + roots[1] = -rtr - A/4; + } + if (r2 >= 0) { + Real rtr = sqrt(r2); + roots[2] = rtr - A/4; + roots[3] = -rtr - A/4; + } + return detail::polish_and_sort(a, b, c, d, e, roots); + } + + // Now split the depressed quartic into two quadratics: + // y^4 + py^2 + qy + r = (y^2 + sy + u)(y^2 - sy + v) = y^4 + (v+u-s^2)y^2 + s(v - u)y + uv + // So p = v+u-s^2, q = s(v - u), r = uv. + // Then (v+u)^2 - (v-u)^2 = 4uv = 4r = (p+s^2)^2 - q^2/s^2. + // Multiply through by s^2 to get s^2(p+s^2)^2 - q^2 - 4rs^2 = 0, which is a cubic in s^2. + // Then we let z = s^2, to get + // z^3 + 2pz^2 + (p^2 - 4r)z - q^2 = 0. + auto z_roots = cubic_roots(Real(1), 2*p, p*p - 4*r, -q*q); + // z = s^2, so s = sqrt(z). + // No real roots: + if (z_roots.back() <= 0) { + return roots; + } + Real s = sqrt(z_roots.back()); + + // s is nonzero, because we took care of the biquadratic case. + Real v = (p + s*s + q/s)/2; + Real u = v - q/s; + // Now solve y^2 + sy + u = 0: + auto [root0, root1] = quadratic_roots(Real(1), s, u); + + // Now solve y^2 - sy + v = 0: + auto [root2, root3] = quadratic_roots(Real(1), -s, v); + roots[0] = root0; + roots[1] = root1; + roots[2] = root2; + roots[3] = root3; + + for (auto& r : roots) { + r -= A/4; + } + return detail::polish_and_sort(a, b, c, d, e, roots); +} + +} +#endif diff --git a/reporting/performance/quartic_roots_performance.cpp b/reporting/performance/quartic_roots_performance.cpp new file mode 100644 index 000000000..fe2f3ffbc --- /dev/null +++ b/reporting/performance/quartic_roots_performance.cpp @@ -0,0 +1,41 @@ +// (C) Copyright Nick Thompson 2021. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include +#include +#include + +using boost::math::tools::quartic_roots; + +template +void QuarticRoots(benchmark::State& state) +{ + std::random_device rd; + auto seed = rd(); + // This seed generates 3 real roots: + //uint32_t seed = 416683252; + std::mt19937_64 mt(seed); + std::uniform_real_distribution unif(-10, 10); + + Real a = unif(mt); + Real b = unif(mt); + Real c = unif(mt); + Real d = unif(mt); + Real e = unif(mt); + for (auto _ : state) + { + auto roots = quartic_roots(a,b,c,d, e); + benchmark::DoNotOptimize(roots[0]); + } +} + +BENCHMARK_TEMPLATE(QuarticRoots, float); +BENCHMARK_TEMPLATE(QuarticRoots, double); +BENCHMARK_TEMPLATE(QuarticRoots, long double); + +BENCHMARK_MAIN(); diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 47c05fb54..f6b633c86 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -993,6 +993,7 @@ test-suite misc : [ run anderson_darling_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ] [ run ljung_box_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ] [ run cubic_roots_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ] + [ run quartic_roots_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ] [ run test_t_test.cpp : : : [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : BOOST_MATH_TEST_FLOAT128 "-Bstatic -lquadmath -Bdynamic" ] [ requires cxx11_hdr_forward_list cxx11_hdr_atomic cxx11_hdr_thread cxx11_hdr_tuple cxx11_hdr_future cxx11_sfinae_expr ] ] [ run test_z_test.cpp : : : [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : BOOST_MATH_TEST_FLOAT128 "-Bstatic -lquadmath -Bdynamic" ] [ requires cxx11_hdr_forward_list cxx11_hdr_atomic cxx11_hdr_thread cxx11_hdr_tuple cxx11_hdr_future cxx11_sfinae_expr ] ] [ run bivariate_statistics_test.cpp : : : [ requires cxx11_hdr_forward_list cxx11_hdr_atomic cxx11_hdr_thread cxx11_hdr_tuple cxx11_hdr_future cxx11_sfinae_expr ] ] diff --git a/test/engel_expansion_test.cpp b/test/engel_expansion_test.cpp index 6ad4d3305..0298b2d72 100644 --- a/test/engel_expansion_test.cpp +++ b/test/engel_expansion_test.cpp @@ -7,6 +7,7 @@ #include "math_unit_test.hpp" #include +#include #include #ifdef BOOST_HAS_FLOAT128 #include diff --git a/test/math_unit_test.hpp b/test/math_unit_test.hpp index 77a26fe02..b5d9755a9 100644 --- a/test/math_unit_test.hpp +++ b/test/math_unit_test.hpp @@ -13,13 +13,29 @@ #include #include #include -#include - +#if defined __has_include +# if __has_include() +#define BOOST_MATH_HAS_CXX_ABI 1 +# include +# endif +#endif namespace boost { namespace math { namespace test { namespace detail { static std::atomic global_error_count{0}; static std::atomic total_ulp_distance{0}; + + inline std::string demangle(char const * name) + { + int status = 0; + std::size_t size = 0; + std::string s = name; +#if BOOST_MATH_HAS_CXX_ABI + return abi::__cxa_demangle( name, NULL, &size, &status ); +#else + return name; +#endif + } } template @@ -49,7 +65,7 @@ bool check_mollified_close(Real expected, Real computed, Real tol, std::string c std::ios_base::fmtflags f( std::cerr.flags() ); std::cerr << std::setprecision(3); std::cerr << "\033[0;31mError at " << filename << ":" << function << ":" << line << ":\n" - << " \033[0m Mollified relative error in " << boost::core::demangle(typeid(Real).name())<< " precision is " << mollified_relative_error + << " \033[0m Mollified relative error in " << detail::demangle(typeid(Real).name())<< " precision is " << mollified_relative_error << ", which exceeds " << tol << ", error/tol = " << mollified_relative_error/tol << ".\n" << std::setprecision(std::numeric_limits::max_digits10) << std::showpos << " Expected: " << std::defaultfloat << std::fixed << expected << std::hexfloat << " = " << expected << "\n" @@ -77,8 +93,8 @@ bool check_ulp_close(PreciseReal expected1, Real computed, size_t ulps, std::str if (sizeof(PreciseReal) < sizeof(Real)) { std::ostringstream err; err << "\n\tThe expected number must be computed in higher (or equal) precision than the number being tested.\n"; - err << "\tType of expected is " << boost::core::demangle(typeid(PreciseReal).name()) << ", which occupies " << sizeof(PreciseReal) << " bytes.\n"; - err << "\tType of computed is " << boost::core::demangle(typeid(Real).name()) << ", which occupies " << sizeof(Real) << " bytes.\n"; + err << "\tType of expected is " << detail::demangle(typeid(PreciseReal).name()) << ", which occupies " << sizeof(PreciseReal) << " bytes.\n"; + err << "\tType of computed is " << detail::demangle(typeid(Real).name()) << ", which occupies " << sizeof(Real) << " bytes.\n"; throw std::logic_error(err.str()); } } @@ -105,7 +121,7 @@ bool check_ulp_close(PreciseReal expected1, Real computed, size_t ulps, std::str std::ios_base::fmtflags f( std::cerr.flags() ); std::cerr << std::setprecision(3); std::cerr << "\033[0;31mError at " << filename << ":" << function << ":" << line << ":\n" - << " \033[0m ULP distance in " << boost::core::demangle(typeid(Real).name())<< " precision is " << dist + << " \033[0m ULP distance in " << detail::demangle(typeid(Real).name())<< " precision is " << dist << ", which exceeds " << ulps; if (ulps > 0) { @@ -161,7 +177,7 @@ bool check_le(Real lesser, Real greater, std::string const & filename, std::stri std::ios_base::fmtflags f( std::cerr.flags() ); std::cerr << std::setprecision(3); std::cerr << "\033[0;31mError at " << filename << ":" << function << ":" << line << ":\n" - << " \033[0m Condition " << lesser << " \u2264 " << greater << " is violated in " << boost::core::demangle(typeid(Real).name()) << " precision.\n"; + << " \033[0m Condition " << lesser << " \u2264 " << greater << " is violated in " << detail::demangle(typeid(Real).name()) << " precision.\n"; std::cerr << std::setprecision(std::numeric_limits::max_digits10) << std::showpos << " \"Lesser\" : " << std::defaultfloat << std::fixed << lesser << " = " << std::scientific << lesser << std::hexfloat << " = " << lesser << "\n" << " \"Greater\": " << std::defaultfloat << std::fixed << greater << " = " << std::scientific << greater << std::hexfloat << " = " << greater << "\n" @@ -214,7 +230,7 @@ bool check_conditioned_error(Real abscissa, PreciseReal expected1, PreciseReal e std::cerr << "\033[0;31mError at " << filename << ":" << function << ":" << line << ":\n"; std::cerr << std::setprecision(std::numeric_limits::max_digits10) << std::showpos; std::cerr << "\033[0m Error at abscissa " << std::defaultfloat << std::fixed << abscissa << " = " << std::hexfloat << abscissa << "\n"; - std::cerr << " Given that the expected value is zero, the computed value in " << boost::core::demangle(typeid(Real).name()) << " precision must satisfy |f(x)| <= " << tol << ".\n"; + std::cerr << " Given that the expected value is zero, the computed value in " << detail::demangle(typeid(Real).name()) << " precision must satisfy |f(x)| <= " << tol << ".\n"; std::cerr << " But the computed value is " << std::defaultfloat << std::fixed << computed << std::hexfloat << " = " << computed << "\n"; std::cerr.flags(f); ++detail::global_error_count; @@ -240,7 +256,7 @@ bool check_conditioned_error(Real abscissa, PreciseReal expected1, PreciseReal e std::cerr << "\033[0;31mError at " << filename << ":" << function << ":" << line << "\n"; std::cerr << std::setprecision(std::numeric_limits::max_digits10); std::cerr << "\033[0m The relative error at abscissa x = " << std::defaultfloat << std::fixed << abscissa << " = " << std::hexfloat << abscissa - << " in " << boost::core::demangle(typeid(Real).name()) << " precision is " << std::scientific << relative_error << "\n" + << " in " << detail::demangle(typeid(Real).name()) << " precision is " << std::scientific << relative_error << "\n" << " This exceeds the tolerance " << tol << "\n" << std::showpos << " Expected: " << std::defaultfloat << std::fixed << expected << " = " << std::scientific << expected << std::hexfloat << " = " << expected << "\n" @@ -288,7 +304,7 @@ bool check_absolute_error(PreciseReal expected1, Real computed, Real acceptable_ std::cerr << std::setprecision(3); std::cerr << "\033[0;31mError at " << filename << ":" << function << ":" << line << "\n"; std::cerr << std::setprecision(std::numeric_limits::max_digits10); - std::cerr << "\033[0m The absolute error in " << boost::core::demangle(typeid(Real).name()) << " precision is " << std::scientific << error << "\n" + std::cerr << "\033[0m The absolute error in " << detail::demangle(typeid(Real).name()) << " precision is " << std::scientific << error << "\n" << " This exceeds the acceptable error " << acceptable_error << "\n" << std::showpos << " Expected: " << std::defaultfloat << std::fixed << expected << " = " << std::scientific << expected << std::hexfloat << " = " << expected << "\n" diff --git a/test/quartic_roots_test.cpp b/test/quartic_roots_test.cpp new file mode 100644 index 000000000..bc456ce8c --- /dev/null +++ b/test/quartic_roots_test.cpp @@ -0,0 +1,137 @@ +/* + * Copyright Nick Thompson, 2021 + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0. (See accompanying file + * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + */ + +#include "math_unit_test.hpp" +#include +#include +#ifdef BOOST_HAS_FLOAT128 +#include +using boost::multiprecision::float128; +#endif + +using boost::math::tools::quartic_roots; +using std::cbrt; +using std::sqrt; + +template +void test_zero_coefficients() +{ + Real a = 0; + Real b = 0; + Real c = 0; + Real d = 0; + Real e = 0; + auto roots = quartic_roots(a,b,c,d,e); + CHECK_EQUAL(roots[0], Real(0)); + CHECK_EQUAL(roots[1], Real(0)); + CHECK_EQUAL(roots[2], Real(0)); + CHECK_EQUAL(roots[3], Real(0)); + + b = 1; + e = 1; + // x^3 + 1 = 0: + roots = quartic_roots(a,b,c,d,e); + CHECK_EQUAL(roots[0], Real(-1)); + CHECK_NAN(roots[1]); + CHECK_NAN(roots[2]); + CHECK_NAN(roots[3]); + e = -1; + // x^3 - 1 = 0: + roots = quartic_roots(a,b,c,d,e); + CHECK_EQUAL(roots[0], Real(1)); + CHECK_NAN(roots[1]); + CHECK_NAN(roots[2]); + CHECK_NAN(roots[3]); + + e = -2; + // x^3 - 2 = 0 + roots = quartic_roots(a,b,c,d,e); + CHECK_ULP_CLOSE(roots[0], cbrt(Real(2)), 2); + CHECK_NAN(roots[1]); + CHECK_NAN(roots[2]); + CHECK_NAN(roots[3]); + + // x^4 -1 = 0 + // x = \pm 1: + roots = quartic_roots(1, 0, 0, 0, -1); + CHECK_ULP_CLOSE(Real(-1), roots[0], 3); + CHECK_ULP_CLOSE(Real(1), roots[1], 3); + CHECK_NAN(roots[2]); + CHECK_NAN(roots[3]); + + // x^4 - 2 = 0 \implies x = \pm sqrt(sqrt(2)) + roots = quartic_roots(1,0,0,0,-2); + CHECK_ULP_CLOSE(-sqrt(sqrt(Real(2))), roots[0], 3); + CHECK_ULP_CLOSE(sqrt(sqrt(Real(2))), roots[1], 3); + CHECK_NAN(roots[2]); + CHECK_NAN(roots[3]); + + + // x(x-1)(x-2)(x-3) = x^4 - 6x^3 + 11x^2 - 6x + roots = quartic_roots(Real(1), Real(-6), Real(11), Real(-6), Real(0)); + CHECK_ULP_CLOSE(roots[0], Real(0), 2); + CHECK_ULP_CLOSE(roots[1], Real(1), 2); + CHECK_ULP_CLOSE(roots[2], Real(2), 2); + CHECK_ULP_CLOSE(roots[3], Real(3), 2); + + // (x-1)(x-2)(x-3)(x-4) = x^4 - 10x^3 + 35x^2 - (2*3*4 + 1*3*4 + 1*2*4 + 1*2*3)x + 1*2*3*4 + roots = quartic_roots(1, -10, 35, -24 - 12 - 8 - 6, 1*2*3*4); + CHECK_ULP_CLOSE(Real(1), roots[0], 2); + CHECK_ULP_CLOSE(Real(2), roots[1], 2); + CHECK_ULP_CLOSE(Real(3), roots[2], 2); + CHECK_ULP_CLOSE(Real(4), roots[3], 2); + + // Double root: + // (x+1)^2(x-2)(x-3) = x^4 - 3x^3 -3x^2 + 7x + 6 + // Note: This test is unstable wrt to perturbations! + roots = quartic_roots(Real(1), Real(-3), Real(-3), Real(7), Real(6)); + CHECK_ULP_CLOSE(Real(-1), roots[0], 2); + CHECK_ULP_CLOSE(Real(-1), roots[1], 2); + CHECK_ULP_CLOSE(Real(2), roots[2], 2); + CHECK_ULP_CLOSE(Real(3), roots[3], 2); + + + std::uniform_real_distribution dis(-2,2); + std::mt19937 gen(12343); + // Expected roots + std::array r; + int trials = 10; + for (int i = 0; i < trials; ++i) { + // Mathematica: + // Expand[(x - r0)*(x - r1)*(x - r2)*(x-r3)] + // r0 r1 r2 r3 - (r0 r1 r2 + r0 r1 r3 + r0 r2 r3 + r1r2r3)x + // + (r0 r1 + r0 r2 + r0 r3 + r1 r2 + r1r3 + r2 r3)x^2 - (r0 + r1 + r2 + r3) x^3 + x^4 + for (auto & root : r) { + root = static_cast(dis(gen)); + } + std::sort(r.begin(), r.end()); + Real a = 1; + Real b = -(r[0] + r[1] + r[2] + r[3]); + Real c = r[0]*r[1] + r[0]*r[2] + r[0]*r[3] + r[1]*r[2] + r[1]*r[3] + r[2]*r[3]; + Real d = -(r[0]*r[1]*r[2] + r[0]*r[1]*r[3] + r[0]*r[2]*r[3] + r[1]*r[2]*r[3]); + Real e = r[0]*r[1]*r[2]*r[3]; + + auto roots = quartic_roots(a, b, c, d, e); + // I could check the condition number here, but this is fine right? + CHECK_ULP_CLOSE(r[0], roots[0], 160); + CHECK_ULP_CLOSE(r[1], roots[1], 260); + CHECK_ULP_CLOSE(r[2], roots[2], 160); + CHECK_ULP_CLOSE(r[3], roots[3], 160); + } +} + + +int main() +{ + test_zero_coefficients(); + test_zero_coefficients(); +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + test_zero_coefficients(); +#endif + + return boost::math::test::report_errors(); +} diff --git a/test/simple_continued_fraction_test.cpp b/test/simple_continued_fraction_test.cpp index 637916496..fedb4e5a3 100644 --- a/test/simple_continued_fraction_test.cpp +++ b/test/simple_continued_fraction_test.cpp @@ -8,6 +8,7 @@ #include "math_unit_test.hpp" #include #include +#include #ifdef BOOST_HAS_FLOAT128 #include using boost::multiprecision::float128; diff --git a/test/test_z_test.cpp b/test/test_z_test.cpp index 1acbde8a0..3f4ad9940 100644 --- a/test/test_z_test.cpp +++ b/test/test_z_test.cpp @@ -13,6 +13,7 @@ #include using quad = boost::multiprecision::cpp_bin_float_quad; +using std::sqrt; template void test_one_sample_z() @@ -56,7 +57,7 @@ void test_two_sample_z() Real computed_statistic = std::get<0>(temp); Real computed_pvalue = std::get<1>(temp); CHECK_ULP_CLOSE(Real(1), computed_statistic, 5); - CHECK_MOLLIFIED_CLOSE(Real(0), computed_pvalue, 5*std::numeric_limits::epsilon()); + CHECK_MOLLIFIED_CLOSE(Real(0), computed_pvalue, sqrt(std::numeric_limits::epsilon())); } template