boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Quartic roots. (#718)

Browse Source
This commit is contained in:
Nick
2022-01-02 17:58:09 -08:00
committed by GitHub
parent 4e74c1251e
commit 60d54e565f
10 changed files with 400 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files

40
doc/roots/quartic_roots.qbk Normal file
View File

@@ -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 <boost/math/roots/quartic_roots.hpp>
namespace boost::math::tools {
// Solves ax⁴ + bx³ + cx² + dx + e = 0.
std::array<Real,3> 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<Real, 4>`, 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]

1
doc/roots/roots_overview.qbk
View File

@@ -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]

150
include/boost/math/tools/quartic_roots.hpp Normal file
View File

@@ -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 <array>
#include <cmath>
#include <boost/math/tools/cubic_roots.hpp>
namespace boost::math::tools {
namespace detail {
// Make sure the nans are always at the back of the array:
template<typename Real>
bool comparator(Real r1, Real r2) {
using std::isnan;
if (isnan(r2)) { return true; }
return r1 < r2;
}
template<typename Real>
std::array<Real, 4> polish_and_sort(Real a, Real b, Real c, Real d, Real e, std::array<Real, 4>& 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<Real>::min)())
{
r -= 2*f*df/denom;
}
}
std::sort(roots.begin(), roots.end(), detail::comparator<Real>);
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<typename Real>
std::array<Real, 4> quartic_roots(Real a, Real b, Real c, Real d, Real e) {
using std::abs;
using std::sqrt;
auto nan = std::numeric_limits<Real>::quiet_NaN();
std::array<Real, 4> roots{nan, nan, nan, nan};
if (abs(a) <= (std::numeric_limits<Real>::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<Real>::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<Real>::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<Real>::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

41
reporting/performance/quartic_roots_performance.cpp Normal file
View File

@@ -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 <random>
#include <array>
#include <vector>
#include <iostream>
#include <benchmark/benchmark.h>
#include <boost/math/tools/quartic_roots.hpp>
using boost::math::tools::quartic_roots;
template<class Real>
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<Real> 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();

1
test/Jamfile.v2
View File

@@ -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" : <define>BOOST_MATH_TEST_FLOAT128 <linkflags>"-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" : <define>BOOST_MATH_TEST_FLOAT128 <linkflags>"-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 ] ]

1
test/engel_expansion_test.cpp
View File

@@ -7,6 +7,7 @@
#include "math_unit_test.hpp"
#include <boost/math/tools/engel_expansion.hpp>
#include <boost/core/demangle.hpp>
#include <boost/math/constants/constants.hpp>
#ifdef BOOST_HAS_FLOAT128
#include <boost/multiprecision/float128.hpp>

36
test/math_unit_test.hpp
View File

@@ -13,13 +13,29 @@
#include <boost/math/tools/assert.hpp>
#include <boost/math/special_functions/next.hpp>
#include <boost/math/special_functions/trunc.hpp>
#include <boost/core/demangle.hpp>
#if defined __has_include
# if __has_include(<cxxabi.h>)
#define BOOST_MATH_HAS_CXX_ABI 1
# include <cxxabi.h>
# endif
#endif
namespace boost { namespace math { namespace test {
namespace detail {
static std::atomic<int64_t> global_error_count{0};
static std::atomic<int64_t> 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<class Real>
@@ -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<Real>::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<Real>::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<Real>::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<Real>::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<Real>::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"

137
test/quartic_roots_test.cpp Normal file
View File

@@ -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 <random>
#include <boost/math/tools/quartic_roots.hpp>
#ifdef BOOST_HAS_FLOAT128
#include <boost/multiprecision/float128.hpp>
using boost::multiprecision::float128;
#endif
using boost::math::tools::quartic_roots;
using std::cbrt;
using std::sqrt;
template<class Real>
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<Real>(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<Real>(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<Real>(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<Real> dis(-2,2);
std::mt19937 gen(12343);
// Expected roots
std::array<Real, 4> 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<Real>(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<float>();
test_zero_coefficients<double>();
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
test_zero_coefficients<long double>();
#endif
return boost::math::test::report_errors();
}

1
test/simple_continued_fraction_test.cpp
View File

@@ -8,6 +8,7 @@
#include "math_unit_test.hpp"
#include <boost/math/tools/simple_continued_fraction.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/core/demangle.hpp>
#ifdef BOOST_HAS_FLOAT128
#include <boost/multiprecision/float128.hpp>
using boost::multiprecision::float128;

3
test/test_z_test.cpp
View File

@@ -13,6 +13,7 @@
#include <utility>
using quad = boost::multiprecision::cpp_bin_float_quad;
using std::sqrt;
template<typename Real>
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<Real>::epsilon());
CHECK_MOLLIFIED_CLOSE(Real(0), computed_pvalue, sqrt(std::numeric_limits<Real>::epsilon()));
}
template<typename Z>