diff --git a/doc/concepts.qbk b/doc/concepts.qbk index f3588326b..1cee3d4c3 100644 --- a/doc/concepts.qbk +++ b/doc/concepts.qbk @@ -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 `` and include that header /before/ you include any of the special function headers. diff --git a/doc/internals_overview.qbk b/doc/internals_overview.qbk new file mode 100644 index 000000000..767e4a7de --- /dev/null +++ b/doc/internals_overview.qbk @@ -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). +] diff --git a/doc/math.qbk b/doc/math.qbk index d9558ce58..966659777 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -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] diff --git a/doc/minima.qbk b/doc/minima.qbk new file mode 100644 index 000000000..9bf8b0b3a --- /dev/null +++ b/doc/minima.qbk @@ -0,0 +1,66 @@ +[section:minima Locating Function Minima] + +[caution __caution ] + +[h4 synopsis] + +`` +#include +`` + + template + std::pair brent_find_minima(F f, T min, T max, int digits); + + template + std::pair 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). +] + diff --git a/doc/minimax.qbk b/doc/minimax.qbk new file mode 100644 index 000000000..fa209c97b --- /dev/null +++ b/doc/minimax.qbk @@ -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). +] diff --git a/doc/relative_error.qbk b/doc/relative_error.qbk index 9e0ff91c1..888564d21 100644 --- a/doc/relative_error.qbk +++ b/doc/relative_error.qbk @@ -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] diff --git a/doc/test_data.qbk b/doc/test_data.qbk index de5051f5d..b9dca408a 100644 --- a/doc/test_data.qbk +++ b/doc/test_data.qbk @@ -66,24 +66,24 @@ graphing program (or spreadsheet) for further analysis. template test_data& insert(F func, const parameter_info& arg1, const parameter_info& arg2, const parameter_info& 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 @@ -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 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 + parameter_info 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 + parameter_info 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 + parameter_info 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 + bool get_user_parameter_info(parameter_info& 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 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] diff --git a/include/boost/math/tools/test_data.hpp b/include/boost/math/tools/test_data.hpp index 62938c59e..7cd50a2cd 100644 --- a/include/boost/math/tools/test_data.hpp +++ b/include/boost/math/tools/test_data.hpp @@ -379,7 +379,7 @@ void test_data::create_test_points(std::set& points, const parameter_info< { random_type r = gen(); power_type p = ldexp(static_cast(r), power); - points.insert(truncate_to_float(real_cast(p))); + points.insert(truncate_to_float(real_cast(arg1.z1 + p))); } } break; diff --git a/minimax/f.cpp b/minimax/f.cpp index b35f93dc7..dd0122781 100644 --- a/minimax/f.cpp +++ b/minimax/f.cpp @@ -10,10 +10,9 @@ #include -template -T f0(T x, int variant) +NTL::RR f(const NTL::RR& x, int variant) { - static const T tiny = boost::math::tools::min_value() * 64; + static const NTL::RR tiny = boost::math::tools::min_value() * 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("0.42278433509846713939348790991759756895784066406008") / 3; + return boost::lexical_cast("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("0.57721566490153286060651209008240243104215933593992"); - NTL::RR r2 = boost::lexical_cast("0.42278433509846713939348790991759756895784066406008"); + NTL::RR r1 = boost::lexical_cast("0.57721566490153286060651209008240243104215933593992"); + NTL::RR r2 = boost::lexical_cast("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("0.57721566490153286060651209008240243104215933593992"); - NTL::RR r2 = boost::lexical_cast("0.42278433509846713939348790991759756895784066406008"); + NTL::RR r1 = boost::lexical_cast("0.57721566490153286060651209008240243104215933593992"); + NTL::RR r2 = boost::lexical_cast("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() / 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() / 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("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("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(x); - T z = 1/x; - double zr = boost::math::tools::real_cast(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("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("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(x, variant); -} +}*/ void show_extra( const boost::math::tools::polynomial& n, diff --git a/tools/rational_tests.cpp b/tools/rational_tests.cpp new file mode 100644 index 000000000..57b691743 --- /dev/null +++ b/tools/rational_tests.cpp @@ -0,0 +1,362 @@ + +#include +#include +#include +#include +#include + +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 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 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(0.125), " << i << "),\n" + " static_cast(" << r1 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_polynomial(n" << i << "c, static_cast(0.25), " << i << "),\n" + " static_cast(" << r2 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_polynomial(n" << i << "c, static_cast(0.75), " << i << "),\n" + " static_cast(" << r3 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_polynomial(n" << i << "c, static_cast(1.0f - 1.0f/64.0f), " << i << "),\n" + " static_cast(" << r4 << "L),\n" + " tolerance);\n\n"; + + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_polynomial(n" << i << "c, static_cast(0.125)),\n" + " static_cast(" << r1 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_polynomial(n" << i << "c, static_cast(0.25)),\n" + " static_cast(" << r2 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_polynomial(n" << i << "c, static_cast(0.75)),\n" + " static_cast(" << r3 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_polynomial(n" << i << "c, static_cast(1.0f - 1.0f/64.0f)),\n" + " static_cast(" << r4 << "L),\n" + " tolerance);\n\n"; + + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_polynomial(n" << i << "a, static_cast(0.125)),\n" + " static_cast(" << r1 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_polynomial(n" << i << "a, static_cast(0.25)),\n" + " static_cast(" << r2 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_polynomial(n" << i << "a, static_cast(0.75)),\n" + " static_cast(" << r3 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_polynomial(n" << i << "a, static_cast(1.0f - 1.0f/64.0f)),\n" + " static_cast(" << 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(0.125), " << i << "),\n" + " static_cast(" << r1 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_even_polynomial(n" << i << "c, static_cast(0.25), " << i << "),\n" + " static_cast(" << r2 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_even_polynomial(n" << i << "c, static_cast(0.75), " << i << "),\n" + " static_cast(" << r3 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_even_polynomial(n" << i << "c, static_cast(1.0f - 1.0f/64.0f), " << i << "),\n" + " static_cast(" << r4 << "L),\n" + " tolerance);\n\n"; + + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_even_polynomial(n" << i << "c, static_cast(0.125)),\n" + " static_cast(" << r1 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_even_polynomial(n" << i << "c, static_cast(0.25)),\n" + " static_cast(" << r2 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_even_polynomial(n" << i << "c, static_cast(0.75)),\n" + " static_cast(" << r3 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_even_polynomial(n" << i << "c, static_cast(1.0f - 1.0f/64.0f)),\n" + " static_cast(" << r4 << "L),\n" + " tolerance);\n\n"; + + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_even_polynomial(n" << i << "a, static_cast(0.125)),\n" + " static_cast(" << r1 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_even_polynomial(n" << i << "a, static_cast(0.25)),\n" + " static_cast(" << r2 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_even_polynomial(n" << i << "a, static_cast(0.75)),\n" + " static_cast(" << r3 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_even_polynomial(n" << i << "a, static_cast(1.0f - 1.0f/64.0f)),\n" + " static_cast(" << 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(0.125), " << i << "),\n" + " static_cast(" << r1 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_odd_polynomial(n" << i << "c, static_cast(0.25), " << i << "),\n" + " static_cast(" << r2 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_odd_polynomial(n" << i << "c, static_cast(0.75), " << i << "),\n" + " static_cast(" << r3 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_odd_polynomial(n" << i << "c, static_cast(1.0f - 1.0f/64.0f), " << i << "),\n" + " static_cast(" << r4 << "L),\n" + " tolerance);\n\n"; + + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_odd_polynomial(n" << i << "c, static_cast(0.125)),\n" + " static_cast(" << r1 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_odd_polynomial(n" << i << "c, static_cast(0.25)),\n" + " static_cast(" << r2 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_odd_polynomial(n" << i << "c, static_cast(0.75)),\n" + " static_cast(" << r3 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_odd_polynomial(n" << i << "c, static_cast(1.0f - 1.0f/64.0f)),\n" + " static_cast(" << r4 << "L),\n" + " tolerance);\n\n"; + + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_odd_polynomial(n" << i << "a, static_cast(0.125)),\n" + " static_cast(" << r1 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_odd_polynomial(n" << i << "a, static_cast(0.25)),\n" + " static_cast(" << r2 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_odd_polynomial(n" << i << "a, static_cast(0.75)),\n" + " static_cast(" << r3 << "L),\n" + " tolerance);\n"; + std::cout << + " BOOST_CHECK_CLOSE(\n" + " boost::math::tools::evaluate_odd_polynomial(n" << i << "a, static_cast(1.0f - 1.0f/64.0f)),\n" + " static_cast(" << 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 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(0.125), " << i << "),\n" + " static_cast(" << 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(0.25), " << i << "),\n" + " static_cast(" << 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(0.75), " << i << "),\n" + " static_cast(" << 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(1.0f - 1.0f/64.0f), " << i << "),\n" + " static_cast(" << 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(0.125)),\n" + " static_cast(" << 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(0.25)),\n" + " static_cast(" << 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(0.75)),\n" + " static_cast(" << 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(1.0f - 1.0f/64.0f)),\n" + " static_cast(" << 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(0.125)),\n" + " static_cast(" << 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(0.25)),\n" + " static_cast(" << 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(0.75)),\n" + " static_cast(" << 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(1.0f - 1.0f/64.0f)),\n" + " static_cast(" << r4/r4d << "L),\n" + " tolerance);\n\n"; + } + + return 0; +} + + +