mirror of
https://github.com/boostorg/math.git
synced 2026-01-19 04:22:09 +00:00
Minimax and distribution_explorer to examples. (#620)
* Minimax and distribution_explorer to examples. * Move minimax to tools.
This commit is contained in:
Nick
GitHub
46
tools/minimax/Jamfile.v2
Normal file
46
tools/minimax/Jamfile.v2
Normal file
@@ -0,0 +1,46 @@
|
||||
# Copyright John Maddock 2010
|
||||
# Copyright Paul A. Bristow 2018
|
||||
# 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.
|
||||
# \math_toolkit\libs\math\minimax\jamfile.v2
|
||||
# Runs minimax using multiprecision, (rather than gmp and mpfr)
|
||||
|
||||
# bring in the rules for testing.
|
||||
import modules ;
|
||||
import path ;
|
||||
|
||||
project
|
||||
: requirements
|
||||
<toolset>gcc:<cxxflags>-Wno-missing-braces
|
||||
<toolset>darwin:<cxxflags>-Wno-missing-braces
|
||||
<toolset>acc:<cxxflags>+W2068,2461,2236,4070,4069
|
||||
<toolset>intel-win:<cxxflags>-nologo
|
||||
<toolset>intel-win:<linkflags>-nologo
|
||||
<toolset>msvc:<warnings>all
|
||||
<toolset>msvc:<asynch-exceptions>on
|
||||
<toolset>msvc:<cxxflags>/wd4996
|
||||
<toolset>msvc:<cxxflags>/wd4512
|
||||
<toolset>msvc:<cxxflags>/wd4610
|
||||
<toolset>msvc:<cxxflags>/wd4510
|
||||
<toolset>msvc:<cxxflags>/wd4127
|
||||
<toolset>msvc:<cxxflags>/wd4701 # needed for lexical cast - temporary.
|
||||
<link>static
|
||||
<toolset>borland:<runtime-link>static
|
||||
<include>../../..
|
||||
<define>BOOST_ALL_NO_LIB=1
|
||||
<define>BOOST_UBLAS_UNSUPPORTED_COMPILER=0
|
||||
<include>.
|
||||
<include>../include_private
|
||||
#<include>$(ntl-path)/include
|
||||
;
|
||||
|
||||
#lib mpfr : gmp : <name>mpfr ;
|
||||
|
||||
#lib gmp : : <name>gmp ;
|
||||
|
||||
# exe minimax : f.cpp main.cpp gmp mpfr ;
|
||||
exe minimax : f.cpp main.cpp ;
|
||||
|
||||
install bin : minimax ;
|
||||
|
||||
404
tools/minimax/f.cpp
Normal file
404
tools/minimax/f.cpp
Normal file
@@ -0,0 +1,404 @@
|
||||
// (C) Copyright John Maddock 2006.
|
||||
// Use, modification and distribution are subject to the
|
||||
// Boost Software License, Version 1.0. (See accompanying file
|
||||
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#define L22
|
||||
//#include "../tools/ntl_rr_lanczos.hpp"
|
||||
//#include "../tools/ntl_rr_digamma.hpp"
|
||||
#include "multiprecision.hpp"
|
||||
#include <boost/math/tools/polynomial.hpp>
|
||||
#include <boost/math/special_functions.hpp>
|
||||
#include <boost/math/special_functions/zeta.hpp>
|
||||
#include <boost/math/special_functions/expint.hpp>
|
||||
#include <boost/math/special_functions/lambert_w.hpp>
|
||||
|
||||
#include <cmath>
|
||||
|
||||
|
||||
mp_type f(const mp_type& x, int variant)
|
||||
{
|
||||
static const mp_type tiny = boost::math::tools::min_value<mp_type>() * 64;
|
||||
switch(variant)
|
||||
{
|
||||
case 0:
|
||||
{
|
||||
mp_type x_ = sqrt(x == 0 ? 1e-80 : x);
|
||||
return boost::math::erf(x_) / x_;
|
||||
}
|
||||
case 1:
|
||||
{
|
||||
mp_type x_ = 1 / x;
|
||||
return boost::math::erfc(x_) * x_ / exp(-x_ * x_);
|
||||
}
|
||||
case 2:
|
||||
{
|
||||
return boost::math::erfc(x) * x / exp(-x * x);
|
||||
}
|
||||
case 3:
|
||||
{
|
||||
mp_type 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:
|
||||
//
|
||||
// lgamma(x) = (x-2) * (x + 1) * (c + R(x - 2))
|
||||
//
|
||||
// Works well at 80-bit long double precision, but doesn't
|
||||
// stretch to 128-bit precision.
|
||||
//
|
||||
if(x == 0)
|
||||
{
|
||||
return boost::lexical_cast<mp_type>("0.42278433509846713939348790991759756895784066406008") / 3;
|
||||
}
|
||||
return boost::math::lgamma(x+2) / (x * (x+3));
|
||||
case 5:
|
||||
{
|
||||
//
|
||||
// lgamma in the range [1,2], use:
|
||||
//
|
||||
// lgamma(x) = (x - 1) * (x - 2) * (c + R(x - 1))
|
||||
//
|
||||
// works well over [1, 1.5] but not near 2 :-(
|
||||
//
|
||||
mp_type r1 = boost::lexical_cast<mp_type>("0.57721566490153286060651209008240243104215933593992");
|
||||
mp_type r2 = boost::lexical_cast<mp_type>("0.42278433509846713939348790991759756895784066406008");
|
||||
if(x == 0)
|
||||
{
|
||||
return r1;
|
||||
}
|
||||
if(x == 1)
|
||||
{
|
||||
return r2;
|
||||
}
|
||||
return boost::math::lgamma(x+1) / (x * (x - 1));
|
||||
}
|
||||
case 6:
|
||||
{
|
||||
//
|
||||
// lgamma in the range [1.5,2], use:
|
||||
//
|
||||
// lgamma(x) = (2 - x) * (1 - x) * (c + R(2 - x))
|
||||
//
|
||||
// works well over [1.5, 2] but not near 1 :-(
|
||||
//
|
||||
mp_type r1 = boost::lexical_cast<mp_type>("0.57721566490153286060651209008240243104215933593992");
|
||||
mp_type r2 = boost::lexical_cast<mp_type>("0.42278433509846713939348790991759756895784066406008");
|
||||
if(x == 0)
|
||||
{
|
||||
return r2;
|
||||
}
|
||||
if(x == 1)
|
||||
{
|
||||
return r1;
|
||||
}
|
||||
return boost::math::lgamma(2-x) / (x * (x - 1));
|
||||
}
|
||||
case 7:
|
||||
{
|
||||
//
|
||||
// erf_inv in range [0, 0.5]
|
||||
//
|
||||
mp_type y = x;
|
||||
if(y == 0)
|
||||
y = boost::math::tools::epsilon<mp_type>() / 64;
|
||||
return boost::math::erf_inv(y) / (y * (y+10));
|
||||
}
|
||||
case 8:
|
||||
{
|
||||
//
|
||||
// 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.
|
||||
//
|
||||
mp_type y = x;
|
||||
if(y == 0)
|
||||
y = boost::lexical_cast<mp_type>("1e-5000");
|
||||
return sqrt(-2 * log(y)) / boost::math::erfc_inv(y);
|
||||
}
|
||||
case 9:
|
||||
{
|
||||
mp_type x2 = x;
|
||||
if(x2 == 0)
|
||||
x2 = boost::lexical_cast<mp_type>("1e-5000");
|
||||
mp_type y = exp(-x2*x2); // sqrt(-log(x2)) - 5;
|
||||
return boost::math::erfc_inv(y) / x2;
|
||||
}
|
||||
case 10:
|
||||
{
|
||||
//
|
||||
// Digamma over the interval [1,2], set x-offset to 1
|
||||
// and optimise for absolute error over [0,1].
|
||||
//
|
||||
int current_precision = get_working_precision();
|
||||
if(current_precision < 1000)
|
||||
set_working_precision(1000);
|
||||
//
|
||||
// This value for the root of digamma is calculated using our
|
||||
// differentiated lanczos approximation. It agrees with Cody
|
||||
// to ~ 25 digits and to Morris to 35 digits. See:
|
||||
// TOMS ALGORITHM 708 (Didonato and Morris).
|
||||
// and Math. Comp. 27, 123-127 (1973) by Cody, Strecok and Thacher.
|
||||
//
|
||||
//mp_type root = boost::lexical_cast<mp_type>("1.4616321449683623412626595423257213234331845807102825466429633351908372838889871");
|
||||
//
|
||||
// Actually better to calculate the root on the fly, it appears to be more
|
||||
// accurate: convergence is easier with the 1000-bit value, the approximation
|
||||
// produced agrees with functions.mathworld.com values to 35 digits even quite
|
||||
// near the root.
|
||||
//
|
||||
static boost::math::tools::eps_tolerance<mp_type> tol(1000);
|
||||
static std::uintmax_t max_iter = 1000;
|
||||
mp_type (*pdg)(mp_type) = &boost::math::digamma;
|
||||
static const mp_type root = boost::math::tools::bracket_and_solve_root(pdg, mp_type(1.4), mp_type(1.5), true, tol, max_iter).first;
|
||||
|
||||
mp_type x2 = x;
|
||||
double lim = 1e-65;
|
||||
if(fabs(x2 - root) < lim)
|
||||
{
|
||||
//
|
||||
// This is a problem area:
|
||||
// x2-root suffers cancellation error, so does digamma.
|
||||
// That gets compounded again when Remez calculates the error
|
||||
// function. This cludge seems to stop the worst of the problems:
|
||||
//
|
||||
static const mp_type a = boost::math::digamma(root - lim) / -lim;
|
||||
static const mp_type b = boost::math::digamma(root + lim) / lim;
|
||||
mp_type fract = (x2 - root + lim) / (2*lim);
|
||||
mp_type r = (1-fract) * a + fract * b;
|
||||
std::cout << "In root area: " << r;
|
||||
return r;
|
||||
}
|
||||
mp_type result = boost::math::digamma(x2) / (x2 - root);
|
||||
if(current_precision < 1000)
|
||||
set_working_precision(current_precision);
|
||||
return result;
|
||||
}
|
||||
case 11:
|
||||
// expm1:
|
||||
if(x == 0)
|
||||
{
|
||||
static mp_type lim = 1e-80;
|
||||
static mp_type a = boost::math::expm1(-lim);
|
||||
static mp_type b = boost::math::expm1(lim);
|
||||
static mp_type l = (b-a) / (2 * lim);
|
||||
return l;
|
||||
}
|
||||
return boost::math::expm1(x) / x;
|
||||
case 12:
|
||||
// demo, and test case:
|
||||
return exp(x);
|
||||
case 13:
|
||||
// K(k):
|
||||
{
|
||||
return boost::math::ellint_1(x);
|
||||
}
|
||||
case 14:
|
||||
// K(k)
|
||||
{
|
||||
return boost::math::ellint_1(1-x) / log(x);
|
||||
}
|
||||
case 15:
|
||||
// E(k)
|
||||
{
|
||||
// x = 1-k^2
|
||||
mp_type z = 1 - x * log(x);
|
||||
return boost::math::ellint_2(sqrt(1-x)) / z;
|
||||
}
|
||||
case 16:
|
||||
// Bessel I0(x) over [0,16]:
|
||||
{
|
||||
return boost::math::cyl_bessel_i(0, sqrt(x));
|
||||
}
|
||||
case 17:
|
||||
// Bessel I0(x) over [16,INF]
|
||||
{
|
||||
mp_type z = 1 / (mp_type(1)/16 - x);
|
||||
return boost::math::cyl_bessel_i(0, z) * sqrt(z) / exp(z);
|
||||
}
|
||||
case 18:
|
||||
// Zeta over [0, 1]
|
||||
{
|
||||
return boost::math::zeta(1 - x) * x - x;
|
||||
}
|
||||
case 19:
|
||||
// Zeta over [1, n]
|
||||
{
|
||||
return boost::math::zeta(x) - 1 / (x - 1);
|
||||
}
|
||||
case 20:
|
||||
// Zeta over [a, b] : a >> 1
|
||||
{
|
||||
return log(boost::math::zeta(x) - 1);
|
||||
}
|
||||
case 21:
|
||||
// expint[1] over [0,1]:
|
||||
{
|
||||
mp_type tiny = boost::lexical_cast<mp_type>("1e-5000");
|
||||
mp_type z = (x <= tiny) ? tiny : x;
|
||||
return boost::math::expint(1, z) - z + log(z);
|
||||
}
|
||||
case 22:
|
||||
// expint[1] over [1,N],
|
||||
// Note that x varies from [0,1]:
|
||||
{
|
||||
mp_type z = 1 / x;
|
||||
return boost::math::expint(1, z) * exp(z) * z;
|
||||
}
|
||||
case 23:
|
||||
// expin Ei over [0,R]
|
||||
{
|
||||
static const mp_type root =
|
||||
boost::lexical_cast<mp_type>("0.372507410781366634461991866580119133535689497771654051555657435242200120636201854384926049951548942392");
|
||||
mp_type z = x < (std::numeric_limits<long double>::min)() ? (std::numeric_limits<long double>::min)() : x;
|
||||
return (boost::math::expint(z) - log(z / root)) / (z - root);
|
||||
}
|
||||
case 24:
|
||||
// Expint Ei for large x:
|
||||
{
|
||||
static const mp_type root =
|
||||
boost::lexical_cast<mp_type>("0.372507410781366634461991866580119133535689497771654051555657435242200120636201854384926049951548942392");
|
||||
mp_type z = x < (std::numeric_limits<long double>::min)() ? (std::numeric_limits<long double>::max)() : mp_type(1 / x);
|
||||
return (boost::math::expint(z) - z) * z * exp(-z);
|
||||
//return (boost::math::expint(z) - log(z)) * z * exp(-z);
|
||||
}
|
||||
case 25:
|
||||
// Expint Ei for large x:
|
||||
{
|
||||
return (boost::math::expint(x) - x) * x * exp(-x);
|
||||
}
|
||||
case 26:
|
||||
{
|
||||
//
|
||||
// erf_inv in range [0, 0.5]
|
||||
//
|
||||
mp_type y = x;
|
||||
if(y == 0)
|
||||
y = boost::math::tools::epsilon<mp_type>() / 64;
|
||||
y = sqrt(y);
|
||||
return boost::math::erf_inv(y) / (y);
|
||||
}
|
||||
case 28:
|
||||
{
|
||||
// log1p over [-0.5,0.5]
|
||||
mp_type y = x;
|
||||
if(fabs(y) < 1e-100)
|
||||
y = (y == 0) ? 1e-100 : boost::math::sign(y) * 1e-100;
|
||||
return (boost::math::log1p(y) - y + y * y / 2) / (y);
|
||||
}
|
||||
case 29:
|
||||
{
|
||||
// cbrt over [0.5, 1]
|
||||
return boost::math::cbrt(x);
|
||||
}
|
||||
case 30:
|
||||
{
|
||||
// trigamma over [x,y]
|
||||
mp_type y = x;
|
||||
y = sqrt(y);
|
||||
return boost::math::trigamma(x) * (x * x);
|
||||
}
|
||||
case 31:
|
||||
{
|
||||
// trigamma over [x, INF]
|
||||
if(x == 0) return 1;
|
||||
mp_type y = (x == 0) ? (std::numeric_limits<double>::max)() / 2 : mp_type(1/x);
|
||||
return boost::math::trigamma(y) * y;
|
||||
}
|
||||
case 32:
|
||||
{
|
||||
// I0 over [N, INF]
|
||||
// Don't need to go past x = 1/1000 = 1e-3 for double, or
|
||||
// 1/15000 = 0.0006 for long double, start at 1/7.75=0.13
|
||||
mp_type arg = 1 / x;
|
||||
return sqrt(arg) * exp(-arg) * boost::math::cyl_bessel_i(0, arg);
|
||||
}
|
||||
case 33:
|
||||
{
|
||||
// I0 over [0, N]
|
||||
mp_type xx = sqrt(x) * 2;
|
||||
return (boost::math::cyl_bessel_i(0, xx) - 1) / x;
|
||||
}
|
||||
case 34:
|
||||
{
|
||||
// I1 over [0, N]
|
||||
mp_type xx = sqrt(x) * 2;
|
||||
return (boost::math::cyl_bessel_i(1, xx) * 2 / xx - 1 - x / 2) / (x * x);
|
||||
}
|
||||
case 35:
|
||||
{
|
||||
// I1 over [N, INF]
|
||||
mp_type xx = 1 / x;
|
||||
return boost::math::cyl_bessel_i(1, xx) * sqrt(xx) * exp(-xx);
|
||||
}
|
||||
case 36:
|
||||
{
|
||||
// K0 over [0, 1]
|
||||
mp_type xx = sqrt(x);
|
||||
return boost::math::cyl_bessel_k(0, xx) + log(xx) * boost::math::cyl_bessel_i(0, xx);
|
||||
}
|
||||
case 37:
|
||||
{
|
||||
// K0 over [1, INF]
|
||||
mp_type xx = 1 / x;
|
||||
return boost::math::cyl_bessel_k(0, xx) * exp(xx) * sqrt(xx);
|
||||
}
|
||||
case 38:
|
||||
{
|
||||
// K1 over [0, 1]
|
||||
mp_type xx = sqrt(x);
|
||||
return (boost::math::cyl_bessel_k(1, xx) - log(xx) * boost::math::cyl_bessel_i(1, xx) - 1 / xx) / xx;
|
||||
}
|
||||
case 39:
|
||||
{
|
||||
// K1 over [1, INF]
|
||||
mp_type xx = 1 / x;
|
||||
return boost::math::cyl_bessel_k(1, xx) * sqrt(xx) * exp(xx);
|
||||
}
|
||||
// Lambert W0
|
||||
case 40:
|
||||
return boost::math::lambert_w0(x);
|
||||
case 41:
|
||||
{
|
||||
if (x == 0)
|
||||
return 1;
|
||||
return boost::math::lambert_w0(x) / x;
|
||||
}
|
||||
case 42:
|
||||
{
|
||||
static const mp_type e1 = exp(mp_type(-1));
|
||||
return x / -boost::math::lambert_w0(-e1 + x);
|
||||
}
|
||||
case 43:
|
||||
{
|
||||
mp_type xx = 1 / x;
|
||||
return 1 / boost::math::lambert_w0(xx);
|
||||
}
|
||||
case 44:
|
||||
{
|
||||
mp_type ex = exp(x);
|
||||
return boost::math::lambert_w0(ex) - x;
|
||||
}
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
void show_extra(
|
||||
const boost::math::tools::polynomial<mp_type>& n,
|
||||
const boost::math::tools::polynomial<mp_type>& d,
|
||||
const mp_type& x_offset,
|
||||
const mp_type& y_offset,
|
||||
int variant)
|
||||
{
|
||||
switch(variant)
|
||||
{
|
||||
default:
|
||||
// do nothing here...
|
||||
;
|
||||
}
|
||||
}
|
||||
|
||||
650
tools/minimax/main.cpp
Normal file
650
tools/minimax/main.cpp
Normal file
@@ -0,0 +1,650 @@
|
||||
// (C) Copyright John Maddock 2006.
|
||||
// Use, modification and distribution are subject to the
|
||||
// Boost Software License, Version 1.0. (See accompanying file
|
||||
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#define BOOST_TEST_MODULE foobar
|
||||
#define BOOST_UBLAS_TYPE_CHECK_EPSILON (type_traits<real_type>::type_sqrt (boost::math::tools::epsilon <real_type>()))
|
||||
#define BOOST_UBLAS_TYPE_CHECK_MIN (type_traits<real_type>::type_sqrt ( boost::math::tools::min_value<real_type>()))
|
||||
#define BOOST_UBLAS_NDEBUG
|
||||
|
||||
#include "multiprecision.hpp"
|
||||
|
||||
#include <boost/math/tools/remez.hpp>
|
||||
#include <boost/math/tools/test.hpp>
|
||||
#include <boost/math/special_functions/binomial.hpp>
|
||||
#include <boost/spirit/include/classic_core.hpp>
|
||||
#include <boost/spirit/include/classic_actor.hpp>
|
||||
#include <boost/lexical_cast.hpp>
|
||||
#include <iostream>
|
||||
#include <iomanip>
|
||||
#include <string>
|
||||
#include <boost/test/included/unit_test.hpp> // for test_main
|
||||
#include <boost/multiprecision/cpp_bin_float.hpp>
|
||||
|
||||
|
||||
extern mp_type f(const mp_type& x, int variant);
|
||||
extern void show_extra(
|
||||
const boost::math::tools::polynomial<mp_type>& n,
|
||||
const boost::math::tools::polynomial<mp_type>& d,
|
||||
const mp_type& x_offset,
|
||||
const mp_type& y_offset,
|
||||
int variant);
|
||||
|
||||
using namespace boost::spirit::classic;
|
||||
|
||||
mp_type a(0), b(1); // range to optimise over
|
||||
bool rel_error(true);
|
||||
bool pin(false);
|
||||
int orderN(3);
|
||||
int orderD(1);
|
||||
int target_precision = boost::math::tools::digits<long double>();
|
||||
int working_precision = target_precision * 2;
|
||||
bool started(false);
|
||||
int variant(0);
|
||||
int skew(0);
|
||||
int brake(50);
|
||||
mp_type x_offset(0), y_offset(0), x_scale(1);
|
||||
bool auto_offset_y;
|
||||
|
||||
boost::shared_ptr<boost::math::tools::remez_minimax<mp_type> > p_remez;
|
||||
|
||||
mp_type the_function(const mp_type& val)
|
||||
{
|
||||
return f(x_scale * (val + x_offset), variant) + y_offset;
|
||||
}
|
||||
|
||||
void step_some(unsigned count)
|
||||
{
|
||||
try{
|
||||
set_working_precision(working_precision);
|
||||
if(!started)
|
||||
{
|
||||
//
|
||||
// If we have an automatic y-offset calculate it now:
|
||||
//
|
||||
if(auto_offset_y)
|
||||
{
|
||||
mp_type fa, fb, fm;
|
||||
fa = f(x_scale * (a + x_offset), variant);
|
||||
fb = f(x_scale * (b + x_offset), variant);
|
||||
fm = f(x_scale * ((a+b)/2 + x_offset), variant);
|
||||
y_offset = -(fa + fb + fm) / 3;
|
||||
set_output_precision(5);
|
||||
std::cout << "Setting auto-y-offset to " << y_offset << std::endl;
|
||||
}
|
||||
//
|
||||
// Truncate offsets to float precision:
|
||||
//
|
||||
x_offset = round_to_precision(x_offset, 20);
|
||||
y_offset = round_to_precision(y_offset, 20);
|
||||
//
|
||||
// Construct new Remez state machine:
|
||||
//
|
||||
p_remez.reset(new boost::math::tools::remez_minimax<mp_type>(
|
||||
&the_function,
|
||||
orderN, orderD,
|
||||
a, b,
|
||||
pin,
|
||||
rel_error,
|
||||
skew,
|
||||
working_precision));
|
||||
std::cout << "Max error in interpolated form: " << std::setprecision(3) << std::scientific << boost::math::tools::real_cast<double>(p_remez->max_error()) << std::endl;
|
||||
//
|
||||
// Signal that we've started:
|
||||
//
|
||||
started = true;
|
||||
}
|
||||
unsigned i;
|
||||
for(i = 0; i < count; ++i)
|
||||
{
|
||||
std::cout << "Stepping..." << std::endl;
|
||||
p_remez->set_brake(brake);
|
||||
mp_type r = p_remez->iterate();
|
||||
set_output_precision(3);
|
||||
std::cout
|
||||
<< "Maximum Deviation Found: " << std::setprecision(3) << std::scientific << boost::math::tools::real_cast<double>(p_remez->max_error()) << std::endl
|
||||
<< "Expected Error Term: " << std::setprecision(3) << std::scientific << boost::math::tools::real_cast<double>(p_remez->error_term()) << std::endl
|
||||
<< "Maximum Relative Change in Control Points: " << std::setprecision(3) << std::scientific << boost::math::tools::real_cast<double>(r) << std::endl;
|
||||
}
|
||||
}
|
||||
catch(const std::exception& e)
|
||||
{
|
||||
std::cout << "Step failed with exception: " << e.what() << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
void step(const char*, const char*)
|
||||
{
|
||||
step_some(1);
|
||||
}
|
||||
|
||||
void show(const char*, const char*)
|
||||
{
|
||||
set_working_precision(working_precision);
|
||||
if(started)
|
||||
{
|
||||
boost::math::tools::polynomial<mp_type> n = p_remez->numerator();
|
||||
boost::math::tools::polynomial<mp_type> d = p_remez->denominator();
|
||||
std::vector<mp_type> cn = n.chebyshev();
|
||||
std::vector<mp_type> cd = d.chebyshev();
|
||||
int prec = 2 + (target_precision * 3010LL)/10000;
|
||||
std::cout << std::scientific << std::setprecision(prec);
|
||||
set_output_precision(prec);
|
||||
boost::numeric::ublas::vector<mp_type> v = p_remez->zero_points();
|
||||
|
||||
std::cout << " Zeros = {\n";
|
||||
unsigned i;
|
||||
for(i = 0; i < v.size(); ++i)
|
||||
{
|
||||
std::cout << " " << v[i] << std::endl;
|
||||
}
|
||||
std::cout << " }\n";
|
||||
|
||||
v = p_remez->chebyshev_points();
|
||||
std::cout << " Chebeshev Control Points = {\n";
|
||||
for(i = 0; i < v.size(); ++i)
|
||||
{
|
||||
std::cout << " " << v[i] << std::endl;
|
||||
}
|
||||
std::cout << " }\n";
|
||||
|
||||
std::cout << "X offset: " << x_offset << std::endl;
|
||||
std::cout << "X scale: " << x_scale << std::endl;
|
||||
std::cout << "Y offset: " << y_offset << std::endl;
|
||||
|
||||
std::cout << "P = {";
|
||||
for(i = 0; i < n.size(); ++i)
|
||||
{
|
||||
std::cout << " " << n[i] << "L," << std::endl;
|
||||
}
|
||||
std::cout << " }\n";
|
||||
|
||||
std::cout << "Q = {";
|
||||
for(i = 0; i < d.size(); ++i)
|
||||
{
|
||||
std::cout << " " << d[i] << "L," << std::endl;
|
||||
}
|
||||
std::cout << " }\n";
|
||||
|
||||
std::cout << "CP = {";
|
||||
for(i = 0; i < cn.size(); ++i)
|
||||
{
|
||||
std::cout << " " << cn[i] << "L," << std::endl;
|
||||
}
|
||||
std::cout << " }\n";
|
||||
|
||||
std::cout << "CQ = {";
|
||||
for(i = 0; i < cd.size(); ++i)
|
||||
{
|
||||
std::cout << " " << cd[i] << "L," << std::endl;
|
||||
}
|
||||
std::cout << " }\n";
|
||||
|
||||
show_extra(n, d, x_offset, y_offset, variant);
|
||||
}
|
||||
else
|
||||
{
|
||||
std::cerr << "Nothing to display" << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
void do_graph(unsigned points)
|
||||
{
|
||||
set_working_precision(working_precision);
|
||||
mp_type step = (b - a) / (points - 1);
|
||||
mp_type x = a;
|
||||
while(points > 1)
|
||||
{
|
||||
set_output_precision(10);
|
||||
std::cout << std::setprecision(10) << std::setw(30) << std::left
|
||||
<< boost::lexical_cast<std::string>(x) << the_function(x) << std::endl;
|
||||
--points;
|
||||
x += step;
|
||||
}
|
||||
std::cout << std::setprecision(10) << std::setw(30) << std::left
|
||||
<< boost::lexical_cast<std::string>(b) << the_function(b) << std::endl;
|
||||
}
|
||||
|
||||
void graph(const char*, const char*)
|
||||
{
|
||||
do_graph(3);
|
||||
}
|
||||
|
||||
template <class T>
|
||||
mp_type convert_to_rr(const T& val)
|
||||
{
|
||||
return val;
|
||||
}
|
||||
template <class Backend, boost::multiprecision::expression_template_option ET>
|
||||
mp_type convert_to_rr(const boost::multiprecision::number<Backend, ET>& val)
|
||||
{
|
||||
return boost::lexical_cast<mp_type>(val.str());
|
||||
}
|
||||
|
||||
template <class T>
|
||||
void do_test(T, const char* name)
|
||||
{
|
||||
set_working_precision(working_precision);
|
||||
if(started)
|
||||
{
|
||||
//
|
||||
// We want to test the approximation at fixed precision:
|
||||
// either float, double or long double. Begin by getting the
|
||||
// polynomials:
|
||||
//
|
||||
boost::math::tools::polynomial<T> n, d;
|
||||
boost::math::tools::polynomial<mp_type> nr, dr;
|
||||
nr = p_remez->numerator();
|
||||
dr = p_remez->denominator();
|
||||
n = nr;
|
||||
d = dr;
|
||||
|
||||
std::vector<mp_type> cn1, cd1;
|
||||
cn1 = nr.chebyshev();
|
||||
cd1 = dr.chebyshev();
|
||||
std::vector<T> cn, cd;
|
||||
for(unsigned i = 0; i < cn1.size(); ++i)
|
||||
{
|
||||
cn.push_back(boost::math::tools::real_cast<T>(cn1[i]));
|
||||
}
|
||||
for(unsigned i = 0; i < cd1.size(); ++i)
|
||||
{
|
||||
cd.push_back(boost::math::tools::real_cast<T>(cd1[i]));
|
||||
}
|
||||
//
|
||||
// We'll test at the Chebeshev control points which is where
|
||||
// (in theory) the largest deviation should occur. For good
|
||||
// measure we'll test at the zeros as well:
|
||||
//
|
||||
boost::numeric::ublas::vector<mp_type>
|
||||
zeros(p_remez->zero_points()),
|
||||
cheb(p_remez->chebyshev_points());
|
||||
|
||||
mp_type max_error(0), cheb_max_error(0);
|
||||
|
||||
//
|
||||
// Do the tests at the zeros:
|
||||
//
|
||||
std::cout << "Starting tests at " << name << " precision...\n";
|
||||
std::cout << "Absissa Error (Poly) Error (Cheb)\n";
|
||||
for(unsigned i = 0; i < zeros.size(); ++i)
|
||||
{
|
||||
mp_type true_result = the_function(zeros[i]);
|
||||
T absissa = boost::math::tools::real_cast<T>(zeros[i]);
|
||||
mp_type test_result = convert_to_rr(n.evaluate(absissa) / d.evaluate(absissa));
|
||||
mp_type cheb_result = convert_to_rr(boost::math::tools::evaluate_chebyshev(cn, absissa) / boost::math::tools::evaluate_chebyshev(cd, absissa));
|
||||
mp_type err, cheb_err;
|
||||
if(rel_error)
|
||||
{
|
||||
err = boost::math::tools::relative_error(test_result, true_result);
|
||||
cheb_err = boost::math::tools::relative_error(cheb_result, true_result);
|
||||
}
|
||||
else
|
||||
{
|
||||
err = fabs(test_result - true_result);
|
||||
cheb_err = fabs(cheb_result - true_result);
|
||||
}
|
||||
if(err > max_error)
|
||||
max_error = err;
|
||||
if(cheb_err > cheb_max_error)
|
||||
cheb_max_error = cheb_err;
|
||||
std::cout << std::setprecision(6) << std::setw(15) << std::left << absissa
|
||||
<< std::setw(15) << std::left << boost::math::tools::real_cast<T>(err) << boost::math::tools::real_cast<T>(cheb_err) << std::endl;
|
||||
}
|
||||
//
|
||||
// Do the tests at the Chebeshev control points:
|
||||
//
|
||||
for(unsigned i = 0; i < cheb.size(); ++i)
|
||||
{
|
||||
mp_type true_result = the_function(cheb[i]);
|
||||
T absissa = boost::math::tools::real_cast<T>(cheb[i]);
|
||||
mp_type test_result = convert_to_rr(n.evaluate(absissa) / d.evaluate(absissa));
|
||||
mp_type cheb_result = convert_to_rr(boost::math::tools::evaluate_chebyshev(cn, absissa) / boost::math::tools::evaluate_chebyshev(cd, absissa));
|
||||
mp_type err, cheb_err;
|
||||
if(rel_error)
|
||||
{
|
||||
err = boost::math::tools::relative_error(test_result, true_result);
|
||||
cheb_err = boost::math::tools::relative_error(cheb_result, true_result);
|
||||
}
|
||||
else
|
||||
{
|
||||
err = fabs(test_result - true_result);
|
||||
cheb_err = fabs(cheb_result - true_result);
|
||||
}
|
||||
if(err > max_error)
|
||||
max_error = err;
|
||||
std::cout << std::setprecision(6) << std::setw(15) << std::left << absissa
|
||||
<< std::setw(15) << std::left << boost::math::tools::real_cast<T>(err) <<
|
||||
boost::math::tools::real_cast<T>(cheb_err) << std::endl;
|
||||
}
|
||||
std::string msg = "Max Error found at ";
|
||||
msg += name;
|
||||
msg += " precision = ";
|
||||
msg.append(62 - 17 - msg.size(), ' ');
|
||||
std::cout << msg << std::setprecision(6) << "Poly: " << std::setw(20) << std::left
|
||||
<< boost::math::tools::real_cast<T>(max_error) << "Cheb: " << boost::math::tools::real_cast<T>(cheb_max_error) << std::endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
std::cout << "Nothing to test: try converging an approximation first!!!" << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
void test_float(const char*, const char*)
|
||||
{
|
||||
do_test(float(0), "float");
|
||||
}
|
||||
|
||||
void test_double(const char*, const char*)
|
||||
{
|
||||
do_test(double(0), "double");
|
||||
}
|
||||
|
||||
void test_long(const char*, const char*)
|
||||
{
|
||||
do_test((long double)(0), "long double");
|
||||
}
|
||||
|
||||
void test_float80(const char*, const char*)
|
||||
{
|
||||
do_test((boost::multiprecision::cpp_bin_float_double_extended)(0), "float80");
|
||||
}
|
||||
|
||||
void test_float128(const char*, const char*)
|
||||
{
|
||||
do_test((boost::multiprecision::cpp_bin_float_quad)(0), "float128");
|
||||
}
|
||||
|
||||
void test_all(const char*, const char*)
|
||||
{
|
||||
do_test(float(0), "float");
|
||||
do_test(double(0), "double");
|
||||
do_test((long double)(0), "long double");
|
||||
}
|
||||
|
||||
template <class T>
|
||||
void do_test_n(T, const char* name, unsigned count)
|
||||
{
|
||||
set_working_precision(working_precision);
|
||||
if(started)
|
||||
{
|
||||
//
|
||||
// We want to test the approximation at fixed precision:
|
||||
// either float, double or long double. Begin by getting the
|
||||
// polynomials:
|
||||
//
|
||||
boost::math::tools::polynomial<T> n, d;
|
||||
boost::math::tools::polynomial<mp_type> nr, dr;
|
||||
nr = p_remez->numerator();
|
||||
dr = p_remez->denominator();
|
||||
n = nr;
|
||||
d = dr;
|
||||
|
||||
std::vector<mp_type> cn1, cd1;
|
||||
cn1 = nr.chebyshev();
|
||||
cd1 = dr.chebyshev();
|
||||
std::vector<T> cn, cd;
|
||||
for(unsigned i = 0; i < cn1.size(); ++i)
|
||||
{
|
||||
cn.push_back(boost::math::tools::real_cast<T>(cn1[i]));
|
||||
}
|
||||
for(unsigned i = 0; i < cd1.size(); ++i)
|
||||
{
|
||||
cd.push_back(boost::math::tools::real_cast<T>(cd1[i]));
|
||||
}
|
||||
|
||||
mp_type max_error(0), max_cheb_error(0);
|
||||
mp_type step = (b - a) / count;
|
||||
|
||||
//
|
||||
// Do the tests at the zeros:
|
||||
//
|
||||
std::cout << "Starting tests at " << name << " precision...\n";
|
||||
std::cout << "Absissa Error (poly) Error (Cheb)\n";
|
||||
for(mp_type x = a; x <= b; x += step)
|
||||
{
|
||||
mp_type true_result = the_function(x);
|
||||
//std::cout << true_result << std::endl;
|
||||
T absissa = boost::math::tools::real_cast<T>(x);
|
||||
mp_type test_result = convert_to_rr(n.evaluate(absissa) / d.evaluate(absissa));
|
||||
//std::cout << test_result << std::endl;
|
||||
mp_type cheb_result = convert_to_rr(boost::math::tools::evaluate_chebyshev(cn, absissa) / boost::math::tools::evaluate_chebyshev(cd, absissa));
|
||||
//std::cout << cheb_result << std::endl;
|
||||
mp_type err, cheb_err;
|
||||
if(rel_error)
|
||||
{
|
||||
err = boost::math::tools::relative_error(test_result, true_result);
|
||||
cheb_err = boost::math::tools::relative_error(cheb_result, true_result);
|
||||
}
|
||||
else
|
||||
{
|
||||
err = fabs(test_result - true_result);
|
||||
cheb_err = fabs(cheb_result - true_result);
|
||||
}
|
||||
if(err > max_error)
|
||||
max_error = err;
|
||||
if(cheb_err > max_cheb_error)
|
||||
max_cheb_error = cheb_err;
|
||||
std::cout << std::setprecision(6) << std::setw(15) << std::left << boost::math::tools::real_cast<double>(absissa)
|
||||
<< (test_result < true_result ? "-" : "") << std::setw(20) << std::left
|
||||
<< boost::math::tools::real_cast<double>(err)
|
||||
<< boost::math::tools::real_cast<double>(cheb_err) << std::endl;
|
||||
}
|
||||
std::string msg = "Max Error found at ";
|
||||
msg += name;
|
||||
msg += " precision = ";
|
||||
//msg.append(62 - 17 - msg.size(), ' ');
|
||||
std::cout << msg << "Poly: " << std::setprecision(6)
|
||||
//<< std::setw(15) << std::left
|
||||
<< boost::math::tools::real_cast<T>(max_error)
|
||||
<< " Cheb: " << boost::math::tools::real_cast<T>(max_cheb_error) << std::endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
std::cout << "Nothing to test: try converging an approximation first!!!" << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
void test_n(unsigned n)
|
||||
{
|
||||
do_test_n(mp_type(), "mp_type", n);
|
||||
}
|
||||
|
||||
void test_float_n(unsigned n)
|
||||
{
|
||||
do_test_n(float(0), "float", n);
|
||||
}
|
||||
|
||||
void test_double_n(unsigned n)
|
||||
{
|
||||
do_test_n(double(0), "double", n);
|
||||
}
|
||||
|
||||
void test_long_n(unsigned n)
|
||||
{
|
||||
do_test_n((long double)(0), "long double", n);
|
||||
}
|
||||
|
||||
void test_float80_n(unsigned n)
|
||||
{
|
||||
do_test_n((boost::multiprecision::cpp_bin_float_double_extended)(0), "float80", n);
|
||||
}
|
||||
|
||||
void test_float128_n(unsigned n)
|
||||
{
|
||||
do_test_n((boost::multiprecision::cpp_bin_float_quad)(0), "float128", n);
|
||||
}
|
||||
|
||||
void rotate(const char*, const char*)
|
||||
{
|
||||
if(p_remez)
|
||||
{
|
||||
p_remez->rotate();
|
||||
}
|
||||
else
|
||||
{
|
||||
std::cerr << "Nothing to rotate" << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
void rescale(const char*, const char*)
|
||||
{
|
||||
if(p_remez)
|
||||
{
|
||||
p_remez->rescale(a, b);
|
||||
}
|
||||
else
|
||||
{
|
||||
std::cerr << "Nothing to rescale" << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
void graph_poly(const char*, const char*)
|
||||
{
|
||||
int i = 50;
|
||||
set_working_precision(working_precision);
|
||||
if(started)
|
||||
{
|
||||
//
|
||||
// We want to test the approximation at fixed precision:
|
||||
// either float, double or long double. Begin by getting the
|
||||
// polynomials:
|
||||
//
|
||||
boost::math::tools::polynomial<mp_type> n, d;
|
||||
n = p_remez->numerator();
|
||||
d = p_remez->denominator();
|
||||
|
||||
mp_type max_error(0);
|
||||
mp_type step = (b - a) / i;
|
||||
|
||||
std::cout << "Evaluating Numerator...\n";
|
||||
mp_type val;
|
||||
for(val = a; val <= b; val += step)
|
||||
std::cout << n.evaluate(val) << std::endl;
|
||||
std::cout << "Evaluating Denominator...\n";
|
||||
for(val = a; val <= b; val += step)
|
||||
std::cout << d.evaluate(val) << std::endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
std::cout << "Nothing to test: try converging an approximation first!!!" << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
BOOST_AUTO_TEST_CASE( test_main )
|
||||
{
|
||||
std::string line;
|
||||
real_parser<long double/*mp_type*/ > const rr_p;
|
||||
while(std::getline(std::cin, line))
|
||||
{
|
||||
if(parse(line.c_str(), str_p("quit"), space_p).full)
|
||||
return;
|
||||
if(false == parse(line.c_str(),
|
||||
(
|
||||
|
||||
str_p("range")[assign_a(started, false)] && real_p[assign_a(a)] && real_p[assign_a(b)]
|
||||
||
|
||||
str_p("relative")[assign_a(started, false)][assign_a(rel_error, true)]
|
||||
||
|
||||
str_p("absolute")[assign_a(started, false)][assign_a(rel_error, false)]
|
||||
||
|
||||
str_p("pin")[assign_a(started, false)] && str_p("true")[assign_a(pin, true)]
|
||||
||
|
||||
str_p("pin")[assign_a(started, false)] && str_p("false")[assign_a(pin, false)]
|
||||
||
|
||||
str_p("pin")[assign_a(started, false)] && str_p("1")[assign_a(pin, true)]
|
||||
||
|
||||
str_p("pin")[assign_a(started, false)] && str_p("0")[assign_a(pin, false)]
|
||||
||
|
||||
str_p("pin")[assign_a(started, false)][assign_a(pin, true)]
|
||||
||
|
||||
str_p("order")[assign_a(started, false)] && uint_p[assign_a(orderN)] && uint_p[assign_a(orderD)]
|
||||
||
|
||||
str_p("order")[assign_a(started, false)] && uint_p[assign_a(orderN)]
|
||||
||
|
||||
str_p("target-precision") && uint_p[assign_a(target_precision)]
|
||||
||
|
||||
str_p("working-precision")[assign_a(started, false)] && uint_p[assign_a(working_precision)]
|
||||
||
|
||||
str_p("variant")[assign_a(started, false)] && int_p[assign_a(variant)]
|
||||
||
|
||||
str_p("skew")[assign_a(started, false)] && int_p[assign_a(skew)]
|
||||
||
|
||||
str_p("brake") && int_p[assign_a(brake)]
|
||||
||
|
||||
str_p("step") && int_p[&step_some]
|
||||
||
|
||||
str_p("step")[&step]
|
||||
||
|
||||
str_p("poly")[&graph_poly]
|
||||
||
|
||||
str_p("info")[&show]
|
||||
||
|
||||
str_p("graph") && uint_p[&do_graph]
|
||||
||
|
||||
str_p("graph")[&graph]
|
||||
||
|
||||
str_p("x-offset") && real_p[assign_a(x_offset)]
|
||||
||
|
||||
str_p("x-scale") && real_p[assign_a(x_scale)]
|
||||
||
|
||||
str_p("y-offset") && str_p("auto")[assign_a(auto_offset_y, true)]
|
||||
||
|
||||
str_p("y-offset") && real_p[assign_a(y_offset)][assign_a(auto_offset_y, false)]
|
||||
||
|
||||
str_p("test") && str_p("float80") && uint_p[&test_float80_n]
|
||||
||
|
||||
str_p("test") && str_p("float80")[&test_float80]
|
||||
||
|
||||
str_p("test") && str_p("float128") && uint_p[&test_float128_n]
|
||||
||
|
||||
str_p("test") && str_p("float128")[&test_float128]
|
||||
||
|
||||
str_p("test") && str_p("float") && uint_p[&test_float_n]
|
||||
||
|
||||
str_p("test") && str_p("float")[&test_float]
|
||||
||
|
||||
str_p("test") && str_p("double") && uint_p[&test_double_n]
|
||||
||
|
||||
str_p("test") && str_p("double")[&test_double]
|
||||
||
|
||||
str_p("test") && str_p("long") && uint_p[&test_long_n]
|
||||
||
|
||||
str_p("test") && str_p("long")[&test_long]
|
||||
||
|
||||
str_p("test") && str_p("all")[&test_all]
|
||||
||
|
||||
str_p("test") && uint_p[&test_n]
|
||||
||
|
||||
str_p("rotate")[&rotate]
|
||||
||
|
||||
str_p("rescale") && real_p[assign_a(a)] && real_p[assign_a(b)] && epsilon_p[&rescale]
|
||||
|
||||
), space_p).full)
|
||||
{
|
||||
std::cout << "Unable to parse directive: \"" << line << "\"" << std::endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
std::cout << "Variant = " << variant << std::endl;
|
||||
std::cout << "range = [" << a << "," << b << "]" << std::endl;
|
||||
std::cout << "Relative Error = " << rel_error << std::endl;
|
||||
std::cout << "Pin to Origin = " << pin << std::endl;
|
||||
std::cout << "Order (Num/Denom) = " << orderN << "/" << orderD << std::endl;
|
||||
std::cout << "Target Precision = " << target_precision << std::endl;
|
||||
std::cout << "Working Precision = " << working_precision << std::endl;
|
||||
std::cout << "Skew = " << skew << std::endl;
|
||||
std::cout << "Brake = " << brake << std::endl;
|
||||
std::cout << "X Offset = " << x_offset << std::endl;
|
||||
std::cout << "X scale = " << x_scale << std::endl;
|
||||
std::cout << "Y Offset = ";
|
||||
if(auto_offset_y)
|
||||
std::cout << "Auto (";
|
||||
std::cout << y_offset;
|
||||
if(auto_offset_y)
|
||||
std::cout << ")";
|
||||
std::cout << std::endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
224
tools/minimax/multiprecision.hpp
Normal file
224
tools/minimax/multiprecision.hpp
Normal file
@@ -0,0 +1,224 @@
|
||||
// (C) Copyright John Maddock 2015.
|
||||
// Use, modification and distribution are subject to the
|
||||
// Boost Software License, Version 1.0. (See accompanying file
|
||||
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
|
||||
|
||||
#ifndef BOOST_REMEZ_MULTIPRECISION_HPP
|
||||
#define BOOST_REMEZ_MULTIPRECISION_HPP
|
||||
|
||||
#include <boost/multiprecision/cpp_bin_float.hpp>
|
||||
|
||||
#ifdef USE_NTL
|
||||
#include <boost/math/bindings/rr.hpp>
|
||||
namespace std {
|
||||
using boost::math::ntl::pow;
|
||||
} // workaround for spirit parser.
|
||||
|
||||
typedef boost::math::ntl::RR mp_type;
|
||||
|
||||
inline void set_working_precision(int n)
|
||||
{
|
||||
NTL::RR::SetPrecision(working_precision);
|
||||
}
|
||||
|
||||
inline int get_working_precision()
|
||||
{
|
||||
return mp_type::precision(working_precision);
|
||||
}
|
||||
|
||||
inline void set_output_precision(int n)
|
||||
{
|
||||
NTL::RR::SetOutputPrecision(n);
|
||||
}
|
||||
|
||||
inline mp_type round_to_precision(mp_type m, int bits)
|
||||
{
|
||||
return NTL::RoundToPrecision(m.value(), bits);
|
||||
}
|
||||
|
||||
|
||||
namespace boost {
|
||||
namespace math {
|
||||
namespace tools {
|
||||
|
||||
template <>
|
||||
inline boost::multiprecision::cpp_bin_float_double_extended real_cast<boost::multiprecision::cpp_bin_float_double_extended, mp_type>(mp_type val)
|
||||
{
|
||||
unsigned p = NTL::RR::OutputPrecision();
|
||||
NTL::RR::SetOutputPrecision(20);
|
||||
boost::multiprecision::cpp_bin_float_double_extended r = boost::lexical_cast<boost::multiprecision::cpp_bin_float_double_extended>(val);
|
||||
NTL::RR::SetOutputPrecision(p);
|
||||
return r;
|
||||
}
|
||||
template <>
|
||||
inline boost::multiprecision::cpp_bin_float_quad real_cast<boost::multiprecision::cpp_bin_float_quad, mp_type>(mp_type val)
|
||||
{
|
||||
unsigned p = NTL::RR::OutputPrecision();
|
||||
NTL::RR::SetOutputPrecision(35);
|
||||
boost::multiprecision::cpp_bin_float_quad r = boost::lexical_cast<boost::multiprecision::cpp_bin_float_quad>(val);
|
||||
NTL::RR::SetOutputPrecision(p);
|
||||
return r;
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#elif defined(USE_CPP_BIN_FLOAT_100)
|
||||
|
||||
#include <boost/multiprecision/cpp_bin_float.hpp>
|
||||
|
||||
typedef boost::multiprecision::cpp_bin_float_100 mp_type;
|
||||
|
||||
inline void set_working_precision(int n)
|
||||
{
|
||||
}
|
||||
|
||||
inline void set_output_precision(int n)
|
||||
{
|
||||
std::cout << std::setprecision(n);
|
||||
std::cerr << std::setprecision(n);
|
||||
}
|
||||
|
||||
inline mp_type round_to_precision(mp_type m, int bits)
|
||||
{
|
||||
int i;
|
||||
mp_type f = frexp(m, &i);
|
||||
f = ldexp(f, bits);
|
||||
i -= bits;
|
||||
f = floor(f);
|
||||
return ldexp(f, i);
|
||||
}
|
||||
|
||||
inline int get_working_precision()
|
||||
{
|
||||
return std::numeric_limits<mp_type>::digits;
|
||||
}
|
||||
|
||||
namespace boost {
|
||||
namespace math {
|
||||
namespace tools {
|
||||
|
||||
template <>
|
||||
inline boost::multiprecision::cpp_bin_float_double_extended real_cast<boost::multiprecision::cpp_bin_float_double_extended, mp_type>(mp_type val)
|
||||
{
|
||||
return boost::multiprecision::cpp_bin_float_double_extended(val);
|
||||
}
|
||||
template <>
|
||||
inline boost::multiprecision::cpp_bin_float_quad real_cast<boost::multiprecision::cpp_bin_float_quad, mp_type>(mp_type val)
|
||||
{
|
||||
return boost::multiprecision::cpp_bin_float_quad(val);
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#elif defined(USE_MPFR_100)
|
||||
|
||||
#include <boost/multiprecision/mpfr.hpp>
|
||||
|
||||
typedef boost::multiprecision::mpfr_float_100 mp_type;
|
||||
|
||||
inline void set_working_precision(int n)
|
||||
{
|
||||
}
|
||||
|
||||
inline void set_output_precision(int n)
|
||||
{
|
||||
std::cout << std::setprecision(n);
|
||||
std::cerr << std::setprecision(n);
|
||||
}
|
||||
|
||||
inline mp_type round_to_precision(mp_type m, int bits)
|
||||
{
|
||||
mpfr_prec_round(m.backend().data(), bits, MPFR_RNDN);
|
||||
return m;
|
||||
}
|
||||
|
||||
inline int get_working_precision()
|
||||
{
|
||||
return std::numeric_limits<mp_type>::digits;
|
||||
}
|
||||
|
||||
namespace boost {
|
||||
namespace math {
|
||||
namespace tools {
|
||||
|
||||
template <>
|
||||
inline boost::multiprecision::cpp_bin_float_double_extended real_cast<boost::multiprecision::cpp_bin_float_double_extended, mp_type>(mp_type val)
|
||||
{
|
||||
return boost::multiprecision::cpp_bin_float_double_extended(val);
|
||||
}
|
||||
template <>
|
||||
inline boost::multiprecision::cpp_bin_float_quad real_cast<boost::multiprecision::cpp_bin_float_quad, mp_type>(mp_type val)
|
||||
{
|
||||
return boost::multiprecision::cpp_bin_float_quad(val);
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#else
|
||||
|
||||
#include <boost/multiprecision/mpfr.hpp>
|
||||
|
||||
typedef boost::multiprecision::mpfr_float mp_type;
|
||||
|
||||
inline void set_working_precision(int n)
|
||||
{
|
||||
boost::multiprecision::mpfr_float::default_precision(boost::multiprecision::detail::digits2_2_10(n));
|
||||
}
|
||||
|
||||
inline void set_output_precision(int n)
|
||||
{
|
||||
std::cout << std::setprecision(n);
|
||||
std::cerr << std::setprecision(n);
|
||||
}
|
||||
|
||||
inline mp_type round_to_precision(mp_type m, int bits)
|
||||
{
|
||||
mpfr_prec_round(m.backend().data(), bits, MPFR_RNDN);
|
||||
return m;
|
||||
}
|
||||
|
||||
inline int get_working_precision()
|
||||
{
|
||||
return mp_type::default_precision();
|
||||
}
|
||||
|
||||
namespace boost {
|
||||
namespace math {
|
||||
namespace tools {
|
||||
|
||||
template <>
|
||||
inline boost::multiprecision::cpp_bin_float_double_extended real_cast<boost::multiprecision::cpp_bin_float_double_extended, mp_type>(mp_type val)
|
||||
{
|
||||
return boost::multiprecision::cpp_bin_float_double_extended(val);
|
||||
}
|
||||
template <>
|
||||
inline boost::multiprecision::cpp_bin_float_quad real_cast<boost::multiprecision::cpp_bin_float_quad, mp_type>(mp_type val)
|
||||
{
|
||||
return boost::multiprecision::cpp_bin_float_quad(val);
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
|
||||
#endif // BOOST_REMEZ_MULTIPRECISION_HPP
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user