mirror of
https://github.com/boostorg/math.git
synced 2026-01-19 04:22:09 +00:00
Fix minor bug in test_data.hpp.
Updated docs with minima and minimax docs. Also reorganised to separate out the more experimental components. [SVN r3218]
This commit is contained in:
John Maddock
@@ -13,7 +13,7 @@ and forces conversions to RR to proceed via `long double` rather than `double`.
|
||||
The latter change is essential to accurately measure relative errors between
|
||||
high precision calculations (using NTL::RR) and fixed precision calculations
|
||||
(using `long double`). These occur if, for example, you are generating additional
|
||||
__lanczos's.
|
||||
__lanczos's or rational approximations.
|
||||
|
||||
You will also need to include NTL via the header `<boost/math/tools/ntl.hpp>`
|
||||
and include that header /before/ you include any of the special function headers.
|
||||
|
||||
21
doc/internals_overview.qbk
Normal file
21
doc/internals_overview.qbk
Normal file
@@ -0,0 +1,21 @@
|
||||
[section:internals_overview Overview]
|
||||
|
||||
This section contains tools that are used in the development and testing
|
||||
of this library. These tools have only minimal documentation, and crucially
|
||||
['do not have stable interfaces].
|
||||
|
||||
There is no doubt that these components can be improved, but they are also
|
||||
largely incidental to the main purpose of this library.
|
||||
|
||||
These tools are designed to "just get the job done", and receive minimal
|
||||
documentation here, in the hopes that they will help stimulate further
|
||||
submissions to this library.
|
||||
|
||||
[endsect][/section:internals_overview Overview]
|
||||
|
||||
[/
|
||||
Copyright 2006 John Maddock and Paul A. Bristow.
|
||||
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).
|
||||
]
|
||||
@@ -122,10 +122,16 @@ some typical applications in statistics.
|
||||
[include rational.qbk]
|
||||
[include roots.qbk]
|
||||
[include roots_without_derivatives.qbk]
|
||||
[include minima.qbk]
|
||||
[endsect][/section:toolkit Toolkit]
|
||||
|
||||
[section:intenals Tools for Internal Testing and Development Use]
|
||||
[include internals_overview.qbk]
|
||||
[include polynomial.qbk]
|
||||
[include minimax.qbk]
|
||||
[include relative_error.qbk]
|
||||
[include test_data.qbk]
|
||||
[endsect][/section:toolkit Toolkit]
|
||||
[endsect]
|
||||
|
||||
[section:Using_UDT Use with User Defined Floating-Point Types]
|
||||
[include concepts.qbk]
|
||||
|
||||
66
doc/minima.qbk
Normal file
66
doc/minima.qbk
Normal file
@@ -0,0 +1,66 @@
|
||||
[section:minima Locating Function Minima]
|
||||
|
||||
[caution __caution ]
|
||||
|
||||
[h4 synopsis]
|
||||
|
||||
``
|
||||
#include <boost/math/tools/minima.hpp>
|
||||
``
|
||||
|
||||
template <class F, class T>
|
||||
std::pair<T, T> brent_find_minima(F f, T min, T max, int digits);
|
||||
|
||||
template <class F, class T>
|
||||
std::pair<T, T> brent_find_minima(F f, T min, T max, int bits, boost::uintmax_t& max_iter);
|
||||
|
||||
[h4 Description]
|
||||
|
||||
These two functions locate the minima of the continuous function /f/ using Brent's
|
||||
algorithm. Parameters are:
|
||||
|
||||
[variablelist
|
||||
[[f] [The function to minimise. The function should be smooth over the
|
||||
range \[min,max\], with no maxima occuring in that interval.]]
|
||||
[[min] [The lower endpoint of the range in which to search
|
||||
for the minima.]]
|
||||
[[max] [The upper endpoint of the range in which to search
|
||||
for the minima.]]
|
||||
[[bits] [The number of bits precision to which the minima should be found.
|
||||
Note that in principle, the minima can not be located to greater
|
||||
accuracy than the square root of machine epsilon, therefore if /bits/
|
||||
is set to a value greater than one half of the bits in type T, then
|
||||
the value will be ignored.]]
|
||||
[[max_iter] [The maximum number of iterations to use
|
||||
in the algorithm, if not provided the algorithm will just
|
||||
keep on going until the minima is found.]]
|
||||
]
|
||||
|
||||
[*Returns:] a pair containing the value of the absissa at the minima and the value
|
||||
of f(x) at the minima.
|
||||
|
||||
[h4 Implementation]
|
||||
|
||||
This is a reasonably faithful implementation of Brent's algorithm, refer
|
||||
to:
|
||||
|
||||
Brent, R.P. 1973, Algorithms for Minimization without Derivatives
|
||||
(Englewood Cliffs, NJ: Prentice-Hall), Chapter 5.
|
||||
|
||||
Numerical Recipes in C, The Art of Scientific Computing,
|
||||
Second Edition, William H. Press, Saul A. Teukolsky,
|
||||
William T. Vetterling, and Brian P. Flannery.
|
||||
Cambridge University Press. 1988, 1992.
|
||||
|
||||
An algorithm with guarenteed convergence for finding a zero
|
||||
of a function, R. P. Brent, The Computer Jounal, Vol 44, 1971.
|
||||
|
||||
[endsect][/section:minima Locating Function Minima]
|
||||
|
||||
[/
|
||||
Copyright 2006 John Maddock and Paul A. Bristow.
|
||||
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).
|
||||
]
|
||||
|
||||
149
doc/minimax.qbk
Normal file
149
doc/minimax.qbk
Normal file
@@ -0,0 +1,149 @@
|
||||
[section:minimax Minimax Approximations and the Remez Algorithm]
|
||||
|
||||
The directory libs/math/minimax contains a command line driven
|
||||
program for the generation of mininmax approximations using the Remez
|
||||
algorithm. Both polynomial and rational approximations are supported,
|
||||
although the latter are tricky to converge, and in addition this
|
||||
implementation only converges correctly when the rational function is
|
||||
optimised over a positive domain. It's not clear if this is a limitation
|
||||
of the algorithm or this implementation. No such limitations are present
|
||||
for polynomial approximations which should always converge smoothly.
|
||||
|
||||
It's worth stressing that developing rational approximations to functions
|
||||
is often not an easy task, and one to which many books have been devoted.
|
||||
To use this tool, you will need to have a reasonable grasp of what the Remez
|
||||
algorithm is, and the general form of the approximation you want to achieve.
|
||||
|
||||
The program consists of two parts:
|
||||
|
||||
[variablelist
|
||||
[[main.cpp][Contains the command line parser, and all the calls to the Remez code.]]
|
||||
[[f.cpp][Contains the function to appoximate.]]
|
||||
]
|
||||
|
||||
Therefore to use this tool, you must modify f.cpp to return the function to
|
||||
approximate. The tools supports multiple function approximations within
|
||||
the same compiled program: each as a separate variant:
|
||||
|
||||
NTL::RR f(const NTL::RR& x, int variant);
|
||||
|
||||
Returns the value of the function /variant/ at point /x/. So if you
|
||||
wish you can just add the function to approximate as a new variant
|
||||
after the existing examples.
|
||||
|
||||
In addition to those two files, the program needs to be linked to
|
||||
a [link math_toolkit.Using_UDT.use_NTL patched NTL library to compile].
|
||||
|
||||
Note that the function /f/ must return the rational part of the
|
||||
approximation: for example if you are approximating a function
|
||||
/f(x)/ then it is quite common to use:
|
||||
|
||||
f(x) = g(x)(Y + R(x))
|
||||
|
||||
where /g(x)/ is the dominant part of /f(x)/, /Y/ is some constant, and
|
||||
/R(x)/ is the rational approximation part, usually optimised for a low
|
||||
absolute error compared to |Y|.
|
||||
|
||||
In this case you would define /f/ to return ['f(x)/g(x)] and then set the
|
||||
y-offset of the approximation to /Y/ (see command line options below).
|
||||
|
||||
Many other forms are possible, but in all cases the objective is to
|
||||
split /f(x)/ into a dominant part that you can evaluate easily using
|
||||
standard math functions, and a smooth and slowly changing rational approximation
|
||||
part. Refer to your favorite textbook for more examples.
|
||||
|
||||
Command line options for the program are as follows:
|
||||
|
||||
[variablelist
|
||||
[[variant N][Sets the current function variant to N. This allows multiple functions
|
||||
that are to be approximated to be compiled into the same executable.
|
||||
Defaults to 0.]]
|
||||
[[range a b][Sets the domain for the approximation to the range \[a,b\], defaults
|
||||
to \[0,1\].]]
|
||||
[[relative][Sets the Remez code to optimise for relative error. This is the default
|
||||
at program startup. Note that relative error can only be used
|
||||
if f(x) has no roots over the range being optimised.]]
|
||||
[[absolute][Sets the Remez code to optimise for absolute error.]]
|
||||
[[pin \[true|false\]]["Pins" the code so that the rational approximation
|
||||
passes through the origin. Obviously only set this to
|
||||
/true/ if R(0) must be zero. This is typically used when
|
||||
trying to preserve a root at \[0,0\] while also optimising
|
||||
for relative error.]]
|
||||
[[order N D][Sets the order of the approximation to /N/ in the numerator and /D/
|
||||
in the denominator. If /D/ is zero then the result will be a polynomial
|
||||
approximation. There will be N+D+2 coefficients in total, the first
|
||||
coefficient of the numerator is zero if /pin/ was set to true, and the
|
||||
first coefficient of the denominator is always one.]]
|
||||
[[working-precision N][Sets the working precision of NTL::RR to /N/ binary digits. Defaults to 250.]]
|
||||
[[target-precision N][Sets the precision of printed output to /N/ binary digits:
|
||||
set to the same number of digits as the type that will be used to
|
||||
evaluate the approximation. Defaults to 53 (for double precision).]]
|
||||
[[skew val]["Skews" the initial interpolated control points towards one
|
||||
end or the other of the range. Positive values skew the
|
||||
initial control points towards the left hand side of the
|
||||
range, and negative values towards the right hand side.
|
||||
If an approximation won't converge (a common situation)
|
||||
try adjusting the skew parameter until the first step yields
|
||||
the smallest possible error. /val/ should be in the range
|
||||
\[-100,+100\], the default is zero.]]
|
||||
[[brake val][Sets a brake on each step so that the change in the
|
||||
control points is braked by /val%/. Defaults to 50,
|
||||
try a higher value if an approximation won't converge,
|
||||
or a lower value to get speedier convergence.]]
|
||||
[[x-offset val][Sets the x-offset to /val/: the approximation will
|
||||
be generated for `f(x + X) + Y` where /X/ is the x-offset
|
||||
and /Y/ is the y-offset. Defaults to zero. To avoid
|
||||
rounding errors, take care to specify a value that can
|
||||
be exactly represented as a floating point number.]]
|
||||
[[y-offset val][Sets the y-offset to /val/: the approximation will
|
||||
be generated for `f(x + X) + Y` where /X/ is the x-offset
|
||||
and /Y/ is the y-offset. Defaults to zero. To avoid
|
||||
rounding errors, take care to specify a value that can
|
||||
be exactly represented as a floating point number.]]
|
||||
[[y-offset auto][Sets the y-offset to the average value of f(x)
|
||||
evaluated at the two endpoints of the range plus the midpoint
|
||||
of the range. The calculated value is deliberately truncated
|
||||
to /float/ precision (and should be stored as a /float/
|
||||
in your code). The approximation will
|
||||
be generated for `f(x + X) + Y` where /X/ is the x-offset
|
||||
and /Y/ is the y-offset. Defaults to zero.]]
|
||||
[[graph N][Prints N evaluations of f(x) at evenly spaced points over the
|
||||
range being optimised. If unspecified then /N/ defaults
|
||||
to 3. Use to check that f(x) is indeed smooth over the range
|
||||
of interest.]]
|
||||
[[step N][Performs /N/ steps, or one step if /N/ is unspecified.
|
||||
After each step prints: the peek error at the extrema of
|
||||
the error function of the approximation,
|
||||
the theoretical error term solved for on the last step,
|
||||
and the maximum relative change in the location of the
|
||||
Chebeshev control points. The approximation is converged on the
|
||||
minimax solution when the two error terms are (approximately)
|
||||
equal, and the change in the control points has decreased to
|
||||
a suitably small value.]]
|
||||
[[test \[float|double|long\]][Tests the current approximation at float,
|
||||
double, or long double precision. Useful to check for rounding
|
||||
errors in evaluating the approximation at fixed precision.
|
||||
Tests are conducted at the extrema of the error function of the
|
||||
approximation, and at the zeros of the error function.]]
|
||||
[[rescale a b][Takes the current Chebeshev control points, and rescales them
|
||||
over a new interval \[a,b\]. Sometimes this can be used to obtain
|
||||
starting control points for an approximation that can not otherwise be
|
||||
converged.]]
|
||||
[[rotate][Moves one term from the numerator to the denominator, but keeps the
|
||||
Chebeshev control points the same. Sometimes this can be used to obtain
|
||||
starting control points for an approximation that can not otherwise be
|
||||
converged.]]
|
||||
[[info][Prints out the current approximation: the location of the zeros of the
|
||||
error function, the location of the Chebeshev control points, the
|
||||
x and y offsets, and of course the coefficients of the polynomials.]]
|
||||
]
|
||||
|
||||
|
||||
[endsect][/section:minimax Minimax Approximations and the Remez Algorithm]
|
||||
|
||||
[/
|
||||
Copyright 2006 John Maddock and Paul A. Bristow.
|
||||
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).
|
||||
]
|
||||
@@ -37,8 +37,24 @@ This function is used for testing a function against tabulated test data.
|
||||
|
||||
The return type contains statistical data on the relative errors (max, mean,
|
||||
variance, and the number of test cases etc), as well as the row of test data that
|
||||
caused the largest relative error. Type test_result is work in progress, refer
|
||||
to the header for more details.
|
||||
caused the largest relative error. Public members of type test_result are:
|
||||
|
||||
[variablelist
|
||||
[[`unsigned worst()const;`][Returns the row at which the worst error occured.]]
|
||||
[[`T min()const;`][Returns the smallest relative error found.]]
|
||||
[[`T max()const;`][Returns the largest relative error found.]]
|
||||
[[`T mean()const;`][Returns the mean error found.]]
|
||||
[[`boost::uintmax_t count()const;`][Returns the number of test cases.]]
|
||||
[[`T variance()const;`][Returns the variance of the errors found.]]
|
||||
[[`T variance1()const;`][Returns the unbiased variance of the errors found.]]
|
||||
[[`T rms()const`][Returns the Root Mean Square, or quadratic mean of the errors.]]
|
||||
[[`test_result& operator+=(const test_result& t)`][Combines two test_result's into
|
||||
a single result.]]
|
||||
]
|
||||
|
||||
The template parameter of test_result, is the same type as the values in the two
|
||||
dimentional array passed to function /test/, roughly that's
|
||||
`A::value_type::value_type`.
|
||||
|
||||
Parameter /a/ is a matrix of test data: and must be a standard library Sequence type,
|
||||
that contains another Sequence type:
|
||||
@@ -56,8 +72,10 @@ from a row of test data in /a/. Typically type F2 is created with Boost.Lambda:
|
||||
the example below.
|
||||
|
||||
If the function under test returns a non-finite value when a finite result is
|
||||
expected, or if a gross error is found, then a message is sent to `std::cerr`.
|
||||
This is mainly a debugging/development aid (and a good place for a breakpoint).
|
||||
expected, or if a gross error is found, then a message is sent to `std::cerr`,
|
||||
and a call to BOOST_ERROR() made (which means that including this header requires
|
||||
you use Boost.Test). This is mainly a debugging/development aid
|
||||
(and a good place for a breakpoint).
|
||||
|
||||
[h4 Example]
|
||||
|
||||
|
||||
@@ -66,24 +66,24 @@ graphing program (or spreadsheet) for further analysis.
|
||||
template <class F>
|
||||
test_data& insert(F func, const parameter_info<T>& arg1, const parameter_info<T>& arg2, const parameter_info<T>& arg3);
|
||||
|
||||
void clear(){ m_data.clear(); }
|
||||
void clear();
|
||||
|
||||
// access:
|
||||
iterator begin() { return m_data.begin(); }
|
||||
iterator end() { return m_data.end(); }
|
||||
const_iterator begin()const { return m_data.begin(); }
|
||||
const_iterator end()const { return m_data.end(); }
|
||||
bool operator==(const test_data& d)const{ return m_data == d.m_data; }
|
||||
bool operator!=(const test_data& d)const{ return m_data != d.m_data; }
|
||||
void swap(test_data& other){ m_data.swap(other.m_data); }
|
||||
size_type size()const{ return m_data.size(); }
|
||||
size_type max_size()const{ return m_data.max_size(); }
|
||||
bool empty()const{ return m_data.empty(); }
|
||||
iterator begin();
|
||||
iterator end();
|
||||
const_iterator begin()const;
|
||||
const_iterator end()const;
|
||||
bool operator==(const test_data& d)const;
|
||||
bool operator!=(const test_data& d)const;
|
||||
void swap(test_data& other);
|
||||
size_type size()const;
|
||||
size_type max_size()const;
|
||||
bool empty()const;
|
||||
|
||||
bool operator < (const test_data& dat)const{ return m_data < dat.m_data; }
|
||||
bool operator <= (const test_data& dat)const{ return m_data <= dat.m_data; }
|
||||
bool operator > (const test_data& dat)const{ return m_data > dat.m_data; }
|
||||
bool operator >= (const test_data& dat)const{ return m_data >= dat.m_data; }
|
||||
bool operator < (const test_data& dat)const;
|
||||
bool operator <= (const test_data& dat)const;
|
||||
bool operator > (const test_data& dat)const;
|
||||
bool operator >= (const test_data& dat)const;
|
||||
};
|
||||
|
||||
template <class charT, class traits, class T>
|
||||
@@ -136,7 +136,7 @@ one could do this like so:
|
||||
|
||||
// insert 500 points at uniform intervals between just after -3 and 100:
|
||||
double (*pf)(double) = boost::math::lgamma;
|
||||
data.insert(pf, make_periodic_param(0.00001, 100.0, 500));
|
||||
data.insert(pf, make_periodic_param(-3.0 + 0.00001, 100.0, 500));
|
||||
|
||||
// print out in csv format:
|
||||
write_csv(std::cout, data, ", ");
|
||||
@@ -288,7 +288,9 @@ when the table is large.
|
||||
|
||||
Alternatively, lets say we want to profile a continued fraction for
|
||||
convergence and error. As an example, we'll use the continued fraction
|
||||
for the upper incomplete gamma function:
|
||||
for the upper incomplete gamma function, the following function object
|
||||
returns the next a[sub N ] and b[sub N ] of the continued fraction
|
||||
each time it's invoked:
|
||||
|
||||
template <class T>
|
||||
struct upper_incomplete_gamma_fract
|
||||
@@ -383,7 +385,42 @@ of a and z.
|
||||
[h4 reference]
|
||||
[#test_data_reference]
|
||||
|
||||
TODO...
|
||||
Most of this tool has been descibed already in the examples above, we'll
|
||||
just add the following notes on the non-member functions:
|
||||
|
||||
template <class T>
|
||||
parameter_info<T> make_random_param(T start_range, T end_range, int n_points);
|
||||
|
||||
Tells class test_data to test /n_points/ random values in the range
|
||||
[start_range,end_range].
|
||||
|
||||
template <class T>
|
||||
parameter_info<T> make_periodic_param(T start_range, T end_range, int n_points);
|
||||
|
||||
Tells class test_data to test /n_points/ evenly spaced values in the range
|
||||
[start_range,end_range].
|
||||
|
||||
template <class T>
|
||||
parameter_info<T> make_power_param(T basis, int start_exponent, int end_exponent);
|
||||
|
||||
Tells class test_data to test points of the form [^/basis + 2[super expon]/] for each
|
||||
/expon/ in the range [start_exponent, end_exponent].
|
||||
|
||||
template <class T>
|
||||
bool get_user_parameter_info(parameter_info<T>& info, const char* param_name);
|
||||
|
||||
Prompts the user for the parameter range and form to use.
|
||||
|
||||
Finally, if we don't want the parameter to be included in the output, we can tell
|
||||
test_data by setting it a "dummy paramter":
|
||||
|
||||
parameter_info<double> p = make_random_param(2.0, 5.0, 10);
|
||||
p.type |= dummy_param;
|
||||
|
||||
This is useful when the functor used transforms the parameter in some way
|
||||
before passing it to the function under test, usually the functor will then
|
||||
return both the transformed input and the result in a tuple, so there's no
|
||||
need for the original pseudo-parameter to be included in program output.
|
||||
|
||||
[endsect][/section:test_data Graphing, Profiling, and Generating Test Data for Special Functions]
|
||||
|
||||
|
||||
@@ -379,7 +379,7 @@ void test_data<T>::create_test_points(std::set<T>& points, const parameter_info<
|
||||
{
|
||||
random_type r = gen();
|
||||
power_type p = ldexp(static_cast<power_type>(r), power);
|
||||
points.insert(truncate_to_float(real_cast<float>(p)));
|
||||
points.insert(truncate_to_float(real_cast<float>(arg1.z1 + p)));
|
||||
}
|
||||
}
|
||||
break;
|
||||
|
||||
@@ -10,10 +10,9 @@
|
||||
|
||||
#include <cmath>
|
||||
|
||||
template <class T>
|
||||
T f0(T x, int variant)
|
||||
NTL::RR f(const NTL::RR& x, int variant)
|
||||
{
|
||||
static const T tiny = boost::math::tools::min_value<T>() * 64;
|
||||
static const NTL::RR tiny = boost::math::tools::min_value<NTL::RR>() * 64;
|
||||
switch(variant)
|
||||
{
|
||||
case 0:
|
||||
@@ -23,9 +22,12 @@ T f0(T x, int variant)
|
||||
case 2:
|
||||
return boost::math::erf(x) / x - 1.125;
|
||||
case 3:
|
||||
if(x == 0)
|
||||
x += tiny;
|
||||
return boost::math::lgamma(x+2) / x - 0.5;
|
||||
{
|
||||
NTL::RR y(x);
|
||||
if(y == 0)
|
||||
y += tiny;
|
||||
return boost::math::lgamma(y+2) / y - 0.5;
|
||||
}
|
||||
case 4:
|
||||
//
|
||||
// lgamma in the range [2,3], use:
|
||||
@@ -37,7 +39,7 @@ T f0(T x, int variant)
|
||||
//
|
||||
if(x == 0)
|
||||
{
|
||||
return boost::lexical_cast<T>("0.42278433509846713939348790991759756895784066406008") / 3;
|
||||
return boost::lexical_cast<NTL::RR>("0.42278433509846713939348790991759756895784066406008") / 3;
|
||||
}
|
||||
return boost::math::lgamma(x+2) / (x * (x+3));
|
||||
case 5:
|
||||
@@ -49,8 +51,8 @@ T f0(T x, int variant)
|
||||
//
|
||||
// works well over [1, 1.5] but not near 2 :-(
|
||||
//
|
||||
NTL::RR r1 = boost::lexical_cast<T>("0.57721566490153286060651209008240243104215933593992");
|
||||
NTL::RR r2 = boost::lexical_cast<T>("0.42278433509846713939348790991759756895784066406008");
|
||||
NTL::RR r1 = boost::lexical_cast<NTL::RR>("0.57721566490153286060651209008240243104215933593992");
|
||||
NTL::RR r2 = boost::lexical_cast<NTL::RR>("0.42278433509846713939348790991759756895784066406008");
|
||||
if(x == 0)
|
||||
{
|
||||
return r1;
|
||||
@@ -70,8 +72,8 @@ T f0(T x, int variant)
|
||||
//
|
||||
// works well over [1.5, 2] but not near 1 :-(
|
||||
//
|
||||
NTL::RR r1 = boost::lexical_cast<T>("0.57721566490153286060651209008240243104215933593992");
|
||||
NTL::RR r2 = boost::lexical_cast<T>("0.42278433509846713939348790991759756895784066406008");
|
||||
NTL::RR r1 = boost::lexical_cast<NTL::RR>("0.57721566490153286060651209008240243104215933593992");
|
||||
NTL::RR r2 = boost::lexical_cast<NTL::RR>("0.42278433509846713939348790991759756895784066406008");
|
||||
if(x == 0)
|
||||
{
|
||||
return r2;
|
||||
@@ -83,48 +85,49 @@ T f0(T x, int variant)
|
||||
return boost::math::lgamma(2-x) / (x * (x - 1));
|
||||
}
|
||||
case 7:
|
||||
//
|
||||
// erf_inv in range [0, 0.5]
|
||||
//
|
||||
if(x == 0)
|
||||
x = boost::math::tools::epsilon<T>() / 64;
|
||||
return boost::math::erf_inv(x) / (x * (x+10));
|
||||
{
|
||||
//
|
||||
// erf_inv in range [0, 0.5]
|
||||
//
|
||||
NTL::RR y = x;
|
||||
if(y == 0)
|
||||
y = boost::math::tools::epsilon<NTL::RR>() / 64;
|
||||
return boost::math::erf_inv(y) / (y * (y+10));
|
||||
}
|
||||
case 8:
|
||||
//
|
||||
// erfc_inv in range [0.25, 0.5]
|
||||
// Use an x-offset of 0.25, and range [0, 0.25]
|
||||
// abs error, auto y-offset.
|
||||
//
|
||||
if(x == 0)
|
||||
x = boost::lexical_cast<T>("1e-5000");
|
||||
return sqrt(-2 * log(x)) / boost::math::erfc_inv(x);
|
||||
{
|
||||
//
|
||||
// erfc_inv in range [0.25, 0.5]
|
||||
// Use an y-offset of 0.25, and range [0, 0.25]
|
||||
// abs error, auto y-offset.
|
||||
//
|
||||
NTL::RR y = x;
|
||||
if(y == 0)
|
||||
y = boost::lexical_cast<NTL::RR>("1e-5000");
|
||||
return sqrt(-2 * log(y)) / boost::math::erfc_inv(y);
|
||||
}
|
||||
case 9:
|
||||
{
|
||||
if(x == 0)
|
||||
x = tiny;
|
||||
double xr = boost::math::tools::real_cast<double>(x);
|
||||
T z = 1/x;
|
||||
double zr = boost::math::tools::real_cast<double>(z);
|
||||
double yr = exp(zr * zr / -2);
|
||||
T y = exp(z * z / -2);
|
||||
return (boost::math::erfc_inv(y)/* - z*/ + log(z) / z) / z;
|
||||
NTL::RR y = 0.5 - x;
|
||||
return boost::math::erfc_inv(y);
|
||||
}
|
||||
case 10:
|
||||
{
|
||||
if(x == 0)
|
||||
x = boost::lexical_cast<T>("1e-5000");
|
||||
T y = exp(-x*x); // sqrt(-log(x)) - 5;
|
||||
return boost::math::erfc_inv(y) / x;
|
||||
NTL::RR x2 = x;
|
||||
if(x2 == 0)
|
||||
x2 = boost::lexical_cast<NTL::RR>("1e-5000");
|
||||
NTL::RR y = exp(-x2*x2); // sqrt(-log(x2)) - 5;
|
||||
return boost::math::erfc_inv(y) / x2;
|
||||
}
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
NTL::RR f(const NTL::RR& x, int variant)
|
||||
{
|
||||
return f0<NTL::RR>(x, variant);
|
||||
}
|
||||
}*/
|
||||
|
||||
void show_extra(
|
||||
const boost::math::tools::polynomial<NTL::RR>& n,
|
||||
|
||||
362
tools/rational_tests.cpp
Normal file
362
tools/rational_tests.cpp
Normal 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;
|
||||
}
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user