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

Updated rational function evaluation with compile time detection of array size.

Browse Source 
Added tests.
Updated Jamfile.


[SVN r3214]
This commit is contained in:
John Maddock
2006-09-22 18:03:57 +00:00
parent 153e797b74
commit 4679d4e421
5 changed files with 2896 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files

89
doc/rational.qbk
View File

@@ -8,29 +8,74 @@
#include <boost/math/tools/rational.hpp>
``
// Polynomials:
template <std::size_t N, class T, class V>
V evaluate_polynomial(const T(&poly)[N], const V& val);
template <std::size_t N, class T, class V>
V evaluate_polynomial(const boost::array<T,N>& poly, const V& val);
template <class T, class U>
U evaluate_polynomial(const T* poly, U z, std::size_t count);
// Even polynomials:
template <std::size_t N, class T, class V>
V evaluate_even_polynomial(const T(&poly)[N], const V& z);
template <std::size_t N, class T, class V>
V evaluate_even_polynomial(const boost::array<T,N>& poly, const V& z);
template <class T, class U>
U evaluate_even_polynomial(const T* poly, U z, std::size_t count);
// Odd polynomials
template <std::size_t N, class T, class V>
V evaluate_odd_polynomial(const T(&a)[N], const V& z);
template <std::size_t N, class T, class V>
V evaluate_odd_polynomial(const boost::array<T,N>& a, const V& z);
template <class T, class U>
U evaluate_odd_polynomial(const T* poly, U z, std::size_t count);
// Rational Functions:
template <std::size_t N, class T, class V>
V evaluate_rational(const T(&a)[N], const T(&b)[N], const V& z);
template <std::size_t N, class T, class V>
V evaluate_rational(const boost::array<T,N>& a, const boost::array<T,N>& b, const V& z);
template <class T, class U, class V>
V evaluate_rational(const T* num, const U* denom, V z, unsigned count);
[h4 Description]
Each of the functions come in three variants: a pair of overloaded functions
where the order of the polynomial or rational function is evaluated at
compile time, and an overload that accepts a runtime variable for the size
of the coefficient array. Generally speaking, compile time evaluation of the
array size results in better type safety, is less prone to programmer errors,
and may result in better optimised code. The polynomial evaluation functions
in particular, are specialised for various array sizes, allowing for
loop unrolling, and one hopes, optimal inline expansion.
template <std::size_t N, class T, class V>
V evaluate_polynomial(const T(&poly)[N], const V& val);
template <std::size_t N, class T, class V>
V evaluate_polynomial(const boost::array<T,N>& poly, const V& val);
template <class T, class U>
U evaluate_polynomial(const T* poly, U z, std::size_t count);
Evaluates the [@http://en.wikipedia.org/wiki/Polynomial polynomial] described by
the coefficients stored in /poly/.
The polynomial most have order /count-1/ with /count/ coefficients.
If the size of the array is specified at runtime, then the polynomial
most have order /count-1/ with /count/ coefficients. Otherwise it has
order /N-1/ with /N/ coefficients.
Coefficients should be stored such that the coefficients for the x^i terms
Coefficients should be stored such that the coefficients for the x[super i ] terms
are in poly[i].
The types of the coefficients and of variable
@@ -38,6 +83,12 @@ The types of the coefficients and of variable
This allows, for example, for the coefficient table
to be a table of integers if this is appropriate.
template <std::size_t N, class T, class V>
V evaluate_even_polynomial(const T(&poly)[N], const V& z);
template <std::size_t N, class T, class V>
V evaluate_even_polynomial(const boost::array<T,N>& poly, const V& z);
template <class T, class U>
U evaluate_even_polynomial(const T* poly, U z, std::size_t count);
@@ -45,23 +96,37 @@ As above, but evaluates an even polynomial: one where all the powers
of /z/ are even numbers. Equivalent to calling
`evaluate_polynomial(poly, z*z, count)`.
template <std::size_t N, class T, class V>
V evaluate_odd_polynomial(const T(&a)[N], const V& z);
template <std::size_t N, class T, class V>
V evaluate_odd_polynomial(const boost::array<T,N>& a, const V& z);
template <class T, class U>
U evaluate_odd_polynomial(const T* poly, U z, std::size_t count);
As above but evaluates a polynomial where all the powers are odd numbers.
Equivalent to `evaluate_polynomial(poly+1, z*z, count-1) * z + poly[0]`.
template <std::size_t N, class T, class V>
V evaluate_rational(const T(&num)[N], const T(&denom)[N], const V& z);
template <std::size_t N, class T, class V>
V evaluate_rational(const boost::array<T,N>& num, const boost::array<T,N>& denom, const V& z);
template <class T, class U, class V>
V evaluate_rational(const T* num, const U* denom, V z, unsigned count);
Evaluates the rational function (the ratio of two polynomials) described by
the coefficients stored in /num/ and /demom/.
Both polynomials most have order /count-1/ with /count/ coefficients.
If the size of the array is specified at runtime then both
polynomials most have order /count-1/ with /count/ coefficients.
Otherwise both polynomials have order /N-1/ with /N/ coefficients.
Array /num/ describes the numerator, and /demon/ the denominator.
Coefficients should be stored such that the coefficients for the x^i terms
Coefficients should be stored such that the coefficients for the x[super i ] terms
are in num[i] and denom[i].
The types of the coefficients and of variable
@@ -69,9 +134,23 @@ The types of the coefficients and of variable
This allows, for example, for one or both of the coefficient tables
to be a table of integers if this is appropriate.
These functions are designed to safely evaluate the result, even when the value
/z/ is very large. As such they do not take advantage of compile time array
sizes to make any optimisations. These functions are best reserved for situations
where /z/ may be large: if you can be sure that numerical overflow will not occur
then polynomial evaluation with compile-time array sizes may offer slightly
better performance.
[h4 Implementation]
Evaluation is by [@http://en.wikipedia.org/wiki/Horner_algorithm Horners method]:
Polynomials are evaluated by
[@http://en.wikipedia.org/wiki/Horner_algorithm Horners method].
If the array size is known at
compile time then the functions dispatch to size-specific implementations
that unroll the evaluation loop.
Rational evaluation is by
[@http://en.wikipedia.org/wiki/Horner_algorithm Horners method]:
with the two polynomials being evaluated
in parallel to make the most of the processors floating-point pipeline.
If /v/ is greater than one, then the polynomials are evaluated in reverse

99
include/boost/math/tools/rational.hpp
View File

@@ -6,14 +6,32 @@
#ifndef BOOST_MATH_TOOLS_RATIONAL_HPP
#define BOOST_MATH_TOOLS_RATIONAL_HPP
#include <boost/array.hpp>
namespace boost{ namespace math{ namespace tools{
//
// Forward declaration to keep two phase lookup happy:
//
template <class T, class U>
U evaluate_polynomial(const T* poly, U const& z, std::size_t count);
namespace detail{
//
// These inline functions evaluate polynomials whose size is
// known at compile time - there's no need for a for-loop
// just an inline application of Horners rule.
//
template <class T, class V>
inline V evaluate_polynomial_c_imp(const T* a, const V&, const mpl::int_<0>*)
{
return static_cast<V>(0);
}
template <class T, class V>
inline V evaluate_polynomial_c_imp(const T* a, const V&, const mpl::int_<1>*)
{
return a[0];
return static_cast<V>(a[0]);
}
template <class T, class V>
@@ -72,16 +90,15 @@ inline V evaluate_polynomial_c_imp(const T* a, const V& val, const Tag*)
} // namespace detail
template <std::size_t N, class T, class V>
inline V evaluate_polynomial(const T(&a)[N], const V& val)
{
typedef mpl::int_<N> tag_type;
return detail::evaluate_polynomial_c_imp(static_cast<const T*>(a), val, static_cast<tag_type const*>(0));
}
//
// Polynomial evaluation with runtime size.
// This requires a for-loop which may be more expensive than
// the loop expanded versions above:
//
template <class T, class U>
U evaluate_polynomial(const T* poly, U const& z, std::size_t count)
{
BOOST_ASSERT(count > 0);
U sum = static_cast<U>(poly[count - 1]);
for(int i = static_cast<int>(count) - 2; i >= 0; --i)
{
@@ -90,19 +107,73 @@ U evaluate_polynomial(const T* poly, U const& z, std::size_t count)
}
return sum;
}
//
// Compile time sized polynomials, just inline forwarders to the
// implementations above:
//
template <std::size_t N, class T, class V>
inline V evaluate_polynomial(const T(&a)[N], const V& val)
{
typedef mpl::int_<N> tag_type;
return detail::evaluate_polynomial_c_imp(static_cast<const T*>(a), val, static_cast<tag_type const*>(0));
}
template <std::size_t N, class T, class V>
inline V evaluate_polynomial(const boost::array<T,N>& a, const V& val)
{
typedef mpl::int_<N> tag_type;
return detail::evaluate_polynomial_c_imp(static_cast<const T*>(a.data()), val, static_cast<tag_type const*>(0));
}
//
// Even polynomials are trivial: just square the argument!
//
template <class T, class U>
inline U evaluate_even_polynomial(const T* poly, U z, std::size_t count)
{
return evaluate_polynomial(poly, z*z, count);
}
template <std::size_t N, class T, class V>
inline V evaluate_even_polynomial(const T(&a)[N], const V& z)
{
return evaluate_polynomial(a, z*z);
}
template <std::size_t N, class T, class V>
inline V evaluate_even_polynomial(const boost::array<T,N>& a, const V& z)
{
return evaluate_polynomial(a, z*z);
}
//
// Odd polynomials come next:
//
template <class T, class U>
inline U evaluate_odd_polynomial(const T* poly, U z, std::size_t count)
{
return poly[0] + z * evaluate_polynomial(poly+1, z*z, count-1);
}
template <std::size_t N, class T, class V>
inline V evaluate_odd_polynomial(const T(&a)[N], const V& z)
{
typedef mpl::int_<N-1> tag_type;
return a[0] + z * detail::evaluate_polynomial_c_imp(static_cast<const T*>(a) + 1, z*z, static_cast<tag_type const*>(0));
}
template <std::size_t N, class T, class V>
inline V evaluate_odd_polynomial(const boost::array<T,N>& a, const V& z)
{
typedef mpl::int_<N-1> tag_type;
return a[0] + z * detail::evaluate_polynomial_c_imp(static_cast<const T*>(a.data()) + 1, z*z, static_cast<tag_type const*>(0));
}
//
// Rational functions: numerator and denominator must be
// equal in size. These always have a for-loop and so may be less
// efficient than evaluating a pair of polynomials. However, there
// are some tricks we can use to prevent overflow that might otherwise
// occur in polynomial evaluation, if z is large. This is important
// in our Lanczos code for example.
//
template <class T, class U, class V>
V evaluate_rational(const T* num, const U* denom, const V& z_, std::size_t count)
{
@@ -136,6 +207,18 @@ V evaluate_rational(const T* num, const U* denom, const V& z_, std::size_t count
return s1 / s2;
}
template <std::size_t N, class T, class V>
inline V evaluate_rational(const T(&a)[N], const T(&b)[N], const V& z)
{
return evaluate_rational(a, b, z, N);
}
template <std::size_t N, class T, class V>
inline V evaluate_rational(const boost::array<T,N>& a, const boost::array<T,N>& b, const V& z)
{
return evaluate_rational(a.data(), b.data(), z, N);
}
} // namespace tools
} // namespace math
} // namespace boost

2
test/Jamfile.v2
View File

@@ -32,10 +32,10 @@ run test_igamma_inva.cpp ;
run test_minima.cpp ;
run test_normal.cpp ;
run test_promotion.cpp ;
run test_rationals.cpp ;
run test_remez.cpp ;
run test_roots.cpp ;
run test_students_t.cpp ;
run test_tgamma_ratio.cpp ;
run test_toms748_solve.cpp ;

2358
test/test_rationals.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

362
tools/generate_rational_test.cpp Normal file
View File

@@ -0,0 +1,362 @@
#include <boost/math/tools/ntl.hpp>
#include <boost/tr1/random.hpp>
#include <boost/math/tools/rational.hpp>
#include <iostream>
#include <fstream>
int main()
{
using namespace boost::math;
using namespace boost::math::tools;
NTL::RR::SetPrecision(500);
NTL::RR::SetOutputPrecision(40);
std::tr1::mt19937 rnd;
std::tr1::variate_generator<
std::tr1::mt19937,
std::tr1::uniform_int<> > gen(rnd, std::tr1::uniform_int<>(1, 12));
for(unsigned i = 1; i < 12; ++i)
{
std::vector<int> coef;
for(unsigned j = 0; j < i; ++j)
{
coef.push_back(gen());
}
std::cout <<
" //\n"
" // Polynomials of order " << i-1 << "\n"
" //\n"
" static const U n" << i << "c[" << i << "] = { ";
for(unsigned j = 0; j < i; ++j)
{
if(j)
std::cout << ", ";
std::cout << coef[j];
}
std::cout << " };\n";
std::cout <<
" static const boost::array<U, " << i << "> n" << i << "a = { ";
for(unsigned j = 0; j < i; ++j)
{
if(j)
std::cout << ", ";
std::cout << coef[j];
}
std::cout << " };\n";
NTL::RR r1 = boost::math::tools::evaluate_polynomial(&coef[0], NTL::RR(0.125), i);
NTL::RR r2 = boost::math::tools::evaluate_polynomial(&coef[0], NTL::RR(0.25), i);
NTL::RR r3 = boost::math::tools::evaluate_polynomial(&coef[0], NTL::RR(0.75), i);
NTL::RR r4 = boost::math::tools::evaluate_polynomial(&coef[0], NTL::RR(1) - NTL::RR(1) / 64, i);
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_polynomial(n" << i << "c, static_cast<T>(0.125), " << i << "),\n"
" static_cast<T>(" << r1 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_polynomial(n" << i << "c, static_cast<T>(0.25), " << i << "),\n"
" static_cast<T>(" << r2 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_polynomial(n" << i << "c, static_cast<T>(0.75), " << i << "),\n"
" static_cast<T>(" << r3 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_polynomial(n" << i << "c, static_cast<T>(1.0f - 1.0f/64.0f), " << i << "),\n"
" static_cast<T>(" << r4 << "L),\n"
" tolerance);\n\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_polynomial(n" << i << "c, static_cast<T>(0.125)),\n"
" static_cast<T>(" << r1 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_polynomial(n" << i << "c, static_cast<T>(0.25)),\n"
" static_cast<T>(" << r2 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_polynomial(n" << i << "c, static_cast<T>(0.75)),\n"
" static_cast<T>(" << r3 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_polynomial(n" << i << "c, static_cast<T>(1.0f - 1.0f/64.0f)),\n"
" static_cast<T>(" << r4 << "L),\n"
" tolerance);\n\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_polynomial(n" << i << "a, static_cast<T>(0.125)),\n"
" static_cast<T>(" << r1 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_polynomial(n" << i << "a, static_cast<T>(0.25)),\n"
" static_cast<T>(" << r2 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_polynomial(n" << i << "a, static_cast<T>(0.75)),\n"
" static_cast<T>(" << r3 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_polynomial(n" << i << "a, static_cast<T>(1.0f - 1.0f/64.0f)),\n"
" static_cast<T>(" << r4 << "L),\n"
" tolerance);\n\n";
r1 = boost::math::tools::evaluate_even_polynomial(&coef[0], NTL::RR(0.125), i);
r2 = boost::math::tools::evaluate_even_polynomial(&coef[0], NTL::RR(0.25), i);
r3 = boost::math::tools::evaluate_even_polynomial(&coef[0], NTL::RR(0.75), i);
r4 = boost::math::tools::evaluate_even_polynomial(&coef[0], NTL::RR(1) - NTL::RR(1) / 64, i);
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_even_polynomial(n" << i << "c, static_cast<T>(0.125), " << i << "),\n"
" static_cast<T>(" << r1 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_even_polynomial(n" << i << "c, static_cast<T>(0.25), " << i << "),\n"
" static_cast<T>(" << r2 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_even_polynomial(n" << i << "c, static_cast<T>(0.75), " << i << "),\n"
" static_cast<T>(" << r3 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_even_polynomial(n" << i << "c, static_cast<T>(1.0f - 1.0f/64.0f), " << i << "),\n"
" static_cast<T>(" << r4 << "L),\n"
" tolerance);\n\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_even_polynomial(n" << i << "c, static_cast<T>(0.125)),\n"
" static_cast<T>(" << r1 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_even_polynomial(n" << i << "c, static_cast<T>(0.25)),\n"
" static_cast<T>(" << r2 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_even_polynomial(n" << i << "c, static_cast<T>(0.75)),\n"
" static_cast<T>(" << r3 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_even_polynomial(n" << i << "c, static_cast<T>(1.0f - 1.0f/64.0f)),\n"
" static_cast<T>(" << r4 << "L),\n"
" tolerance);\n\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_even_polynomial(n" << i << "a, static_cast<T>(0.125)),\n"
" static_cast<T>(" << r1 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_even_polynomial(n" << i << "a, static_cast<T>(0.25)),\n"
" static_cast<T>(" << r2 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_even_polynomial(n" << i << "a, static_cast<T>(0.75)),\n"
" static_cast<T>(" << r3 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_even_polynomial(n" << i << "a, static_cast<T>(1.0f - 1.0f/64.0f)),\n"
" static_cast<T>(" << r4 << "L),\n"
" tolerance);\n\n";
if(i > 1)
{
r1 = boost::math::tools::evaluate_odd_polynomial(&coef[0], NTL::RR(0.125), i);
r2 = boost::math::tools::evaluate_odd_polynomial(&coef[0], NTL::RR(0.25), i);
r3 = boost::math::tools::evaluate_odd_polynomial(&coef[0], NTL::RR(0.75), i);
r4 = boost::math::tools::evaluate_odd_polynomial(&coef[0], NTL::RR(1) - NTL::RR(1) / 64, i);
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_odd_polynomial(n" << i << "c, static_cast<T>(0.125), " << i << "),\n"
" static_cast<T>(" << r1 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_odd_polynomial(n" << i << "c, static_cast<T>(0.25), " << i << "),\n"
" static_cast<T>(" << r2 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_odd_polynomial(n" << i << "c, static_cast<T>(0.75), " << i << "),\n"
" static_cast<T>(" << r3 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_odd_polynomial(n" << i << "c, static_cast<T>(1.0f - 1.0f/64.0f), " << i << "),\n"
" static_cast<T>(" << r4 << "L),\n"
" tolerance);\n\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_odd_polynomial(n" << i << "c, static_cast<T>(0.125)),\n"
" static_cast<T>(" << r1 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_odd_polynomial(n" << i << "c, static_cast<T>(0.25)),\n"
" static_cast<T>(" << r2 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_odd_polynomial(n" << i << "c, static_cast<T>(0.75)),\n"
" static_cast<T>(" << r3 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_odd_polynomial(n" << i << "c, static_cast<T>(1.0f - 1.0f/64.0f)),\n"
" static_cast<T>(" << r4 << "L),\n"
" tolerance);\n\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_odd_polynomial(n" << i << "a, static_cast<T>(0.125)),\n"
" static_cast<T>(" << r1 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_odd_polynomial(n" << i << "a, static_cast<T>(0.25)),\n"
" static_cast<T>(" << r2 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_odd_polynomial(n" << i << "a, static_cast<T>(0.75)),\n"
" static_cast<T>(" << r3 << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_odd_polynomial(n" << i << "a, static_cast<T>(1.0f - 1.0f/64.0f)),\n"
" static_cast<T>(" << r4 << "L),\n"
" tolerance);\n\n";
}
r1 = boost::math::tools::evaluate_polynomial(&coef[0], NTL::RR(0.125), i);
r2 = boost::math::tools::evaluate_polynomial(&coef[0], NTL::RR(0.25), i);
r3 = boost::math::tools::evaluate_polynomial(&coef[0], NTL::RR(0.75), i);
r4 = boost::math::tools::evaluate_polynomial(&coef[0], NTL::RR(1) - NTL::RR(1) / 64, i);
coef.clear();
for(unsigned j = 0; j < i; ++j)
{
coef.push_back(gen());
}
std::cout <<
" //\n"
" // Rational functions of order " << i-1 << "\n"
" //\n"
" static const U d" << i << "c[" << i << "] = { ";
for(unsigned j = 0; j < i; ++j)
{
if(j)
std::cout << ", ";
std::cout << coef[j];
}
std::cout << " };\n";
std::cout <<
" static const boost::array<U, " << i << "> d" << i << "a = { ";
for(unsigned j = 0; j < i; ++j)
{
if(j)
std::cout << ", ";
std::cout << coef[j];
}
std::cout << " };\n";
NTL::RR r1d = boost::math::tools::evaluate_polynomial(&coef[0], NTL::RR(0.125), i);
NTL::RR r2d = boost::math::tools::evaluate_polynomial(&coef[0], NTL::RR(0.25), i);
NTL::RR r3d = boost::math::tools::evaluate_polynomial(&coef[0], NTL::RR(0.75), i);
NTL::RR r4d = boost::math::tools::evaluate_polynomial(&coef[0], NTL::RR(1) - NTL::RR(1) / 64, i);
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_rational(n" << i << "c, d" << i << "c, static_cast<T>(0.125), " << i << "),\n"
" static_cast<T>(" << r1/r1d << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_rational(n" << i << "c, d" << i << "c, static_cast<T>(0.25), " << i << "),\n"
" static_cast<T>(" << r2/r2d << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_rational(n" << i << "c, d" << i << "c, static_cast<T>(0.75), " << i << "),\n"
" static_cast<T>(" << r3/r3d << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_rational(n" << i << "c, d" << i << "c, static_cast<T>(1.0f - 1.0f/64.0f), " << i << "),\n"
" static_cast<T>(" << r4/r4d << "L),\n"
" tolerance);\n\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_rational(n" << i << "c, d" << i << "c, static_cast<T>(0.125)),\n"
" static_cast<T>(" << r1/r1d << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_rational(n" << i << "c, d" << i << "c, static_cast<T>(0.25)),\n"
" static_cast<T>(" << r2/r2d << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_rational(n" << i << "c, d" << i << "c, static_cast<T>(0.75)),\n"
" static_cast<T>(" << r3/r3d << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_rational(n" << i << "c, d" << i << "c, static_cast<T>(1.0f - 1.0f/64.0f)),\n"
" static_cast<T>(" << r4/r4d << "L),\n"
" tolerance);\n\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_rational(n" << i << "a, d" << i << "a, static_cast<T>(0.125)),\n"
" static_cast<T>(" << r1/r1d << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_rational(n" << i << "a, d" << i << "a, static_cast<T>(0.25)),\n"
" static_cast<T>(" << r2/r2d << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_rational(n" << i << "a, d" << i << "a, static_cast<T>(0.75)),\n"
" static_cast<T>(" << r3/r3d << "L),\n"
" tolerance);\n";
std::cout <<
" BOOST_CHECK_CLOSE(\n"
" boost::math::tools::evaluate_rational(n" << i << "a, d" << i << "a, static_cast<T>(1.0f - 1.0f/64.0f)),\n"
" static_cast<T>(" << r4/r4d << "L),\n"
" tolerance);\n\n";
}
return 0;
}