mirror of
https://github.com/boostorg/math.git
synced 2026-01-19 04:22:09 +00:00
Merge branch 'legendre_stieltjes' of https://github.com/NAThompson/math into stieltjes
This commit is contained in:
jzmaddock
71
doc/sf/legendre_stieltjes.qbk
Normal file
71
doc/sf/legendre_stieltjes.qbk
Normal file
@@ -0,0 +1,71 @@
|
||||
[section:legendre_stieltjes Legendre-Stieltjes Polynomials]
|
||||
|
||||
[h4 Synopsis]
|
||||
|
||||
``
|
||||
#include <boost/math/special_functions/legendre_stieltjes.hpp>
|
||||
``
|
||||
|
||||
namespace boost{ namespace math{
|
||||
|
||||
template <class T>
|
||||
class legendre_stieltjes
|
||||
{
|
||||
public:
|
||||
legendre_stieltjes(size_t m);
|
||||
|
||||
Real norm_sq() const;
|
||||
|
||||
Real operator()(Real x) const;
|
||||
|
||||
Real prime(Real x) const;
|
||||
|
||||
std::vector<Real> zeros() const;
|
||||
}
|
||||
|
||||
}}
|
||||
|
||||
[h4 Description]
|
||||
|
||||
The Legendre-Stieltjes polynomials are a family of polynomials used to generate Gauss-Konrod quadrature formulas.
|
||||
Gauss-Konrod quadratures are algorithms which extend a Gaussian quadrature in such a way that all abscissas
|
||||
are reused when computed a higher-order estimate of the integral, allowing efficient calculation of an error estimate.
|
||||
The Legendre-Stieltjes polynomials assist with this task because their zeros /interlace/ the zeros of the Legendre polynomials,
|
||||
meaning that between any two zeros of a Legendre polynomial of degree n, there exists a zero of the Legendre-Stieltjes polynomial
|
||||
of degree n+1.
|
||||
|
||||
The Legendre-Stieltjes polynomials /E/[sub n+1] are defined by the property that they have /n/ vanishing moments against the oscillatory measure P[sub n], i.e., \u222B[sub -1][super 1] E[sub n+1](x)P[sub n](x) x[super k] dx = 0 for /k = 0, 1, ..., n/.
|
||||
The first few are
|
||||
|
||||
* E[sub 1](x) = P[sub 1](x)
|
||||
* E[sub 2](x) = P[sub 2](x) - 2P[sub 0](x)/5
|
||||
* E[sub 3](x) = P[sub 3](x) - 9P[sub 1](x)/14
|
||||
* E[sub 4](x) = P[sub 4](x) - 20P[sub 2](x)/27 + 14P[sub 0](x)/891
|
||||
* E[sub 5](x) = P[sub 5](x) - 35P[sub 3](x)/44 + 135P[sub 1](x)/12584
|
||||
|
||||
where P[sub i] are the Legendre polynomials.
|
||||
The scaling follows [@http://www.ams.org/journals/mcom/1968-22-104/S0025-5718-68-99866-9/S0025-5718-68-99866-9.pdf Patterson],
|
||||
who expanded the Legendre-Stieltjes polynomials in a Legendre series and took the coefficient of the highest-order Legendre polynomial in the series to be unity.
|
||||
|
||||
The Legendre-Stieltjes polynomials do not satisfy three-term recurrence relations or have a particulary simple representation.
|
||||
Hence the constructor call determines what, in fact, the polynomial is.
|
||||
Once the constructor comes back, the polynomial can be evaluated via the Legendre series.
|
||||
|
||||
Example usage:
|
||||
|
||||
// Call to the constructor determines the coefficients in the Legendre expansion
|
||||
legendre_stieltjes<double> E(12);
|
||||
// Evaluate the polynomial at a point:
|
||||
double x = E(0.3);
|
||||
// Evaluate the derivative at a point:
|
||||
double x_p = E.prime(0.3);
|
||||
// Use the norm_sq to change between scalings, if desired:
|
||||
double norm = std::sqrt(E.norm_sq());
|
||||
|
||||
[endsect]
|
||||
[/
|
||||
Copyright 2017 Nick Thompson
|
||||
Distributed under 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).
|
||||
]
|
||||
100
example/legendre_stieltjes_example.cpp
Normal file
100
example/legendre_stieltjes_example.cpp
Normal file
@@ -0,0 +1,100 @@
|
||||
// Copyright Nick Thompson 2017.
|
||||
// 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 <iostream>
|
||||
#include <string>
|
||||
#include <boost/math/constants/constants.hpp>
|
||||
#include <boost/multiprecision/cpp_bin_float.hpp>
|
||||
#include <boost/multiprecision/cpp_dec_float.hpp>
|
||||
#include <boost/math/special_functions/legendre.hpp>
|
||||
#include <boost/math/special_functions/legendre_stieltjes.hpp>
|
||||
|
||||
using boost::math::legendre_p;
|
||||
using boost::math::legendre_p_zeros;
|
||||
using boost::math::legendre_p_prime;
|
||||
using boost::math::legendre_stieltjes;
|
||||
using boost::multiprecision::cpp_bin_float_quad;
|
||||
using boost::multiprecision::cpp_dec_float_100;
|
||||
|
||||
template<class Real>
|
||||
void gauss_konrod_rule(size_t order)
|
||||
{
|
||||
std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
|
||||
std::cout << std::fixed;
|
||||
auto gauss_nodes = boost::math::legendre_p_zeros<Real>(order);
|
||||
auto E = legendre_stieltjes<Real>(order + 1);
|
||||
std::vector<Real> gauss_weights(gauss_nodes.size(), std::numeric_limits<Real>::quiet_NaN());
|
||||
std::vector<Real> gauss_konrod_weights(gauss_nodes.size(), std::numeric_limits<Real>::quiet_NaN());
|
||||
for (size_t i = 0; i < gauss_nodes.size(); ++i)
|
||||
{
|
||||
Real node = gauss_nodes[i];
|
||||
Real lp = legendre_p_prime<Real>(order, node);
|
||||
gauss_weights[i] = 2/( (1-node*node)*lp*lp);
|
||||
// P_n(x) = (2n)!/(2^n (n!)^2) pi_n(x), where pi_n is the monic Legendre polynomial.
|
||||
gauss_konrod_weights[i] = gauss_weights[i] + static_cast<Real>(2)/(static_cast<Real>(order+1)*legendre_p_prime(order, node)*E(node));
|
||||
}
|
||||
|
||||
std::cout << "Gauss Nodes:\n";
|
||||
for (auto const & node : gauss_nodes)
|
||||
{
|
||||
std::cout << node << "\n";
|
||||
}
|
||||
|
||||
std::cout << "Gauss Weights:\n";
|
||||
for (auto const & weight : gauss_weights)
|
||||
{
|
||||
std::cout << weight << "\n";
|
||||
}
|
||||
|
||||
std::cout << "Gauss-Konrod weights: \n";
|
||||
for (auto const & w : gauss_konrod_weights)
|
||||
{
|
||||
std::cout << w << "\n";
|
||||
}
|
||||
|
||||
auto konrod_nodes = E.zeros();
|
||||
std::vector<Real> konrod_weights(konrod_nodes.size());
|
||||
for (size_t i = 0; i < konrod_weights.size(); ++i)
|
||||
{
|
||||
Real node = konrod_nodes[i];
|
||||
konrod_weights[i] = static_cast<Real>(2)/(static_cast<Real>(order+1)*legendre_p(order, node)*E.prime(node));
|
||||
}
|
||||
|
||||
std::cout << "Konrod nodes:\n";
|
||||
for (auto node : konrod_nodes)
|
||||
{
|
||||
std::cout << node << "\n";
|
||||
}
|
||||
|
||||
std::cout << "Konrod weights: \n";
|
||||
for (auto const & w : gauss_konrod_weights)
|
||||
{
|
||||
std::cout << w << "\n";
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
std::cout << "Gauss-Konrod 7-15 Rule:\n";
|
||||
gauss_konrod_rule<cpp_dec_float_100>(7);
|
||||
|
||||
std::cout << "\n\nGauss-Konrod 10-21 Rule:\n";
|
||||
gauss_konrod_rule<cpp_dec_float_100>(10);
|
||||
|
||||
std::cout << "\n\nGauss-Konrod 15-31 Rule:\n";
|
||||
gauss_konrod_rule<cpp_dec_float_100>(15);
|
||||
|
||||
std::cout << "\n\nGauss-Konrod 20-41 Rule:\n";
|
||||
gauss_konrod_rule<cpp_dec_float_100>(20);
|
||||
|
||||
std::cout << "\n\nGauss-Konrod 25-51 Rule:\n";
|
||||
gauss_konrod_rule<cpp_dec_float_100>(25);
|
||||
|
||||
std::cout << "\n\nGauss-Konrod 30-61 Rule:\n";
|
||||
gauss_konrod_rule<cpp_dec_float_100>(30);
|
||||
|
||||
}
|
||||
235
include/boost/math/special_functions/legendre_stieltjes.hpp
Normal file
235
include/boost/math/special_functions/legendre_stieltjes.hpp
Normal file
@@ -0,0 +1,235 @@
|
||||
// Copyright Nick Thompson 2017.
|
||||
// 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_SPECIAL_LEGENDRE_STIELTJES_HPP
|
||||
#define BOOST_MATH_SPECIAL_LEGENDRE_STIELTJES_HPP
|
||||
|
||||
/*
|
||||
* Constructs the Legendre-Stieltjes polynomial of degree m.
|
||||
* The Legendre-Stieltjes polynomials are used to create extensions for Gaussian quadratures,
|
||||
* commonly called "Gauss-Konrod" quadratures.
|
||||
*
|
||||
* References:
|
||||
* Patterson, TNL. "The optimum addition of points to quadrature formulae." Mathematics of Computation 22.104 (1968): 847-856.
|
||||
*/
|
||||
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
#include <boost/math/tools/roots.hpp>
|
||||
#include <boost/multiprecision/cpp_int.hpp>
|
||||
|
||||
namespace boost{
|
||||
namespace math{
|
||||
|
||||
template<class Real>
|
||||
class legendre_stieltjes
|
||||
{
|
||||
public:
|
||||
legendre_stieltjes(size_t m)
|
||||
{
|
||||
if (m == 0)
|
||||
{
|
||||
throw std::domain_error("The Legendre-Stieltjes polynomial is defined for order m > 0.\n");
|
||||
}
|
||||
m_m = m;
|
||||
int64_t n = m - 1;
|
||||
int64_t q;
|
||||
int64_t r;
|
||||
bool odd = n & 1;
|
||||
if (odd)
|
||||
{
|
||||
q = 1;
|
||||
r = (n-1)/2 + 2;
|
||||
}
|
||||
else
|
||||
{
|
||||
q = 0;
|
||||
r = n/2 + 1;
|
||||
}
|
||||
m_a.resize(r + 1);
|
||||
// We'll keep the ones-based indexing at the cost of storing a superfluous element
|
||||
// so that we can follow Patterson's notation exactly.
|
||||
m_a[r] = static_cast<Real>(1);
|
||||
// Make sure using the zero index is a bug:
|
||||
m_a[0] = std::numeric_limits<Real>::quiet_NaN();
|
||||
|
||||
for (int64_t k = 1; k < r; ++k)
|
||||
{
|
||||
Real ratio = 1;
|
||||
m_a[r - k] = 0;
|
||||
for (int64_t i = r + 1 - k; i <= r; ++i)
|
||||
{
|
||||
// See Patterson, equation 12
|
||||
int64_t num = (n - q + 2*(i + k - 1))*(n + q + 2*(k - i + 1))*(n-1-q+2*(i-k))*(2*(k+i-1) -1 -q -n);
|
||||
int64_t den = (n - q + 2*(i - k))*(2*(k + i - 1) - q - n)*(n + 1 + q + 2*(k - i))*(n - 1 - q + 2*(i + k));
|
||||
ratio *= static_cast<Real>(num)/static_cast<Real>(den);
|
||||
m_a[r - k] -= ratio*m_a[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
Real norm_sq() const
|
||||
{
|
||||
Real t = 0;
|
||||
bool odd = m_m & 1;
|
||||
for (size_t i = 1; i < m_a.size(); ++i)
|
||||
{
|
||||
if(odd)
|
||||
{
|
||||
t += 2*m_a[i]*m_a[i]/static_cast<Real>(4*i-1);
|
||||
}
|
||||
else
|
||||
{
|
||||
t += 2*m_a[i]*m_a[i]/static_cast<Real>(4*i-3);
|
||||
}
|
||||
}
|
||||
return t;
|
||||
}
|
||||
|
||||
|
||||
Real operator()(Real x) const
|
||||
{
|
||||
// Trivial implementation:
|
||||
// Em += m_a[i]*legendre_p(2*i - 1, x); m odd
|
||||
// Em += m_a[i]*legendre_p(2*i - 2, x); m even
|
||||
size_t r = m_a.size() - 1;
|
||||
Real p0 = 1;
|
||||
Real p1 = x;
|
||||
|
||||
Real Em;
|
||||
bool odd = m_m & 1;
|
||||
if (odd)
|
||||
{
|
||||
Em = m_a[1]*p1;
|
||||
}
|
||||
else
|
||||
{
|
||||
Em = m_a[1]*p0;
|
||||
}
|
||||
|
||||
unsigned n = 1;
|
||||
for (size_t i = 2; i <= r; ++i)
|
||||
{
|
||||
std::swap(p0, p1);
|
||||
p1 = boost::math::legendre_next(n, x, p0, p1);
|
||||
++n;
|
||||
if (!odd)
|
||||
{
|
||||
Em += m_a[i]*p1;
|
||||
}
|
||||
std::swap(p0, p1);
|
||||
p1 = boost::math::legendre_next(n, x, p0, p1);
|
||||
++n;
|
||||
if(odd)
|
||||
{
|
||||
Em += m_a[i]*p1;
|
||||
}
|
||||
}
|
||||
return Em;
|
||||
}
|
||||
|
||||
|
||||
Real prime(Real x) const
|
||||
{
|
||||
Real Em_prime = 0;
|
||||
|
||||
for (size_t i = 1; i < m_a.size(); ++i)
|
||||
{
|
||||
if(m_m & 1)
|
||||
{
|
||||
Em_prime += m_a[i]*detail::legendre_p_prime_imp(2*i - 1, x, policies::policy<>());
|
||||
}
|
||||
else
|
||||
{
|
||||
Em_prime += m_a[i]*detail::legendre_p_prime_imp(2*i - 2, x, policies::policy<>());
|
||||
}
|
||||
}
|
||||
return Em_prime;
|
||||
}
|
||||
|
||||
std::vector<Real> zeros() const
|
||||
{
|
||||
using boost::math::constants::half;
|
||||
|
||||
std::vector<Real> stieltjes_zeros;
|
||||
std::vector<Real> legendre_zeros = legendre_p_zeros<Real>(m_m - 1);
|
||||
int k;
|
||||
if (m_m & 1)
|
||||
{
|
||||
stieltjes_zeros.resize(legendre_zeros.size() + 1, std::numeric_limits<Real>::quiet_NaN());
|
||||
stieltjes_zeros[0] = 0;
|
||||
k = 1;
|
||||
}
|
||||
else
|
||||
{
|
||||
stieltjes_zeros.resize(legendre_zeros.size(), std::numeric_limits<Real>::quiet_NaN());
|
||||
k = 0;
|
||||
}
|
||||
|
||||
while (k < stieltjes_zeros.size())
|
||||
{
|
||||
Real lower_bound;
|
||||
Real upper_bound;
|
||||
if (m_m & 1)
|
||||
{
|
||||
lower_bound = legendre_zeros[k - 1];
|
||||
if (k == legendre_zeros.size())
|
||||
{
|
||||
upper_bound = 1;
|
||||
}
|
||||
else
|
||||
{
|
||||
upper_bound = legendre_zeros[k];
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
lower_bound = legendre_zeros[k];
|
||||
if (k == legendre_zeros.size() - 1)
|
||||
{
|
||||
upper_bound = 1;
|
||||
}
|
||||
else
|
||||
{
|
||||
upper_bound = legendre_zeros[k+1];
|
||||
}
|
||||
}
|
||||
|
||||
// The root bracketing is not very tight; to keep weird stuff from happening
|
||||
// in the Newton's method, let's tighten up the tolerance using a few bisections.
|
||||
boost::math::tools::eps_tolerance<Real> tol(6);
|
||||
auto g = [&](Real t) { return this->operator()(t); };
|
||||
auto p = boost::math::tools::bisect(g, lower_bound, upper_bound, tol);
|
||||
|
||||
Real x_nk_guess = p.first + (p.second - p.first)*half<Real>();
|
||||
boost::uintmax_t number_of_iterations = 500;
|
||||
|
||||
auto f = [&] (Real x) { Real Pn = this->operator()(x);
|
||||
Real Pn_prime = this->prime(x);
|
||||
return std::pair<Real, Real>(Pn, Pn_prime); };
|
||||
|
||||
const Real x_nk = boost::math::tools::newton_raphson_iterate(f, x_nk_guess,
|
||||
p.first, p.second,
|
||||
2*std::numeric_limits<Real>::digits10,
|
||||
number_of_iterations);
|
||||
|
||||
BOOST_ASSERT(p.first < x_nk);
|
||||
BOOST_ASSERT(x_nk < p.second);
|
||||
stieltjes_zeros[k] = x_nk;
|
||||
++k;
|
||||
}
|
||||
return stieltjes_zeros;
|
||||
}
|
||||
|
||||
private:
|
||||
// Coefficients of Legendre expansion
|
||||
std::vector<Real> m_a;
|
||||
int m_m;
|
||||
};
|
||||
|
||||
}}
|
||||
#endif
|
||||
141
test/legendre_stieltjes_test.cpp
Normal file
141
test/legendre_stieltjes_test.cpp
Normal file
@@ -0,0 +1,141 @@
|
||||
// Copyright Nick Thompson 2017.
|
||||
// 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)
|
||||
|
||||
#define BOOST_TEST_MODULE LegendreStieltjesTest
|
||||
#define BOOST_TEST_DYN_LINK
|
||||
#include <boost/test/unit_test.hpp>
|
||||
#include <boost/math/special_functions/legendre.hpp>
|
||||
#include <boost/math/special_functions/legendre_stieltjes.hpp>
|
||||
#include <boost/math/constants/constants.hpp>
|
||||
#include <boost/multiprecision/cpp_bin_float.hpp>
|
||||
|
||||
|
||||
using boost::math::legendre_stieltjes;
|
||||
using boost::math::legendre_p;
|
||||
using boost::multiprecision::cpp_bin_float_quad;
|
||||
|
||||
|
||||
template<class Real>
|
||||
void test_legendre_stieltjes()
|
||||
{
|
||||
std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
|
||||
using std::sqrt;
|
||||
using std::abs;
|
||||
using boost::math::constants::third;
|
||||
using boost::math::constants::half;
|
||||
|
||||
Real tol = std::numeric_limits<Real>::epsilon();
|
||||
legendre_stieltjes<Real> ls1(1);
|
||||
legendre_stieltjes<Real> ls2(2);
|
||||
legendre_stieltjes<Real> ls3(3);
|
||||
legendre_stieltjes<Real> ls4(4);
|
||||
legendre_stieltjes<Real> ls5(5);
|
||||
legendre_stieltjes<Real> ls8(8);
|
||||
Real x = -1;
|
||||
while(x <= 1)
|
||||
{
|
||||
BOOST_CHECK_CLOSE_FRACTION(ls1(x), x, tol);
|
||||
BOOST_CHECK_CLOSE_FRACTION(ls1.prime(x), 1, tol);
|
||||
|
||||
Real p2 = legendre_p(2, x);
|
||||
BOOST_CHECK_CLOSE_FRACTION(ls2(x), p2 - 2/static_cast<Real>(5), tol);
|
||||
BOOST_CHECK_CLOSE_FRACTION(ls2.prime(x), 3*x, tol);
|
||||
|
||||
Real p3 = legendre_p(3, x);
|
||||
BOOST_CHECK_CLOSE_FRACTION(ls3(x), p3 - 9*x/static_cast<Real>(14), 100*tol);
|
||||
BOOST_CHECK_CLOSE_FRACTION(ls3.prime(x), 15*x*x*half<Real>() -3*half<Real>()-9/static_cast<Real>(14), 100*tol);
|
||||
|
||||
Real p4 = legendre_p(4, x);
|
||||
//-20P_2(x)/27 + 14P_0(x)/891
|
||||
Real E4 = p4 - 20*p2/static_cast<Real>(27) + 14/static_cast<Real>(891);
|
||||
BOOST_CHECK_CLOSE_FRACTION(ls4(x), E4, 250*tol);
|
||||
BOOST_CHECK_CLOSE_FRACTION(ls4.prime(x), 35*x*(9*x*x -5)/static_cast<Real>(18), 250*tol);
|
||||
|
||||
Real p5 = legendre_p(5, x);
|
||||
Real E5 = p5 - 35*p3/static_cast<Real>(44) + 135*x/static_cast<Real>(12584);
|
||||
BOOST_CHECK_CLOSE_FRACTION(ls5(x), E5, 29000*tol);
|
||||
Real E5prime = (315*(123 + 143*x*x*(11*x*x-9)))/static_cast<Real>(12584);
|
||||
BOOST_CHECK_CLOSE_FRACTION(ls5.prime(x), E5prime, 29000*tol);
|
||||
x += 1/static_cast<Real>(1 << 9);
|
||||
}
|
||||
|
||||
// Test norm:
|
||||
// E_1 = x
|
||||
Real expected_norm_sq = 2*third<Real>();
|
||||
BOOST_CHECK_CLOSE_FRACTION(expected_norm_sq, ls1.norm_sq(), tol);
|
||||
|
||||
// E_2 = P[sub 2](x) - 2P[sup 0](x)/5
|
||||
expected_norm_sq = 2/static_cast<Real>(5) + 8/static_cast<Real>(25);
|
||||
BOOST_CHECK_CLOSE_FRACTION(expected_norm_sq, ls2.norm_sq(), tol);
|
||||
|
||||
// E_3 = P[sub 3](x) - 9P[sub 1]/14
|
||||
expected_norm_sq = 2/static_cast<Real>(7) + 9*9*2*third<Real>()/static_cast<Real>(14*14);
|
||||
BOOST_CHECK_CLOSE_FRACTION(expected_norm_sq, ls3.norm_sq(), tol);
|
||||
|
||||
// E_4 = P[sub 4](x) -20P[sub 2](x)/27 + 14P[sub 0](x)/891
|
||||
expected_norm_sq = static_cast<Real>(2)/static_cast<Real>(9) + static_cast<Real>(20*20*2)/static_cast<Real>(27*27*5) + 14*14*2/static_cast<Real>(891*891);
|
||||
BOOST_CHECK_CLOSE_FRACTION(expected_norm_sq, ls4.norm_sq(), tol);
|
||||
|
||||
// E_5 = P[sub 5](x) - 35P[sub 3](x)/44 + 135P[sub 1](x)/12584
|
||||
expected_norm_sq = 2/static_cast<Real>(11) + (35*35/static_cast<Real>(44*44))*(2/static_cast<Real>(7)) + (135*135/static_cast<Real>(12584*12584))*2*third<Real>();
|
||||
BOOST_CHECK_CLOSE_FRACTION(expected_norm_sq, ls5.norm_sq(), tol);
|
||||
|
||||
// Only zero of E1 is 0:
|
||||
std::vector<Real> zeros = ls1.zeros();
|
||||
BOOST_CHECK(zeros.size() == 1);
|
||||
BOOST_CHECK_SMALL(zeros[0], tol);
|
||||
BOOST_CHECK_SMALL(ls1(zeros[0]), tol);
|
||||
|
||||
zeros = ls2.zeros();
|
||||
BOOST_CHECK(zeros.size() == 1);
|
||||
BOOST_CHECK_CLOSE_FRACTION(zeros[0], sqrt(3/static_cast<Real>(5)), tol);
|
||||
BOOST_CHECK_SMALL(ls2(zeros[0]), tol);
|
||||
|
||||
zeros = ls3.zeros();
|
||||
BOOST_CHECK(zeros.size() == 2);
|
||||
BOOST_CHECK_SMALL(zeros[0], tol);
|
||||
BOOST_CHECK_CLOSE_FRACTION(zeros[1], sqrt(6/static_cast<Real>(7)), tol);
|
||||
|
||||
|
||||
zeros = ls4.zeros();
|
||||
BOOST_CHECK(zeros.size() == 2);
|
||||
Real expected = sqrt( (55 - 2*sqrt(static_cast<Real>(330)))/static_cast<Real>(11) )/static_cast<Real>(3);
|
||||
BOOST_CHECK_CLOSE_FRACTION(zeros[0], expected, tol);
|
||||
|
||||
expected = sqrt( (55 + 2*sqrt(static_cast<Real>(330)))/static_cast<Real>(11) )/static_cast<Real>(3);
|
||||
BOOST_CHECK_CLOSE_FRACTION(zeros[1], expected, 10*tol);
|
||||
|
||||
|
||||
zeros = ls5.zeros();
|
||||
BOOST_CHECK(zeros.size() == 3);
|
||||
BOOST_CHECK_SMALL(zeros[0], tol);
|
||||
|
||||
expected = sqrt( ( 195 - sqrt(static_cast<Real>(6045)) )/static_cast<Real>(286));
|
||||
BOOST_CHECK_CLOSE_FRACTION(zeros[1], expected, tol);
|
||||
|
||||
expected = sqrt( ( 195 + sqrt(static_cast<Real>(6045)) )/static_cast<Real>(286));
|
||||
BOOST_CHECK_CLOSE_FRACTION(zeros[2], expected, tol);
|
||||
|
||||
|
||||
for (size_t i = 6; i < 50; ++i)
|
||||
{
|
||||
legendre_stieltjes<Real> En(i);
|
||||
zeros = En.zeros();
|
||||
for(auto const & zero : zeros)
|
||||
{
|
||||
BOOST_CHECK_SMALL(En(zero), 50*tol);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
BOOST_AUTO_TEST_CASE(LegendreStieltjesZeros)
|
||||
{
|
||||
test_legendre_stieltjes<double>();
|
||||
test_legendre_stieltjes<long double>();
|
||||
test_legendre_stieltjes<cpp_bin_float_quad>();
|
||||
//test_legendre_stieltjes<boost::multiprecision::cpp_bin_float_100>();
|
||||
}
|
||||
Reference in New Issue
Block a user