boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 16:32:10 +00:00
Code Issues Projects Releases Wiki Activity

Remove all references to tr1 components.

Browse Source 
Update data generators to use Boost.Multiprecision.
This commit is contained in:
jzmaddock
2014-09-28 18:09:27 +01:00
parent bd807c74e1
commit 9860071f84
36 changed files with 536 additions and 599 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
test/test_carlson.hpp
View File

@@ -11,7 +11,7 @@
#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/array.hpp>
#include <boost/tr1/random.hpp>
#include <boost/random.hpp>
#include "functor.hpp"
#include "handle_test_result.hpp"
@@ -150,8 +150,8 @@ void test_spots(T, const char* type_name)
// Sanity/consistency checks from Numerical Computation of Real or Complex
// Elliptic Integrals, B. C. Carlson: http://arxiv.org/abs/math.CA/9409227
std::tr1::mt19937 ran;
std::tr1::uniform_real<float> ur(0, 1000);
boost::mt19937 ran;
boost::uniform_real<float> ur(0, 1000);
T eps40 = 40 * tools::epsilon<T>();
for(unsigned i = 0; i < 1000; ++i)

14
tools/Jamfile.v2
View File

@@ -10,8 +10,6 @@
import modules ;
import path ;
local ntl-path = [ modules.peek : NTL_PATH ] ;
project
: requirements
<toolset>gcc:<cxxflags>-Wno-missing-braces
@@ -34,22 +32,12 @@ project
<define>BOOST_UBLAS_UNSUPPORTED_COMPILER=0
<include>.
<include>../include_private
<include>$(ntl-path)/include
;
if $(ntl-path)
{
lib ntl : [ GLOB $(ntl-path)/src : *.cpp ] ;
}
else
{
lib ntl ;
}
for local source in [ glob *_data.cpp ] generate_test_values.cpp igamma_temme_large_coef.cpp lanczos_generator.cpp factorial_tables.cpp generate_rational_test.cpp
{
exe $(source:B) : $(source) ntl ;
exe $(source:B) : $(source) ;
install $(source:B)_bin : $(source:B) : <location>bin ;
}

16
tools/bessel_data.cpp
View File

@@ -5,17 +5,16 @@
//
// Computes test data for the various bessel functions using
// archived - deliberately naive - version of the code.
// We'll rely on the high precision of boost::math::ntl::RR to get us out of
// We'll rely on the high precision of mp_t to get us out of
// trouble and not worry about how long the calculations take.
// This provides a reasonably independent set of test data to
// compare against newly added asymptotic expansions etc.
//
#include <fstream>
#include <boost/math/bindings/rr.hpp>
#include <boost/math/tools/test_data.hpp>
#include <boost/math/special_functions/bessel.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
using namespace boost::math;
@@ -267,11 +266,10 @@ int main(int argc, char* argv[])
std::cout << sph_bessel_j_bare(0., 0.1185395751953125e4) << std::endl;
std::cout << sph_bessel_j_bare(22., 0.6540834903717041015625) << std::endl;
parameter_info<boost::math::ntl::RR> arg1, arg2;
test_data<boost::math::ntl::RR> data;
std::cout << std::setprecision(40) << std::scientific;
boost::math::ntl::RR::SetPrecision(1000);
boost::math::ntl::RR::SetOutputPrecision(40);
parameter_info<mp_t> arg1, arg2;
test_data<mp_t> data;
int functype = 0;
std::string letter = "J";
@@ -315,7 +313,7 @@ int main(int argc, char* argv[])
do{
get_user_parameter_info(arg1, "v");
get_user_parameter_info(arg2, "x");
boost::math::ntl::RR (*fp)(boost::math::ntl::RR, boost::math::ntl::RR);
mp_t (*fp)(mp_t, mp_t);
if(functype == func_J)
fp = cyl_bessel_j_bare;
else if(functype == func_I)
@@ -346,7 +344,7 @@ int main(int argc, char* argv[])
line = "bessel_j_data.ipp";
std::ofstream ofs(line.c_str());
line.erase(line.find('.'));
ofs << std::scientific;
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, line.c_str());
return 0;

14
tools/beta_data.cpp
View File

@@ -3,22 +3,22 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/tools/test_data.hpp>
#include <fstream>
#include "mp_t.hpp"
using namespace boost::math::tools;
struct beta_data_generator
{
boost::math::ntl::RR operator()(boost::math::ntl::RR a, boost::math::ntl::RR b)
mp_t operator()(mp_t a, mp_t b)
{
if(a < b)
throw std::domain_error("");
// very naively calculate spots:
boost::math::ntl::RR g1, g2, g3;
mp_t g1, g2, g3;
int s1, s2, s3;
g1 = boost::math::lgamma(a, &s1);
g2 = boost::math::lgamma(b, &s2);
@@ -33,11 +33,8 @@ struct beta_data_generator
int main()
{
boost::math::ntl::RR::SetPrecision(1000);
boost::math::ntl::RR::SetOutputPrecision(40);
parameter_info<boost::math::ntl::RR> arg1, arg2;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1, arg2;
test_data<mp_t> data;
std::cout << "Welcome.\n"
"This program will generate spot tests for the beta function:\n"
@@ -63,6 +60,7 @@ int main()
if(line == "")
line = "beta_data.ipp";
std::ofstream ofs(line.c_str());
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, "beta_data");
return 0;

82
tools/carlson_ellint_data.cpp
View File

@@ -4,14 +4,13 @@
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
//#include <boost/math/tools/dual_precision.hpp>
#include <boost/math/tools/test_data.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include <boost/math/special_functions/ellint_3.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include <boost/tr1/random.hpp>
#include <boost/random.hpp>
#include "mp_t.hpp"
float extern_val;
// confuse the compilers optimiser, and force a truncation to float precision:
@@ -21,82 +20,82 @@ float truncate_to_float(float const * pf)
return *pf;
}
boost::math::tuple<boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR> generate_rf_data(boost::math::ntl::RR n)
boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rf_data(mp_t n)
{
static std::tr1::mt19937 r;
std::tr1::uniform_real<float> ur(0, 1);
std::tr1::uniform_int<int> ui(-100, 100);
static boost::mt19937 r;
boost::uniform_real<float> ur(0, 1);
boost::uniform_int<int> ui(-100, 100);
float x = ur(r);
x = ldexp(x, ui(r));
boost::math::ntl::RR xr(truncate_to_float(&x));
mp_t xr(truncate_to_float(&x));
float y = ur(r);
y = ldexp(y, ui(r));
boost::math::ntl::RR yr(truncate_to_float(&y));
mp_t yr(truncate_to_float(&y));
float z = ur(r);
z = ldexp(z, ui(r));
boost::math::ntl::RR zr(truncate_to_float(&z));
mp_t zr(truncate_to_float(&z));
boost::math::ntl::RR result = boost::math::ellint_rf(xr, yr, zr);
mp_t result = boost::math::ellint_rf(xr, yr, zr);
return boost::math::make_tuple(xr, yr, zr, result);
}
boost::math::tuple<boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR> generate_rc_data(boost::math::ntl::RR n)
boost::math::tuple<mp_t, mp_t, mp_t> generate_rc_data(mp_t n)
{
static std::tr1::mt19937 r;
std::tr1::uniform_real<float> ur(0, 1);
std::tr1::uniform_int<int> ui(-100, 100);
static boost::mt19937 r;
boost::uniform_real<float> ur(0, 1);
boost::uniform_int<int> ui(-100, 100);
float x = ur(r);
x = ldexp(x, ui(r));
boost::math::ntl::RR xr(truncate_to_float(&x));
mp_t xr(truncate_to_float(&x));
float y = ur(r);
y = ldexp(y, ui(r));
boost::math::ntl::RR yr(truncate_to_float(&y));
mp_t yr(truncate_to_float(&y));
boost::math::ntl::RR result = boost::math::ellint_rc(xr, yr);
mp_t result = boost::math::ellint_rc(xr, yr);
return boost::math::make_tuple(xr, yr, result);
}
boost::math::tuple<boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR> generate_rj_data(boost::math::ntl::RR n)
boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t> generate_rj_data(mp_t n)
{
static std::tr1::mt19937 r;
std::tr1::uniform_real<float> ur(0, 1);
std::tr1::uniform_real<float> nur(-1, 1);
std::tr1::uniform_int<int> ui(-100, 100);
static boost::mt19937 r;
boost::uniform_real<float> ur(0, 1);
boost::uniform_real<float> nur(-1, 1);
boost::uniform_int<int> ui(-100, 100);
float x = ur(r);
x = ldexp(x, ui(r));
boost::math::ntl::RR xr(truncate_to_float(&x));
mp_t xr(truncate_to_float(&x));
float y = ur(r);
y = ldexp(y, ui(r));
boost::math::ntl::RR yr(truncate_to_float(&y));
mp_t yr(truncate_to_float(&y));
float z = ur(r);
z = ldexp(z, ui(r));
boost::math::ntl::RR zr(truncate_to_float(&z));
mp_t zr(truncate_to_float(&z));
float p = nur(r);
p = ldexp(p, ui(r));
boost::math::ntl::RR pr(truncate_to_float(&p));
mp_t pr(truncate_to_float(&p));
boost::math::ellint_rj(x, y, z, p);
boost::math::ntl::RR result = boost::math::ellint_rj(xr, yr, zr, pr);
mp_t result = boost::math::ellint_rj(xr, yr, zr, pr);
return boost::math::make_tuple(xr, yr, zr, pr, result);
}
boost::math::tuple<boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR> generate_rd_data(boost::math::ntl::RR n)
boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rd_data(mp_t n)
{
static std::tr1::mt19937 r;
std::tr1::uniform_real<float> ur(0, 1);
std::tr1::uniform_int<int> ui(-100, 100);
static boost::mt19937 r;
boost::uniform_real<float> ur(0, 1);
boost::uniform_int<int> ui(-100, 100);
float x = ur(r);
x = ldexp(x, ui(r));
boost::math::ntl::RR xr(truncate_to_float(&x));
mp_t xr(truncate_to_float(&x));
float y = ur(r);
y = ldexp(y, ui(r));
boost::math::ntl::RR yr(truncate_to_float(&y));
mp_t yr(truncate_to_float(&y));
float z = ur(r);
z = ldexp(z, ui(r));
boost::math::ntl::RR zr(truncate_to_float(&z));
mp_t zr(truncate_to_float(&z));
boost::math::ntl::RR result = boost::math::ellint_rd(xr, yr, zr);
mp_t result = boost::math::ellint_rd(xr, yr, zr);
return boost::math::make_tuple(xr, yr, zr, result);
}
@@ -104,11 +103,8 @@ int cpp_main(int argc, char*argv [])
{
using namespace boost::math::tools;
boost::math::ntl::RR::SetOutputPrecision(50);
boost::math::ntl::RR::SetPrecision(1000);
parameter_info<boost::math::ntl::RR> arg1, arg2;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1, arg2;
test_data<mp_t> data;
bool cont;
std::string line;
@@ -121,7 +117,7 @@ int cpp_main(int argc, char*argv [])
std::cout << "Number of points: ";
std::cin >> count;
arg1 = make_periodic_param(boost::math::ntl::RR(0), boost::math::ntl::RR(1), count);
arg1 = make_periodic_param(mp_t(0), mp_t(1), count);
arg1.type |= dummy_param;
//
@@ -142,7 +138,7 @@ int cpp_main(int argc, char*argv [])
line = "ellint_rf_data.ipp";
std::ofstream ofs(line.c_str());
line.erase(line.find('.'));
ofs << std::scientific;
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, line.c_str());
return 0;

16
tools/cbrt_data.cpp
View File

@@ -3,19 +3,18 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
using namespace std;
struct cube_data_generator
{
boost::math::ntl::RR operator()(boost::math::ntl::RR z)
mp_t operator()(mp_t z)
{
boost::math::ntl::RR result = z*z*z;
mp_t result = z*z*z;
// if result is out of range of a float,
// don't include in test data as it messes up our results:
if(result > (std::numeric_limits<float>::max)())
@@ -28,11 +27,8 @@ struct cube_data_generator
int main(int argc, char* argv[])
{
boost::math::ntl::RR::SetPrecision(1000);
boost::math::ntl::RR::SetOutputPrecision(40);
parameter_info<boost::math::ntl::RR> arg1;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1;
test_data<mp_t> data;
std::cout << "Welcome.\n"
"This program will generate spot tests for the cbrt function:\n\n";
@@ -57,7 +53,7 @@ int main(int argc, char* argv[])
if(line == "")
line = "cbrt_data.ipp";
std::ofstream ofs(line.c_str());
ofs << std::scientific;
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, "cbrt_data");
return 0;

15
tools/digamma_data.cpp
View File

@@ -3,12 +3,11 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/math/special_functions/digamma.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
using namespace std;
@@ -20,7 +19,7 @@ float force_truncate(const float* f)
return external_f;
}
float truncate_to_float(boost::math::ntl::RR r)
float truncate_to_float(mp_t r)
{
float f = boost::math::tools::real_cast<float>(r);
return force_truncate(&f);
@@ -28,11 +27,8 @@ float truncate_to_float(boost::math::ntl::RR r)
int cpp_main(int argc, char*argv [])
{
boost::math::ntl::RR::SetPrecision(1000);
boost::math::ntl::RR::SetOutputPrecision(40);
parameter_info<boost::math::ntl::RR> arg1;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1;
test_data<mp_t> data;
bool cont;
std::string line;
@@ -44,7 +40,7 @@ int cpp_main(int argc, char*argv [])
do{
if(0 == get_user_parameter_info(arg1, "z"))
return 1;
data.insert(&boost::math::digamma<boost::math::ntl::RR>, arg1);
data.insert(&boost::math::digamma<mp_t>, arg1);
std::cout << "Any more data [y/n]?";
std::getline(std::cin, line);
@@ -58,6 +54,7 @@ int cpp_main(int argc, char*argv [])
if(line == "")
line = "digamma_data.ipp";
std::ofstream ofs(line.c_str());
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, "digamma_data");
return 0;

14
tools/ellint_e_data.cpp
View File

@@ -3,13 +3,12 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/math/tools/test_data.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include <boost/math/special_functions/ellint_2.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include <boost/tr1/random.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
using namespace boost::math;
@@ -25,11 +24,8 @@ int cpp_main(int argc, char*argv [])
{
using namespace boost::math::tools;
boost::math::ntl::RR::SetOutputPrecision(50);
boost::math::ntl::RR::SetPrecision(1000);
parameter_info<boost::math::ntl::RR> arg1;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1;
test_data<mp_t> data;
bool cont;
std::string line;
@@ -41,7 +37,7 @@ int cpp_main(int argc, char*argv [])
if(0 == get_user_parameter_info(arg1, "phi"))
return 1;
data.insert(&ellint_e_data<boost::math::ntl::RR>, arg1);
data.insert(&ellint_e_data<mp_t>, arg1);
std::cout << "Any more data [y/n]?";
std::getline(std::cin, line);
@@ -56,7 +52,7 @@ int cpp_main(int argc, char*argv [])
line = "ellint_e_data.ipp";
std::ofstream ofs(line.c_str());
line.erase(line.find('.'));
ofs << std::scientific;
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, line.c_str());
return 0;

14
tools/ellint_f_data.cpp
View File

@@ -3,13 +3,12 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/math/tools/test_data.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include <boost/math/special_functions/ellint_1.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include <boost/tr1/random.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
using namespace boost::math;
@@ -35,11 +34,8 @@ int cpp_main(int argc, char*argv [])
{
using namespace boost::math::tools;
boost::math::ntl::RR::SetOutputPrecision(50);
boost::math::ntl::RR::SetPrecision(1000);
parameter_info<boost::math::ntl::RR> arg1, arg2;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1, arg2;
test_data<mp_t> data;
bool cont;
std::string line;
@@ -53,7 +49,7 @@ int cpp_main(int argc, char*argv [])
if(0 == get_user_parameter_info(arg2, "k"))
return 1;
data.insert(&ellint_f_data<boost::math::ntl::RR>, arg1, arg2);
data.insert(&ellint_f_data<mp_t>, arg1, arg2);
std::cout << "Any more data [y/n]?";
std::getline(std::cin, line);
@@ -68,7 +64,7 @@ int cpp_main(int argc, char*argv [])
line = "ellint_f.ipp";
std::ofstream ofs(line.c_str());
line.erase(line.find('.'));
ofs << std::scientific;
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, line.c_str());
return 0;

14
tools/ellint_k_data.cpp
View File

@@ -3,13 +3,12 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/math/tools/test_data.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include <boost/math/special_functions/ellint_1.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include <boost/tr1/random.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
using namespace boost::math;
@@ -25,11 +24,8 @@ int cpp_main(int argc, char*argv [])
{
using namespace boost::math::tools;
boost::math::ntl::RR::SetOutputPrecision(50);
boost::math::ntl::RR::SetPrecision(1000);
parameter_info<boost::math::ntl::RR> arg1;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1;
test_data<mp_t> data;
bool cont;
std::string line;
@@ -41,7 +37,7 @@ int cpp_main(int argc, char*argv [])
if(0 == get_user_parameter_info(arg1, "phi"))
return 1;
data.insert(&ellint_k_data<boost::math::ntl::RR>, arg1);
data.insert(&ellint_k_data<mp_t>, arg1);
std::cout << "Any more data [y/n]?";
std::getline(std::cin, line);
@@ -56,7 +52,7 @@ int cpp_main(int argc, char*argv [])
line = "ellint_k_data.ipp";
std::ofstream ofs(line.c_str());
line.erase(line.find('.'));
ofs << std::scientific;
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, line.c_str());
return 0;

15
tools/ellint_pi2_data.cpp
View File

@@ -3,15 +3,13 @@
// Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
//#include <boost/math/tools/dual_precision.hpp>
#include <boost/math/tools/test_data.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include <boost/math/special_functions/ellint_2.hpp>
#include <boost/math/special_functions/ellint_3.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include <boost/tr1/random.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
using namespace boost::math;
@@ -21,11 +19,8 @@ int cpp_main(int argc, char*argv [])
{
using namespace boost::math::tools;
boost::math::ntl::RR::SetOutputPrecision(50);
boost::math::ntl::RR::SetPrecision(1000);
parameter_info<boost::math::ntl::RR> arg1, arg2;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1, arg2;
test_data<mp_t> data;
bool cont;
std::string line;
@@ -39,7 +34,7 @@ int cpp_main(int argc, char*argv [])
if(0 == get_user_parameter_info(arg2, "k"))
return 1;
boost::math::ntl::RR (*fp)(boost::math::ntl::RR, boost::math::ntl::RR) = &ellint_3;
mp_t (*fp)(mp_t, mp_t) = &ellint_3;
data.insert(fp, arg2, arg1);
std::cout << "Any more data [y/n]?";
@@ -55,7 +50,7 @@ int cpp_main(int argc, char*argv [])
line = "ellint_pi2_data.ipp";
std::ofstream ofs(line.c_str());
line.erase(line.find('.'));
ofs << std::scientific;
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, line.c_str());
return 0;

24
tools/ellint_pi3_data.cpp
View File

@@ -6,14 +6,13 @@
// or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
//#include <boost/math/tools/dual_precision.hpp>
#include <boost/math/tools/test_data.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include <boost/math/special_functions/ellint_3.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include <boost/tr1/random.hpp>
#include <boost/random.hpp>
#include "mp_t.hpp"
float extern_val;
// confuse the compilers optimiser, and force a truncation to float precision:
@@ -23,13 +22,13 @@ float truncate_to_float(float const * pf)
return *pf;
}
boost::math::tuple<boost::math::ntl::RR, boost::math::ntl::RR> generate_data(boost::math::ntl::RR n, boost::math::ntl::RR phi)
boost::math::tuple<mp_t, mp_t> generate_data(mp_t n, mp_t phi)
{
static std::tr1::mt19937 r;
std::tr1::uniform_real<float> ui(0, 1);
static boost::mt19937 r;
boost::uniform_real<float> ui(0, 1);
float k = ui(r);
boost::math::ntl::RR kr(truncate_to_float(&k));
boost::math::ntl::RR result = boost::math::ellint_3(kr, n, phi);
mp_t kr(truncate_to_float(&k));
mp_t result = boost::math::ellint_3(kr, n, phi);
return boost::math::make_tuple(kr, result);
}
@@ -37,11 +36,8 @@ int cpp_main(int argc, char*argv [])
{
using namespace boost::math::tools;
boost::math::ntl::RR::SetOutputPrecision(50);
boost::math::ntl::RR::SetPrecision(1000);
parameter_info<boost::math::ntl::RR> arg1, arg2;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1, arg2;
test_data<mp_t> data;
bool cont;
std::string line;
@@ -70,7 +66,7 @@ int cpp_main(int argc, char*argv [])
line = "ellint_pi3_data.ipp";
std::ofstream ofs(line.c_str());
line.erase(line.find('.'));
ofs << std::scientific;
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, line.c_str());
return 0;

45
tools/erf_data.cpp
View File

@@ -3,15 +3,12 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/special_functions/erf.hpp> // for inverses
#include <boost/math/constants/constants.hpp>
#include <boost/math/tools/test.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
using namespace std;
@@ -23,7 +20,7 @@ float force_truncate(const float* f)
return external_f;
}
float truncate_to_float(boost::math::ntl::RR r)
float truncate_to_float(mp_t r)
{
float f = boost::math::tools::real_cast<float>(r);
return force_truncate(&f);
@@ -31,7 +28,7 @@ float truncate_to_float(boost::math::ntl::RR r)
struct erf_data_generator
{
boost::math::tuple<boost::math::ntl::RR, boost::math::ntl::RR> operator()(boost::math::ntl::RR z)
boost::math::tuple<mp_t, mp_t> operator()(mp_t z)
{
// very naively calculate spots using the gamma function at high precision:
int sign = 1;
@@ -40,9 +37,9 @@ struct erf_data_generator
sign = -1;
z = -z;
}
boost::math::ntl::RR g1, g2;
g1 = boost::math::tgamma_lower(boost::math::ntl::RR(0.5), z * z);
g1 /= sqrt(boost::math::constants::pi<boost::math::ntl::RR>());
mp_t g1, g2;
g1 = boost::math::tgamma_lower(mp_t(0.5), z * z);
g1 /= sqrt(boost::math::constants::pi<mp_t>());
g1 *= sign;
if(z < 0.5)
@@ -51,8 +48,8 @@ struct erf_data_generator
}
else
{
g2 = boost::math::tgamma(boost::math::ntl::RR(0.5), z * z);
g2 /= sqrt(boost::math::constants::pi<boost::math::ntl::RR>());
g2 = boost::math::tgamma(mp_t(0.5), z * z);
g2 /= sqrt(boost::math::constants::pi<mp_t>());
}
if(sign < 1)
g2 = 2 - g2;
@@ -117,24 +114,21 @@ void asymptotic_limit(int Bits)
<< terms << " terms." << std::endl;
}
boost::math::tuple<boost::math::ntl::RR, boost::math::ntl::RR> erfc_inv(boost::math::ntl::RR r)
boost::math::tuple<mp_t, mp_t> erfc_inv(mp_t r)
{
boost::math::ntl::RR x = exp(-r * r);
x = NTL::RoundToPrecision(x.value(), 64);
mp_t x = exp(-r * r);
x = x.convert_to<double>();
std::cout << x << " ";
boost::math::ntl::RR result = boost::math::erfc_inv(x);
mp_t result = boost::math::erfc_inv(x);
std::cout << result << std::endl;
return boost::math::make_tuple(x, result);
}
int cpp_main(int argc, char*argv [])
int main(int argc, char*argv [])
{
boost::math::ntl::RR::SetPrecision(1000);
boost::math::ntl::RR::SetOutputPrecision(40);
parameter_info<boost::math::ntl::RR> arg1;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1;
test_data<mp_t> data;
bool cont;
std::string line;
@@ -152,18 +146,18 @@ int cpp_main(int argc, char*argv [])
}
else if(strcmp(argv[1], "--erf_inv") == 0)
{
boost::math::ntl::RR (*f)(boost::math::ntl::RR);
mp_t (*f)(mp_t);
f = boost::math::erf_inv;
std::cout << "Welcome.\n"
"This program will generate spot tests for the inverse erf function:\n";
std::cout << "Enter the number of data points: ";
int points;
std::cin >> points;
data.insert(f, make_random_param(boost::math::ntl::RR(-1), boost::math::ntl::RR(1), points));
data.insert(f, make_random_param(mp_t(-1), mp_t(1), points));
}
else if(strcmp(argv[1], "--erfc_inv") == 0)
{
boost::math::tuple<boost::math::ntl::RR, boost::math::ntl::RR> (*f)(boost::math::ntl::RR);
boost::math::tuple<mp_t, mp_t> (*f)(mp_t);
f = erfc_inv;
std::cout << "Welcome.\n"
"This program will generate spot tests for the inverse erfc function:\n";
@@ -173,7 +167,7 @@ int cpp_main(int argc, char*argv [])
std::cout << "Enter the number of data points: ";
int points;
std::cin >> points;
parameter_info<boost::math::ntl::RR> arg = make_random_param(boost::math::ntl::RR(0), boost::math::ntl::RR(max_val), points);
parameter_info<mp_t> arg = make_random_param(mp_t(0), mp_t(max_val), points);
arg.type |= dummy_param;
data.insert(f, arg);
}
@@ -202,6 +196,7 @@ int cpp_main(int argc, char*argv [])
if(line == "")
line = "erf_data.ipp";
std::ofstream ofs(line.c_str());
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, "erf_data");
return 0;

17
tools/expint_data.cpp
View File

@@ -3,22 +3,21 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/math/special_functions/expint.hpp>
#include <boost/math/constants/constants.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
struct expint_data_generator
{
boost::math::ntl::RR operator()(boost::math::ntl::RR a, boost::math::ntl::RR b)
mp_t operator()(mp_t a, mp_t b)
{
unsigned n = boost::math::tools::real_cast<unsigned>(a);
std::cout << n << " " << b << " ";
boost::math::ntl::RR result = boost::math::expint(n, b);
mp_t result = boost::math::expint(n, b);
std::cout << result << std::endl;
return result;
}
@@ -27,14 +26,11 @@ struct expint_data_generator
int main()
{
boost::math::ntl::RR::SetPrecision(1000);
boost::math::ntl::RR::SetOutputPrecision(40);
boost::math::expint(1, 0.06227754056453704833984375);
std::cout << boost::math::expint(1, boost::math::ntl::RR(0.5)) << std::endl;
std::cout << boost::math::expint(1, mp_t(0.5)) << std::endl;
parameter_info<boost::math::ntl::RR> arg1, arg2;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1, arg2;
test_data<mp_t> data;
std::cout << "Welcome.\n"
"This program will generate spot tests for the expint function:\n"
@@ -60,6 +56,7 @@ int main()
if(line == "")
line = "expint_data.ipp";
std::ofstream ofs(line.c_str());
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, "expint_data");
return 0;

13
tools/expint_i_data.cpp
View File

@@ -3,25 +3,21 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/math/special_functions/expint.hpp>
#include <boost/math/constants/constants.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
int main()
{
boost::math::ntl::RR::SetPrecision(1000);
boost::math::ntl::RR::SetOutputPrecision(40);
parameter_info<mp_t> arg1;
test_data<mp_t> data;
parameter_info<boost::math::ntl::RR> arg1;
test_data<boost::math::ntl::RR> data;
boost::math::ntl::RR (*f)(boost::math::ntl::RR) = boost::math::expint;
mp_t (*f)(mp_t) = boost::math::expint;
std::cout << "Welcome.\n"
"This program will generate spot tests for the expint Ei function:\n"
@@ -46,6 +42,7 @@ int main()
if(line == "")
line = "expinti_data.ipp";
std::ofstream ofs(line.c_str());
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, "expinti_data");
return 0;

12
tools/factorial_tables.cpp
View File

@@ -3,18 +3,18 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/limits.hpp>
#include <vector>
#include "mp_t.hpp"
void write_table(unsigned max_exponent)
{
boost::math::ntl::RR max = ldexp(boost::math::ntl::RR(1), max_exponent);
mp_t max = ldexp(mp_t(1), (int)max_exponent);
std::vector<boost::math::ntl::RR> factorials;
std::vector<mp_t> factorials;
factorials.push_back(1);
boost::math::ntl::RR f(1);
mp_t f(1);
unsigned i = 1;
while(f < max)
@@ -27,7 +27,7 @@ void write_table(unsigned max_exponent)
//
// now write out the results to cout:
//
std::cout << std::scientific;
std::cout << std::scientific << std::setprecision(40);
std::cout << " static const boost::array<T, " << factorials.size() << "> factorials = {\n";
for(unsigned j = 0; j < factorials.size(); ++j)
std::cout << " " << factorials[j] << "L,\n";
@@ -37,7 +37,5 @@ void write_table(unsigned max_exponent)
int main()
{
boost::math::ntl::RR::SetPrecision(300);
boost::math::ntl::RR::SetOutputPrecision(40);
write_table(16384/*std::numeric_limits<float>::max_exponent*/);
}

27
tools/gamma_P_inva_data.cpp
View File

@@ -3,15 +3,12 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/tools/test.hpp>
#include <boost/lexical_cast.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
@@ -28,7 +25,7 @@ float force_truncate(const float* f)
return external_f;
}
float truncate_to_float(boost::math::ntl::RR r)
float truncate_to_float(mp_t r)
{
float f = boost::math::tools::real_cast<float>(r);
return force_truncate(&f);
@@ -36,37 +33,35 @@ float truncate_to_float(boost::math::ntl::RR r)
struct gamma_inverse_generator_a
{
boost::math::tuple<boost::math::ntl::RR, boost::math::ntl::RR> operator()(const boost::math::ntl::RR x, const boost::math::ntl::RR p)
boost::math::tuple<mp_t, mp_t> operator()(const mp_t x, const mp_t p)
{
boost::math::ntl::RR x1 = boost::math::gamma_p_inva(x, p);
boost::math::ntl::RR x2 = boost::math::gamma_q_inva(x, p);
mp_t x1 = boost::math::gamma_p_inva(x, p);
mp_t x2 = boost::math::gamma_q_inva(x, p);
std::cout << "Inverse for " << x << " " << p << std::endl;
return boost::math::make_tuple(x1, x2);
}
};
int cpp_main(int argc, char*argv [])
int main(int argc, char*argv [])
{
boost::math::ntl::RR::SetPrecision(1000);
boost::math::ntl::RR::SetOutputPrecision(100);
bool cont;
std::string line;
parameter_info<boost::math::ntl::RR> arg1, arg2;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1, arg2;
test_data<mp_t> data;
std::cout << "Welcome.\n"
"This program will generate spot tests for the inverse incomplete gamma function:\n"
" gamma_p_inva(a, p) and gamma_q_inva(a, q)\n\n";
arg1 = make_power_param<boost::math::ntl::RR>(boost::math::ntl::RR(0), -4, 24);
arg2 = make_random_param<boost::math::ntl::RR>(boost::math::ntl::RR(0), boost::math::ntl::RR(1), 15);
arg1 = make_power_param<mp_t>(mp_t(0), -4, 24);
arg2 = make_random_param<mp_t>(mp_t(0), mp_t(1), 15);
data.insert(gamma_inverse_generator_a(), arg1, arg2);
line = "igamma_inva_data.ipp";
std::ofstream ofs(line.c_str());
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, "igamma_inva_data");
return 0;

76
tools/generate_rational_test.cpp
View File

@@ -6,11 +6,11 @@
#define BOOST_MATH_POLY_METHOD 0
#define BOOST_MATH_RATIONAL_METHOD 0
#include <boost/math/bindings/rr.hpp>
#include <boost/tr1/random.hpp>
#include <boost/random.hpp>
#include <boost/math/tools/rational.hpp>
#include <iostream>
#include <fstream>
#include "mp_t.hpp"
int main()
{
@@ -18,14 +18,12 @@ int main()
using namespace boost::math::tools;
static const unsigned max_order = 20;
std::cout << std::scientific << std::setprecision(40);
boost::math::ntl::RR::SetPrecision(500);
boost::math::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));
boost::mt19937 rnd;
boost::variate_generator<
boost::mt19937,
boost::uniform_int<> > gen(rnd, boost::uniform_int<>(1, 12));
for(unsigned i = 1; i < max_order; ++i)
{
@@ -57,12 +55,12 @@ int main()
}
std::cout << " };\n";
boost::math::ntl::RR r1 = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(0.125), i);
boost::math::ntl::RR r2 = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(0.25), i);
boost::math::ntl::RR r3 = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(0.75), i);
boost::math::ntl::RR r4 = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(1) - boost::math::ntl::RR(1) / 64, i);
boost::math::ntl::RR r5 = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(6.5), i);
boost::math::ntl::RR r6 = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(10247.25), i);
mp_t r1 = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(0.125), i);
mp_t r2 = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(0.25), i);
mp_t r3 = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(0.75), i);
mp_t r4 = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(1) - mp_t(1) / 64, i);
mp_t r5 = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(6.5), i);
mp_t r6 = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(10247.25), i);
std::cout <<
" BOOST_CHECK_CLOSE(\n"
@@ -163,12 +161,12 @@ int main()
" static_cast<T>(" << r6 << "L),\n"
" tolerance);\n\n";
r1 = boost::math::tools::evaluate_even_polynomial(&coef[0], boost::math::ntl::RR(0.125), i);
r2 = boost::math::tools::evaluate_even_polynomial(&coef[0], boost::math::ntl::RR(0.25), i);
r3 = boost::math::tools::evaluate_even_polynomial(&coef[0], boost::math::ntl::RR(0.75), i);
r4 = boost::math::tools::evaluate_even_polynomial(&coef[0], boost::math::ntl::RR(1) - boost::math::ntl::RR(1) / 64, i);
r5 = boost::math::tools::evaluate_even_polynomial(&coef[0], boost::math::ntl::RR(6.5), i);
r6 = boost::math::tools::evaluate_even_polynomial(&coef[0], boost::math::ntl::RR(10247.25), i);
r1 = boost::math::tools::evaluate_even_polynomial(&coef[0], mp_t(0.125), i);
r2 = boost::math::tools::evaluate_even_polynomial(&coef[0], mp_t(0.25), i);
r3 = boost::math::tools::evaluate_even_polynomial(&coef[0], mp_t(0.75), i);
r4 = boost::math::tools::evaluate_even_polynomial(&coef[0], mp_t(1) - mp_t(1) / 64, i);
r5 = boost::math::tools::evaluate_even_polynomial(&coef[0], mp_t(6.5), i);
r6 = boost::math::tools::evaluate_even_polynomial(&coef[0], mp_t(10247.25), i);
std::cout <<
" BOOST_CHECK_CLOSE(\n"
@@ -271,12 +269,12 @@ int main()
if(i > 1)
{
r1 = boost::math::tools::evaluate_odd_polynomial(&coef[0], boost::math::ntl::RR(0.125), i);
r2 = boost::math::tools::evaluate_odd_polynomial(&coef[0], boost::math::ntl::RR(0.25), i);
r3 = boost::math::tools::evaluate_odd_polynomial(&coef[0], boost::math::ntl::RR(0.75), i);
r4 = boost::math::tools::evaluate_odd_polynomial(&coef[0], boost::math::ntl::RR(1) - boost::math::ntl::RR(1) / 64, i);
r5 = boost::math::tools::evaluate_odd_polynomial(&coef[0], boost::math::ntl::RR(6.5), i);
r6 = boost::math::tools::evaluate_odd_polynomial(&coef[0], boost::math::ntl::RR(10247.25), i);
r1 = boost::math::tools::evaluate_odd_polynomial(&coef[0], mp_t(0.125), i);
r2 = boost::math::tools::evaluate_odd_polynomial(&coef[0], mp_t(0.25), i);
r3 = boost::math::tools::evaluate_odd_polynomial(&coef[0], mp_t(0.75), i);
r4 = boost::math::tools::evaluate_odd_polynomial(&coef[0], mp_t(1) - mp_t(1) / 64, i);
r5 = boost::math::tools::evaluate_odd_polynomial(&coef[0], mp_t(6.5), i);
r6 = boost::math::tools::evaluate_odd_polynomial(&coef[0], mp_t(10247.25), i);
std::cout <<
" BOOST_CHECK_CLOSE(\n"
@@ -378,12 +376,12 @@ int main()
" tolerance);\n\n";
}
r1 = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(0.125), i);
r2 = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(0.25), i);
r3 = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(0.75), i);
r4 = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(1) - boost::math::ntl::RR(1) / 64, i);
r5 = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(6.5), i);
r6 = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(10247.25), i);
r1 = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(0.125), i);
r2 = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(0.25), i);
r3 = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(0.75), i);
r4 = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(1) - mp_t(1) / 64, i);
r5 = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(6.5), i);
r6 = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(10247.25), i);
coef.clear();
for(unsigned j = 0; j < i; ++j)
@@ -412,12 +410,12 @@ int main()
}
std::cout << " };\n";
boost::math::ntl::RR r1d = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(0.125), i);
boost::math::ntl::RR r2d = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(0.25), i);
boost::math::ntl::RR r3d = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(0.75), i);
boost::math::ntl::RR r4d = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(1) - boost::math::ntl::RR(1) / 64, i);
boost::math::ntl::RR r5d = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(6.5), i);
boost::math::ntl::RR r6d = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(10247.25), i);
mp_t r1d = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(0.125), i);
mp_t r2d = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(0.25), i);
mp_t r3d = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(0.75), i);
mp_t r4d = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(1) - mp_t(1) / 64, i);
mp_t r5d = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(6.5), i);
mp_t r6d = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(10247.25), i);
std::cout <<
" BOOST_CHECK_CLOSE(\n"

26
tools/generate_test_values.cpp
View File

@@ -7,25 +7,33 @@
# define NTL_STD_CXX
#endif
#include <NTL/RR.h>
#include <iostream>
#include <iomanip>
#include "mp_t.hpp"
mp_t log1p(mp_t arg)
{
return log(arg + 1);
}
mp_t expm1(mp_t arg)
{
return exp(arg) - 1;
}
int main()
{
NTL::RR r, root_two;
r.SetPrecision(256);
root_two.SetPrecision(256);
mp_t r, root_two;
r = 1.0;
root_two = 2.0;
root_two = NTL::sqrt(root_two);
root_two = sqrt(root_two);
r /= root_two;
NTL::RR lim = NTL::pow(NTL::RR(NTL::INIT_VAL, 2.0), NTL::RR(NTL::INIT_VAL, -128));
NTL::RR::SetOutputPrecision(40);
mp_t lim = pow(mp_t(2), mp_t(-128));
std::cout << std::scientific << std::setprecision(40);
while(r > lim)
{
std::cout << " { " << r << "L, " << NTL::log1p(r) << "L, " << NTL::expm1(r) << "L, }, \n";
std::cout << " { " << -r << "L, " << NTL::log1p(-r) << "L, " << NTL::expm1(-r) << "L, }, \n";
std::cout << " { " << r << "L, " << log1p(r) << "L, " << expm1(r) << "L, }, \n";
std::cout << " { " << -r << "L, " << log1p(-r) << "L, " << expm1(-r) << "L, }, \n";
r /= root_two;
}
return 0;

20
tools/hermite_data.cpp
View File

@@ -3,13 +3,10 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/math/tools/test_data.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include <boost/math/special_functions/hermite.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include <boost/tr1/random.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
using namespace boost::math;
@@ -24,17 +21,14 @@ boost::math::tuple<T, T, T> hermite_data(T n, T x)
return boost::math::make_tuple(n, x, r1);
}
int cpp_main(int argc, char*argv [])
int main(int argc, char*argv [])
{
using namespace boost::math::tools;
boost::math::ntl::RR::SetOutputPrecision(50);
boost::math::ntl::RR::SetPrecision(1000);
parameter_info<mp_t> arg1, arg2, arg3;
test_data<mp_t> data;
parameter_info<boost::math::ntl::RR> arg1, arg2, arg3;
test_data<boost::math::ntl::RR> data;
std::cout << boost::math::hermite(10, static_cast<boost::math::ntl::RR>(1e300)) << std::endl;
std::cout << boost::math::hermite(10, static_cast<mp_t>(1e300)) << std::endl;
bool cont;
std::string line;
@@ -50,7 +44,7 @@ int cpp_main(int argc, char*argv [])
arg1.type |= dummy_param;
arg2.type |= dummy_param;
data.insert(&hermite_data<boost::math::ntl::RR>, arg1, arg2);
data.insert(&hermite_data<mp_t>, arg1, arg2);
std::cout << "Any more data [y/n]?";
std::getline(std::cin, line);
@@ -65,7 +59,7 @@ int cpp_main(int argc, char*argv [])
line = "hermite.ipp";
std::ofstream ofs(line.c_str());
line.erase(line.find('.'));
ofs << std::scientific;
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, line.c_str());
return 0;

45
tools/hypergeometric_dist_data.cpp
View File

@@ -5,32 +5,29 @@
//#define BOOST_MATH_INSTRUMENT
#include <boost/math/bindings/rr.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include <boost/math/distributions/hypergeometric.hpp>
#include <boost/math/special_functions/trunc.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/tools/test.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/random/uniform_int.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
std::tr1::mt19937 rnd;
std::mt19937 rnd;
struct hypergeometric_generator
{
boost::math::tuple<
boost::math::ntl::RR,
boost::math::ntl::RR,
boost::math::ntl::RR,
boost::math::ntl::RR,
boost::math::ntl::RR,
boost::math::ntl::RR,
boost::math::ntl::RR> operator()(boost::math::ntl::RR rN, boost::math::ntl::RR rr, boost::math::ntl::RR rn)
mp_t,
mp_t,
mp_t,
mp_t,
mp_t,
mp_t,
mp_t> operator()(mp_t rN, mp_t rr, mp_t rn)
{
using namespace std;
using namespace boost;
@@ -46,16 +43,16 @@ struct hypergeometric_generator
boost::uniform_int<> ui((std::max)(0, n + r - N), (std::min)(n, r));
int k = ui(rnd);
hypergeometric_distribution<ntl::RR> d(r, n, N);
hypergeometric_distribution<mp_t> d(r, n, N);
ntl::RR p = pdf(d, k);
mp_t p = pdf(d, k);
if((p == 1) || (p == 0))
{
// trivial case, don't clutter up our table with it:
throw std::domain_error("");
}
ntl::RR c = cdf(d, k);
ntl::RR cc = cdf(complement(d, k));
mp_t c = cdf(d, k);
mp_t cc = cdf(complement(d, k));
std::cout << "N = " << N << " r = " << r << " n = " << n << " PDF = " << p << " CDF = " << c << " CCDF = " << cc << std::endl;
@@ -69,21 +66,18 @@ struct hypergeometric_generator
}
};
int cpp_main(int argc, char*argv [])
int main(int argc, char*argv [])
{
boost::math::ntl::RR::SetPrecision(1000);
boost::math::ntl::RR::SetOutputPrecision(100);
std::string line;
parameter_info<boost::math::ntl::RR> arg1, arg2, arg3;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1, arg2, arg3;
test_data<mp_t> data;
std::cout << "Welcome.\n"
"This program will generate spot tests hypergeoemtric distribution:\n";
arg1 = make_power_param(boost::math::ntl::RR(0), 1, 21);
arg2 = make_power_param(boost::math::ntl::RR(0), 1, 21);
arg3 = make_power_param(boost::math::ntl::RR(0), 1, 21);
arg1 = make_power_param(mp_t(0), 1, 21);
arg2 = make_power_param(mp_t(0), 1, 21);
arg3 = make_power_param(mp_t(0), 1, 21);
arg1.type |= dummy_param;
arg2.type |= dummy_param;
@@ -93,6 +87,7 @@ int cpp_main(int argc, char*argv [])
line = "hypergeometric_dist_data2.ipp";
std::ofstream ofs(line.c_str());
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, "hypergeometric_dist_data2");
return 0;

90
tools/ibeta_data.cpp
View File

@@ -3,17 +3,15 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/special_functions/beta.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/tools/test.hpp>
#include <boost/lexical_cast.hpp>
#include <fstream>
#include <map>
#include <boost/math/tools/test_data.hpp>
#include <boost/random.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
using namespace boost::math;
@@ -163,28 +161,28 @@ float force_truncate(const float* f)
return external_f;
}
float truncate_to_float(boost::math::ntl::RR r)
float truncate_to_float(mp_t r)
{
float f = boost::math::tools::real_cast<float>(r);
return force_truncate(&f);
}
std::tr1::mt19937 rnd;
std::tr1::uniform_real<float> ur_a(1.0F, 5.0F);
std::tr1::variate_generator<std::tr1::mt19937, std::tr1::uniform_real<float> > gen(rnd, ur_a);
std::tr1::uniform_real<float> ur_a2(0.0F, 100.0F);
std::tr1::variate_generator<std::tr1::mt19937, std::tr1::uniform_real<float> > gen2(rnd, ur_a2);
boost::mt19937 rnd;
boost::uniform_real<float> ur_a(1.0F, 5.0F);
boost::variate_generator<boost::mt19937, boost::uniform_real<float> > gen(rnd, ur_a);
boost::uniform_real<float> ur_a2(0.0F, 100.0F);
boost::variate_generator<boost::mt19937, boost::uniform_real<float> > gen2(rnd, ur_a2);
struct beta_data_generator
{
boost::math::tuple<boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR> operator()(boost::math::ntl::RR ap, boost::math::ntl::RR bp, boost::math::ntl::RR x_)
boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t, mp_t, mp_t> operator()(mp_t ap, mp_t bp, mp_t x_)
{
float a = truncate_to_float(real_cast<float>(gen() * pow(boost::math::ntl::RR(10), ap)));
float b = truncate_to_float(real_cast<float>(gen() * pow(boost::math::ntl::RR(10), bp)));
float a = truncate_to_float(real_cast<float>(gen() * pow(mp_t(10), ap)));
float b = truncate_to_float(real_cast<float>(gen() * pow(mp_t(10), bp)));
float x = truncate_to_float(real_cast<float>(x_));
std::cout << a << " " << b << " " << x << std::endl;
std::pair<boost::math::ntl::RR, boost::math::ntl::RR> ib_full = ibeta_fraction1(boost::math::ntl::RR(a), boost::math::ntl::RR(b), boost::math::ntl::RR(x));
std::pair<boost::math::ntl::RR, boost::math::ntl::RR> ib_reg = ibeta_fraction1_regular(boost::math::ntl::RR(a), boost::math::ntl::RR(b), boost::math::ntl::RR(x));
std::pair<mp_t, mp_t> ib_full = ibeta_fraction1(mp_t(a), mp_t(b), mp_t(x));
std::pair<mp_t, mp_t> ib_reg = ibeta_fraction1_regular(mp_t(a), mp_t(b), mp_t(x));
return boost::math::make_tuple(a, b, x, ib_full.first, ib_full.second, ib_reg.first, ib_reg.second);
}
};
@@ -192,71 +190,68 @@ struct beta_data_generator
// medium sized values:
struct beta_data_generator_medium
{
boost::math::tuple<boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR> operator()(boost::math::ntl::RR x_)
boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t, mp_t, mp_t> operator()(mp_t x_)
{
boost::math::ntl::RR a = gen2();
boost::math::ntl::RR b = gen2();
boost::math::ntl::RR x = x_;
a = ConvPrec(a.value(), 22);
b = ConvPrec(b.value(), 22);
x = ConvPrec(x.value(), 22);
mp_t a = gen2();
mp_t b = gen2();
mp_t x = x_;
a = ConvPrec(a, 22);
b = ConvPrec(b, 22);
x = ConvPrec(x, 22);
std::cout << a << " " << b << " " << x << std::endl;
//boost::math::ntl::RR exp_beta = boost::math::beta(a, b, x);
std::pair<boost::math::ntl::RR, boost::math::ntl::RR> ib_full = ibeta_fraction1(boost::math::ntl::RR(a), boost::math::ntl::RR(b), boost::math::ntl::RR(x));
//mp_t exp_beta = boost::math::beta(a, b, x);
std::pair<mp_t, mp_t> ib_full = ibeta_fraction1(mp_t(a), mp_t(b), mp_t(x));
/*exp_beta = boost::math::tools::relative_error(ib_full.first, exp_beta);
if(exp_beta > 1e-40)
{
std::cout << exp_beta << std::endl;
}*/
std::pair<boost::math::ntl::RR, boost::math::ntl::RR> ib_reg = ibeta_fraction1_regular(boost::math::ntl::RR(a), boost::math::ntl::RR(b), boost::math::ntl::RR(x));
std::pair<mp_t, mp_t> ib_reg = ibeta_fraction1_regular(mp_t(a), mp_t(b), mp_t(x));
return boost::math::make_tuple(a, b, x, ib_full.first, ib_full.second, ib_reg.first, ib_reg.second);
}
};
struct beta_data_generator_small
{
boost::math::tuple<boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR> operator()(boost::math::ntl::RR x_)
boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t, mp_t, mp_t> operator()(mp_t x_)
{
float a = truncate_to_float(gen2()/10);
float b = truncate_to_float(gen2()/10);
float x = truncate_to_float(real_cast<float>(x_));
std::cout << a << " " << b << " " << x << std::endl;
std::pair<boost::math::ntl::RR, boost::math::ntl::RR> ib_full = ibeta_fraction1(boost::math::ntl::RR(a), boost::math::ntl::RR(b), boost::math::ntl::RR(x));
std::pair<boost::math::ntl::RR, boost::math::ntl::RR> ib_reg = ibeta_fraction1_regular(boost::math::ntl::RR(a), boost::math::ntl::RR(b), boost::math::ntl::RR(x));
std::pair<mp_t, mp_t> ib_full = ibeta_fraction1(mp_t(a), mp_t(b), mp_t(x));
std::pair<mp_t, mp_t> ib_reg = ibeta_fraction1_regular(mp_t(a), mp_t(b), mp_t(x));
return boost::math::make_tuple(a, b, x, ib_full.first, ib_full.second, ib_reg.first, ib_reg.second);
}
};
struct beta_data_generator_int
{
boost::math::tuple<boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR> operator()(boost::math::ntl::RR a, boost::math::ntl::RR b, boost::math::ntl::RR x_)
boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t, mp_t, mp_t> operator()(mp_t a, mp_t b, mp_t x_)
{
float x = truncate_to_float(real_cast<float>(x_));
std::cout << a << " " << b << " " << x << std::endl;
std::pair<boost::math::ntl::RR, boost::math::ntl::RR> ib_full = ibeta_fraction1(a, b, boost::math::ntl::RR(x));
std::pair<boost::math::ntl::RR, boost::math::ntl::RR> ib_reg = ibeta_fraction1_regular(a, b, boost::math::ntl::RR(x));
std::pair<mp_t, mp_t> ib_full = ibeta_fraction1(a, b, mp_t(x));
std::pair<mp_t, mp_t> ib_reg = ibeta_fraction1_regular(a, b, mp_t(x));
return boost::math::make_tuple(a, b, x, ib_full.first, ib_full.second, ib_reg.first, ib_reg.second);
}
};
int cpp_main(int, char* [])
int main(int, char* [])
{
boost::math::ntl::RR::SetPrecision(1000);
boost::math::ntl::RR::SetOutputPrecision(40);
parameter_info<boost::math::ntl::RR> arg1, arg2, arg3, arg4, arg5;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1, arg2, arg3, arg4, arg5;
test_data<mp_t> data;
std::cout << "Welcome.\n"
"This program will generate spot tests for the incomplete beta functions:\n"
" beta(a, b, x) and ibeta(a, b, x)\n\n"
"This is not an interactive program be prepared for a long wait!!!\n\n";
arg1 = make_periodic_param(boost::math::ntl::RR(-5), boost::math::ntl::RR(6), 11);
arg2 = make_periodic_param(boost::math::ntl::RR(-5), boost::math::ntl::RR(6), 11);
arg3 = make_random_param(boost::math::ntl::RR(0.0001), boost::math::ntl::RR(1), 10);
arg4 = make_random_param(boost::math::ntl::RR(0.0001), boost::math::ntl::RR(1), 100 /*500*/);
arg5 = make_periodic_param(boost::math::ntl::RR(1), boost::math::ntl::RR(41), 10);
arg1 = make_periodic_param(mp_t(-5), mp_t(6), 11);
arg2 = make_periodic_param(mp_t(-5), mp_t(6), 11);
arg3 = make_random_param(mp_t(0.0001), mp_t(1), 10);
arg4 = make_random_param(mp_t(0.0001), mp_t(1), 100 /*500*/);
arg5 = make_periodic_param(mp_t(1), mp_t(41), 10);
arg1.type |= dummy_param;
arg2.type |= dummy_param;
@@ -271,16 +266,16 @@ int cpp_main(int, char* [])
//data.insert(beta_data_generator_small(), arg4);
data.insert(beta_data_generator_int(), arg5, arg5, arg3);
test_data<boost::math::ntl::RR>::const_iterator i, j;
test_data<mp_t>::const_iterator i, j;
i = data.begin();
j = data.end();
while(i != j)
{
boost::math::ntl::RR v1 = beta((*i)[0], (*i)[1], (*i)[2]);
boost::math::ntl::RR v2 = relative_error(v1, (*i)[3]);
mp_t v1 = beta((*i)[0], (*i)[1], (*i)[2]);
mp_t v2 = relative_error(v1, (*i)[3]);
std::string s = boost::lexical_cast<std::string>((*i)[3]);
boost::math::ntl::RR v3 = boost::lexical_cast<boost::math::ntl::RR>(s);
boost::math::ntl::RR v4 = relative_error(v3, (*i)[3]);
mp_t v3 = boost::lexical_cast<mp_t>(s);
mp_t v4 = relative_error(v3, (*i)[3]);
if(v2 > 1e-40)
{
std::cout << v2 << std::endl;
@@ -299,6 +294,7 @@ int cpp_main(int, char* [])
if(line == "")
line = "ibeta_data.ipp";
std::ofstream ofs(line.c_str());
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, "ibeta_data");
return 0;

46
tools/ibeta_inv_data.cpp
View File

@@ -3,16 +3,13 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#if 1
#include <boost/math/bindings/rr.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include <boost/math/special_functions/beta.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/tools/test.hpp>
#include <boost/lexical_cast.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include <boost/random.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
@@ -29,52 +26,49 @@ float force_truncate(const float* f)
return external_f;
}
float truncate_to_float(boost::math::ntl::RR r)
float truncate_to_float(mp_t r)
{
float f = boost::math::tools::real_cast<float>(r);
return force_truncate(&f);
}
std::tr1::mt19937 rnd;
std::tr1::uniform_real<float> ur_a(1.0F, 5.0F);
std::tr1::variate_generator<std::tr1::mt19937, std::tr1::uniform_real<float> > gen(rnd, ur_a);
std::tr1::uniform_real<float> ur_a2(0.0F, 100.0F);
std::tr1::variate_generator<std::tr1::mt19937, std::tr1::uniform_real<float> > gen2(rnd, ur_a2);
boost::mt19937 rnd;
boost::uniform_real<float> ur_a(1.0F, 5.0F);
boost::variate_generator<boost::mt19937, boost::uniform_real<float> > gen(rnd, ur_a);
boost::uniform_real<float> ur_a2(0.0F, 100.0F);
boost::variate_generator<boost::mt19937, boost::uniform_real<float> > gen2(rnd, ur_a2);
struct ibeta_inv_data_generator
{
boost::math::tuple<boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR> operator()(boost::math::ntl::RR ap, boost::math::ntl::RR bp, boost::math::ntl::RR x_)
boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t> operator()(mp_t ap, mp_t bp, mp_t x_)
{
float a = truncate_to_float(real_cast<float>(gen() * pow(boost::math::ntl::RR(10), ap)));
float b = truncate_to_float(real_cast<float>(gen() * pow(boost::math::ntl::RR(10), bp)));
float a = truncate_to_float(real_cast<float>(gen() * pow(mp_t(10), ap)));
float b = truncate_to_float(real_cast<float>(gen() * pow(mp_t(10), bp)));
float x = truncate_to_float(real_cast<float>(x_));
std::cout << a << " " << b << " " << x << std::flush;
boost::math::ntl::RR inv = boost::math::ibeta_inv(boost::math::ntl::RR(a), boost::math::ntl::RR(b), boost::math::ntl::RR(x));
mp_t inv = boost::math::ibeta_inv(mp_t(a), mp_t(b), mp_t(x));
std::cout << " " << inv << std::flush;
boost::math::ntl::RR invc = boost::math::ibetac_inv(boost::math::ntl::RR(a), boost::math::ntl::RR(b), boost::math::ntl::RR(x));
mp_t invc = boost::math::ibetac_inv(mp_t(a), mp_t(b), mp_t(x));
std::cout << " " << invc << std::endl;
return boost::math::make_tuple(a, b, x, inv, invc);
}
};
int cpp_main(int argc, char*argv [])
int main(int argc, char*argv [])
{
boost::math::ntl::RR::SetPrecision(1000);
boost::math::ntl::RR::SetOutputPrecision(100);
bool cont;
std::string line;
parameter_info<boost::math::ntl::RR> arg1, arg2, arg3;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1, arg2, arg3;
test_data<mp_t> data;
std::cout << "Welcome.\n"
"This program will generate spot tests for the inverse incomplete beta function:\n"
" ibeta_inv(a, p) and ibetac_inv(a, q)\n\n";
arg1 = make_periodic_param(boost::math::ntl::RR(-5), boost::math::ntl::RR(6), 11);
arg2 = make_periodic_param(boost::math::ntl::RR(-5), boost::math::ntl::RR(6), 11);
arg3 = make_random_param(boost::math::ntl::RR(0.0001), boost::math::ntl::RR(1), 10);
arg1 = make_periodic_param(mp_t(-5), mp_t(6), 11);
arg2 = make_periodic_param(mp_t(-5), mp_t(6), 11);
arg3 = make_random_param(mp_t(0.0001), mp_t(1), 10);
arg1.type |= dummy_param;
arg2.type |= dummy_param;
@@ -84,9 +78,9 @@ int cpp_main(int argc, char*argv [])
line = "ibeta_inv_data.ipp";
std::ofstream ofs(line.c_str());
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, "ibeta_inv_data");
return 0;
}
#endif

48
tools/ibeta_invab_data.cpp
View File

@@ -3,15 +3,13 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include <boost/math/special_functions/beta.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/tools/test.hpp>
#include <boost/lexical_cast.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include <boost/random.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
@@ -28,57 +26,54 @@ float force_truncate(const float* f)
return external_f;
}
float truncate_to_float(boost::math::ntl::RR r)
float truncate_to_float(mp_t r)
{
float f = boost::math::tools::real_cast<float>(r);
return force_truncate(&f);
}
std::tr1::mt19937 rnd;
std::tr1::uniform_real<float> ur_a(1.0F, 5.0F);
std::tr1::variate_generator<std::tr1::mt19937, std::tr1::uniform_real<float> > gen(rnd, ur_a);
std::tr1::uniform_real<float> ur_a2(0.0F, 100.0F);
std::tr1::variate_generator<std::tr1::mt19937, std::tr1::uniform_real<float> > gen2(rnd, ur_a2);
boost::mt19937 rnd;
boost::uniform_real<float> ur_a(1.0F, 5.0F);
boost::variate_generator<boost::mt19937, boost::uniform_real<float> > gen(rnd, ur_a);
boost::uniform_real<float> ur_a2(0.0F, 100.0F);
boost::variate_generator<boost::mt19937, boost::uniform_real<float> > gen2(rnd, ur_a2);
struct ibeta_inv_data_generator
{
boost::math::tuple<boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR> operator()
(boost::math::ntl::RR bp, boost::math::ntl::RR x_, boost::math::ntl::RR p_)
boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t, mp_t, mp_t> operator()
(mp_t bp, mp_t x_, mp_t p_)
{
float b = truncate_to_float(real_cast<float>(gen() * pow(boost::math::ntl::RR(10), bp)));
float b = truncate_to_float(real_cast<float>(gen() * pow(mp_t(10), bp)));
float x = truncate_to_float(real_cast<float>(x_));
float p = truncate_to_float(real_cast<float>(p_));
std::cout << b << " " << x << " " << p << std::flush;
boost::math::ntl::RR inv = boost::math::ibeta_inva(boost::math::ntl::RR(b), boost::math::ntl::RR(x), boost::math::ntl::RR(p));
mp_t inv = boost::math::ibeta_inva(mp_t(b), mp_t(x), mp_t(p));
std::cout << " " << inv << std::flush;
boost::math::ntl::RR invc = boost::math::ibetac_inva(boost::math::ntl::RR(b), boost::math::ntl::RR(x), boost::math::ntl::RR(p));
mp_t invc = boost::math::ibetac_inva(mp_t(b), mp_t(x), mp_t(p));
std::cout << " " << invc << std::endl;
boost::math::ntl::RR invb = boost::math::ibeta_invb(boost::math::ntl::RR(b), boost::math::ntl::RR(x), boost::math::ntl::RR(p));
mp_t invb = boost::math::ibeta_invb(mp_t(b), mp_t(x), mp_t(p));
std::cout << " " << invb << std::flush;
boost::math::ntl::RR invbc = boost::math::ibetac_invb(boost::math::ntl::RR(b), boost::math::ntl::RR(x), boost::math::ntl::RR(p));
mp_t invbc = boost::math::ibetac_invb(mp_t(b), mp_t(x), mp_t(p));
std::cout << " " << invbc << std::endl;
return boost::math::make_tuple(b, x, p, inv, invc, invb, invbc);
}
};
int cpp_main(int argc, char*argv [])
int main(int argc, char*argv [])
{
boost::math::ntl::RR::SetPrecision(1000);
boost::math::ntl::RR::SetOutputPrecision(100);
bool cont;
std::string line;
parameter_info<boost::math::ntl::RR> arg1, arg2, arg3;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1, arg2, arg3;
test_data<mp_t> data;
std::cout << "Welcome.\n"
"This program will generate spot tests for the inverse incomplete beta function:\n"
" ibeta_inva(a, p) and ibetac_inva(a, q)\n\n";
arg1 = make_periodic_param(boost::math::ntl::RR(-5), boost::math::ntl::RR(6), 11);
arg2 = make_random_param(boost::math::ntl::RR(0.0001), boost::math::ntl::RR(1), 10);
arg3 = make_random_param(boost::math::ntl::RR(0.0001), boost::math::ntl::RR(1), 10);
arg1 = make_periodic_param(mp_t(-5), mp_t(6), 11);
arg2 = make_random_param(mp_t(0.0001), mp_t(1), 10);
arg3 = make_random_param(mp_t(0.0001), mp_t(1), 10);
arg1.type |= dummy_param;
arg2.type |= dummy_param;
@@ -88,6 +83,7 @@ int cpp_main(int argc, char*argv [])
line = "ibeta_inva_data.ipp";
std::ofstream ofs(line.c_str());
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, "ibeta_inva_data");
return 0;

49
tools/igamma_data.cpp
View File

@@ -3,15 +3,12 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/tools/test.hpp>
#include <boost/lexical_cast.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
@@ -28,7 +25,7 @@ float force_truncate(const float* f)
return external_f;
}
float truncate_to_float(boost::math::ntl::RR r)
float truncate_to_float(mp_t r)
{
float f = boost::math::tools::real_cast<float>(r);
return force_truncate(&f);
@@ -42,14 +39,14 @@ float truncate_to_float(boost::math::ntl::RR r)
//
struct igamma_data_generator
{
boost::math::tuple<boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR, boost::math::ntl::RR> operator()(boost::math::ntl::RR a, boost::math::ntl::RR x)
boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t> operator()(mp_t a, mp_t x)
{
// very naively calculate spots:
boost::math::ntl::RR z;
mp_t z;
switch((int)real_cast<float>(x))
{
case 1:
z = truncate_to_float((std::min)(boost::math::ntl::RR(1), a/100));
z = truncate_to_float((std::min)(mp_t(1), a/100));
break;
case 2:
z = truncate_to_float(a / 2);
@@ -67,17 +64,17 @@ struct igamma_data_generator
z = truncate_to_float(a * 2);
break;
case 7:
z = truncate_to_float((std::max)(boost::math::ntl::RR(100), a*100));
z = truncate_to_float((std::max)(mp_t(100), a*100));
break;
default:
BOOST_ASSERT(0 == "Can't get here!!");
}
//boost::math::ntl::RR g = boost::math::tgamma(a);
boost::math::ntl::RR lg = boost::math::tgamma_lower(a, z);
boost::math::ntl::RR ug = boost::math::tgamma(a, z);
boost::math::ntl::RR rlg = boost::math::gamma_p(a, z);
boost::math::ntl::RR rug = boost::math::gamma_q(a, z);
//mp_t g = boost::math::tgamma(a);
mp_t lg = boost::math::tgamma_lower(a, z);
mp_t ug = boost::math::tgamma(a, z);
mp_t rlg = boost::math::gamma_p(a, z);
mp_t rug = boost::math::gamma_q(a, z);
return boost::math::make_tuple(z, ug, rug, lg, rlg);
}
@@ -85,10 +82,10 @@ struct igamma_data_generator
struct gamma_inverse_generator
{
boost::math::tuple<boost::math::ntl::RR, boost::math::ntl::RR> operator()(const boost::math::ntl::RR a, const boost::math::ntl::RR p)
boost::math::tuple<mp_t, mp_t> operator()(const mp_t a, const mp_t p)
{
boost::math::ntl::RR x1 = boost::math::gamma_p_inv(a, p);
boost::math::ntl::RR x2 = boost::math::gamma_q_inv(a, p);
mp_t x1 = boost::math::gamma_p_inv(a, p);
mp_t x2 = boost::math::gamma_q_inv(a, p);
std::cout << "Inverse for " << a << " " << p << std::endl;
return boost::math::make_tuple(x1, x2);
}
@@ -96,26 +93,23 @@ struct gamma_inverse_generator
struct gamma_inverse_generator_a
{
boost::math::tuple<boost::math::ntl::RR, boost::math::ntl::RR> operator()(const boost::math::ntl::RR x, const boost::math::ntl::RR p)
boost::math::tuple<mp_t, mp_t> operator()(const mp_t x, const mp_t p)
{
boost::math::ntl::RR x1 = boost::math::gamma_p_inva(x, p);
boost::math::ntl::RR x2 = boost::math::gamma_q_inva(x, p);
mp_t x1 = boost::math::gamma_p_inva(x, p);
mp_t x2 = boost::math::gamma_q_inva(x, p);
std::cout << "Inverse for " << x << " " << p << std::endl;
return boost::math::make_tuple(x1, x2);
}
};
int cpp_main(int argc, char*argv [])
int main(int argc, char*argv [])
{
boost::math::ntl::RR::SetPrecision(1000);
boost::math::ntl::RR::SetOutputPrecision(100);
bool cont;
std::string line;
parameter_info<boost::math::ntl::RR> arg1, arg2;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1, arg2;
test_data<mp_t> data;
if((argc >= 2) && (std::strcmp(argv[1], "-inverse") == 0))
{
@@ -151,7 +145,7 @@ int cpp_main(int argc, char*argv [])
}
else
{
arg2 = make_periodic_param(boost::math::ntl::RR(1), boost::math::ntl::RR(8), 7);
arg2 = make_periodic_param(mp_t(1), mp_t(8), 7);
arg2.type |= boost::math::tools::dummy_param;
std::cout << "Welcome.\n"
@@ -175,6 +169,7 @@ int cpp_main(int argc, char*argv [])
if(line == "")
line = "igamma_data.ipp";
std::ofstream ofs(line.c_str());
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, "igamma_data");
return 0;

38
tools/igamma_temme_large_coef.cpp
View File

@@ -3,16 +3,15 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/math/special_functions/log1p.hpp>
#include <boost/math/special_functions/erf.hpp>
#include <boost/math/constants/constants.hpp>
#include <map>
#include <iostream>
#include <iomanip>
#include "mp_t.hpp"
using namespace std;
using namespace NTL;
using namespace boost::math;
//
@@ -30,21 +29,21 @@ using namespace boost::math;
//
// Alpha:
//
RR alpha(unsigned k)
mp_t alpha(unsigned k)
{
static map<unsigned, RR> data;
static map<unsigned, mp_t> data;
if(data.empty())
{
data[1] = 1;
}
map<unsigned, RR>::const_iterator pos = data.find(k);
map<unsigned, mp_t>::const_iterator pos = data.find(k);
if(pos != data.end())
return (*pos).second;
//
// OK try and calculate the value:
//
RR result = alpha(k-1);
mp_t result = alpha(k-1);
for(unsigned j = 2; j <= k-1; ++j)
{
result -= j * alpha(j) * alpha(k-j+1);
@@ -54,15 +53,15 @@ RR alpha(unsigned k)
return result;
}
boost::math::ntl::RR gamma(unsigned k)
mp_t gamma(unsigned k)
{
static map<unsigned, boost::math::ntl::RR> data;
static map<unsigned, mp_t> data;
map<unsigned, boost::math::ntl::RR>::const_iterator pos = data.find(k);
map<unsigned, mp_t>::const_iterator pos = data.find(k);
if(pos != data.end())
return (*pos).second;
boost::math::ntl::RR result = (k&1) ? -1 : 1;
mp_t result = (k&1) ? -1 : 1;
for(unsigned i = 1; i <= (2 * k + 1); i += 2)
result *= i;
@@ -71,16 +70,16 @@ boost::math::ntl::RR gamma(unsigned k)
return result;
}
boost::math::ntl::RR Coeff(unsigned n, unsigned k)
mp_t Coeff(unsigned n, unsigned k)
{
map<unsigned, map<unsigned, boost::math::ntl::RR> > data;
map<unsigned, map<unsigned, mp_t> > data;
if(data.empty())
data[0][0] = boost::math::ntl::RR(-1) / 3;
data[0][0] = mp_t(-1) / 3;
map<unsigned, map<unsigned, boost::math::ntl::RR> >::const_iterator p1 = data.find(n);
map<unsigned, map<unsigned, mp_t> >::const_iterator p1 = data.find(n);
if(p1 != data.end())
{
map<unsigned, boost::math::ntl::RR>::const_iterator p2 = p1->second.find(k);
map<unsigned, mp_t>::const_iterator p2 = p1->second.find(k);
if(p2 != p1->second.end())
{
return p2->second;
@@ -93,12 +92,12 @@ boost::math::ntl::RR Coeff(unsigned n, unsigned k)
if(k == 0)
{
// special case:
boost::math::ntl::RR result = (n+2) * alpha(n+2);
mp_t result = (n+2) * alpha(n+2);
data[n][k] = result;
return result;
}
// general case:
boost::math::ntl::RR result = gamma(k) * Coeff(n, 0) + (n+2) * Coeff(n+2, k-1);
mp_t result = gamma(k) * Coeff(n, 0) + (n+2) * Coeff(n+2, k-1);
data[n][k] = result;
return result;
}
@@ -152,7 +151,7 @@ void calculate_terms(double sigma, double a, unsigned bits)
cin >> code;
int prec = 2 + (static_cast<double>(bits) * 3010LL)/10000;
RR::SetOutputPrecision(prec);
std::cout << std::scientific << std::setprecision(40);
if(code)
{
@@ -181,9 +180,6 @@ void calculate_terms(double sigma, double a, unsigned bits)
int main()
{
RR::SetOutputPrecision(50);
RR::SetPrecision(1000);
bool cont;
do{
cont = false;

29
tools/inv_hyp_data.cpp
View File

@@ -3,23 +3,20 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/tools/test.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
using namespace std;
struct asinh_data_generator
{
boost::math::ntl::RR operator()(boost::math::ntl::RR z)
mp_t operator()(mp_t z)
{
std::cout << z << " ";
boost::math::ntl::RR result = log(z + sqrt(z * z + 1));
mp_t result = log(z + sqrt(z * z + 1));
std::cout << result << std::endl;
return result;
}
@@ -27,10 +24,10 @@ struct asinh_data_generator
struct acosh_data_generator
{
boost::math::ntl::RR operator()(boost::math::ntl::RR z)
mp_t operator()(mp_t z)
{
std::cout << z << " ";
boost::math::ntl::RR result = log(z + sqrt(z * z - 1));
mp_t result = log(z + sqrt(z * z - 1));
std::cout << result << std::endl;
return result;
}
@@ -38,22 +35,19 @@ struct acosh_data_generator
struct atanh_data_generator
{
boost::math::ntl::RR operator()(boost::math::ntl::RR z)
mp_t operator()(mp_t z)
{
std::cout << z << " ";
boost::math::ntl::RR result = log((z + 1) / (1 - z)) / 2;
mp_t result = log((z + 1) / (1 - z)) / 2;
std::cout << result << std::endl;
return result;
}
};
int cpp_main(int argc, char*argv [])
int main(int argc, char*argv [])
{
boost::math::ntl::RR::SetPrecision(500);
boost::math::ntl::RR::SetOutputPrecision(40);
parameter_info<boost::math::ntl::RR> arg1;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1;
test_data<mp_t> data;
std::ofstream ofs;
bool cont;
@@ -79,6 +73,7 @@ int cpp_main(int argc, char*argv [])
if(line == "")
line = "asinh_data.ipp";
ofs.open(line.c_str());
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, "asinh_data");
data.clear();
@@ -103,6 +98,7 @@ int cpp_main(int argc, char*argv [])
line = "acosh_data.ipp";
ofs.close();
ofs.open(line.c_str());
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, "acosh_data");
data.clear();
@@ -127,6 +123,7 @@ int cpp_main(int argc, char*argv [])
line = "atanh_data.ipp";
ofs.close();
ofs.open(line.c_str());
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, "atanh_data");
return 0;

21
tools/laguerre_data.cpp
View File

@@ -3,14 +3,12 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/math/tools/test_data.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include <boost/math/special_functions/laguerre.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include <boost/tr1/random.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
using namespace boost::math;
@@ -34,20 +32,17 @@ boost::math::tuple<T, T, T, T> laguerre3_data(T n, T m, T x)
return boost::math::make_tuple(n, m, x, r1);
}
int cpp_main(int argc, char*argv [])
int main(int argc, char*argv [])
{
using namespace boost::math::tools;
boost::math::ntl::RR::SetOutputPrecision(50);
boost::math::ntl::RR::SetPrecision(1000);
parameter_info<boost::math::ntl::RR> arg1, arg2, arg3;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1, arg2, arg3;
test_data<mp_t> data;
bool cont;
std::string line;
if(argc < 1)
if(argc < 2)
return 1;
if(strcmp(argv[1], "--laguerre2") == 0)
@@ -60,7 +55,7 @@ int cpp_main(int argc, char*argv [])
arg1.type |= dummy_param;
arg2.type |= dummy_param;
data.insert(&laguerre2_data<boost::math::ntl::RR>, arg1, arg2);
data.insert(&laguerre2_data<mp_t>, arg1, arg2);
std::cout << "Any more data [y/n]?";
std::getline(std::cin, line);
@@ -81,7 +76,7 @@ int cpp_main(int argc, char*argv [])
arg2.type |= dummy_param;
arg3.type |= dummy_param;
data.insert(&laguerre3_data<boost::math::ntl::RR>, arg1, arg2, arg3);
data.insert(&laguerre3_data<mp_t>, arg1, arg2, arg3);
std::cout << "Any more data [y/n]?";
std::getline(std::cin, line);
@@ -98,7 +93,7 @@ int cpp_main(int argc, char*argv [])
line = "laguerre.ipp";
std::ofstream ofs(line.c_str());
line.erase(line.find('.'));
ofs << std::scientific;
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, line.c_str());
return 0;

49
tools/lanczos_generator.cpp
View File

@@ -3,10 +3,8 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/tools/precision.hpp>
#include <boost/math/tools/test.hpp>
#include <boost/math/tools/polynomial.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
@@ -16,7 +14,7 @@
#include <boost/math/special_functions/log1p.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/array.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include "mp_t.hpp"
//
// this is a sort of recursive include, since this file
@@ -3635,14 +3633,14 @@ lanczos_spot_data sweet_spots[] = {
//
// A bunch of helper functions used for calculating the coefficients:
//
boost::math::ntl::RR factorial(unsigned n)
mp_t factorial(unsigned n)
{
static boost::array<boost::array<boost::math::ntl::RR, 1>, 201> result;
static boost::array<boost::array<mp_t, 1>, 201> result;
static bool init = false;
if(!init)
{
unsigned k = 1;
boost::math::ntl::RR fact = 1;
mp_t fact = 1;
do{
result[k-1][0] = fact;
fact *= k++;
@@ -3655,7 +3653,7 @@ boost::math::ntl::RR factorial(unsigned n)
return result[n][0];
unsigned i = (unsigned)result.size()-1;
boost::math::ntl::RR r = result[i][0];
mp_t r = result[i][0];
while(i < n)
r *= ++i;
return r;
@@ -4050,7 +4048,7 @@ T get_max_error(const lanczos_info<T>& dat, R)
continue;
T gam = gamr;
T exp = tests[i][1];
T err = boost::math::tools::relative_error(gam, exp);
T err = relative_error(gam, exp);
if(err > max_error)
max_error = err;
}
@@ -4063,6 +4061,7 @@ T get_max_error(const lanczos_info<T>& dat, R)
template <class T>
void print_code(const lanczos_info<T>& l, const char* name)
{
std::cout << std::scientific << std::setprecision(std::numeric_limits<T>::digits10 + 3);
using namespace std;
lanczos_info<T> l2(l);
T factor = exp(l.r);
@@ -4226,7 +4225,7 @@ void print_code(const lanczos_info<T>& l, const char* name)
//
// Print out the test values:
//
void print_test_values(const std::vector<std::vector<boost::math::ntl::RR> >& v, const char* name, int offset = 1)
void print_test_values(const std::vector<std::vector<mp_t> >& v, const char* name, int offset = 1)
{
std::cout << "#define SC_(x) static_cast<T>(BOOST_JOIN(x, L))\n";
std::cout <<
@@ -4240,9 +4239,9 @@ void print_test_values(const std::vector<std::vector<boost::math::ntl::RR> >& v,
//
// Get the error for a specific approximation, and print out it's code:
//
void calculate_lanczos_spot(int n, boost::math::ntl::RR r, const char* suffix = "")
void calculate_lanczos_spot(int n, mp_t r, const char* suffix = "")
{
lanczos_info<boost::math::ntl::RR> info = generate_lanczos(n, r);
lanczos_info<mp_t> info = generate_lanczos(n, r);
// note error is calculated at high precision:
info.err = get_max_error(info, r);
print_code(info, suffix);
@@ -4255,14 +4254,14 @@ void find_best_lanczos(const char* name, T eps, int max_scan = 100)
{
using namespace std;
lanczos_info<boost::math::ntl::RR> best;
lanczos_info<mp_t> best;
best.err = 100; // best case had better be better than this!
for(int i = 0; i < sizeof(sweet_spots)/sizeof(sweet_spots[0]); ++i)
{
if((sweet_spots[i].err < eps*10) && (sweet_spots[i].N < max_scan))
{
lanczos_info<boost::math::ntl::RR> info = generate_lanczos(sweet_spots[i].N, boost::math::ntl::RR(sweet_spots[i].g));
boost::math::ntl::RR err = get_max_error(info, eps);
lanczos_info<mp_t> info = generate_lanczos(sweet_spots[i].N, mp_t(sweet_spots[i].g));
mp_t err = get_max_error(info, eps);
if(err/eps < 1000)
{
std::cout << sweet_spots[i].N << " " << sweet_spots[i].g << " " << err/eps << std::endl;
@@ -4278,7 +4277,7 @@ void find_best_lanczos(const char* name, T eps, int max_scan = 100)
print_code(best, name);
}
int cpp_main(int argc, char*argv [])
int main(int argc, char*argv [])
{
bool test_double(false), test_long(false), test_float(false), test_quad(false), spots(false), test_data(false);
@@ -4291,7 +4290,7 @@ int cpp_main(int argc, char*argv [])
" -long-double test type long double for the best approximation\n"
" -spots print out the best cases found in previous runs\n"
" -data print out the test data\n" << std::flush;
return 0;
return 1;
}
for(int i = 1; i < argc; ++i)
@@ -4310,9 +4309,6 @@ int cpp_main(int argc, char*argv [])
test_data = true;
}
boost::math::ntl::RR::SetPrecision(1000);
boost::math::ntl::RR::SetOutputPrecision(40);
if(spots)
{
// these are optimal N and R from Pugh:
@@ -4327,7 +4323,6 @@ int cpp_main(int argc, char*argv [])
calculate_lanczos_spot(6, 1.428456135094165802001953125);
calculate_lanczos_spot(13, 6.024680040776729583740234375);
calculate_lanczos_spot(17, 12.2252227365970611572265625);
boost::math::ntl::RR::SetOutputPrecision(90);
calculate_lanczos_spot(24, 20.3209821879863739013671875);
}
if(test_float)
@@ -4351,15 +4346,13 @@ int cpp_main(int argc, char*argv [])
{
std::cout << "Test Data follows:\n\n";
boost::math::ntl::RR::SetPrecision(1000);
std::vector<std::vector<boost::math::ntl::RR> > const & tests = get_test_data<boost::math::ntl::RR>();
std::vector<std::vector<mp_t> > const & tests = get_test_data<mp_t>();
print_test_values(tests, "factorials");
print_test_values(get_test_data_near_1<boost::math::ntl::RR>(), "near_1", 0);
print_test_values(get_test_data_near_2<boost::math::ntl::RR>(), "near_2", 0);
print_test_values(get_test_data_near_x<boost::math::ntl::RR>(boost::math::ntl::RR(0)), "near_0", 0);
print_test_values(get_test_data_near_x<boost::math::ntl::RR>(boost::math::ntl::RR(-10)), "near_m10", 0);
print_test_values(get_test_data_near_x<boost::math::ntl::RR>(boost::math::ntl::RR(-55)), "near_m55", 0);
print_test_values(get_test_data_near_1<mp_t>(), "near_1", 0);
print_test_values(get_test_data_near_2<mp_t>(), "near_2", 0);
print_test_values(get_test_data_near_x<mp_t>(mp_t(0)), "near_0", 0);
print_test_values(get_test_data_near_x<mp_t>(mp_t(-10)), "near_m10", 0);
print_test_values(get_test_data_near_x<mp_t>(mp_t(-55)), "near_m55", 0);
}
return 0;
}

24
tools/legendre_data.cpp
View File

@@ -3,14 +3,12 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/math/tools/test_data.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include <boost/math/special_functions/legendre.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include <boost/tr1/random.hpp>
#include <boost/random.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
using namespace boost::math;
@@ -29,23 +27,20 @@ boost::math::tuple<T, T, T, T> legendre_p_data(T n, T x)
template<class T>
boost::math::tuple<T, T, T, T> assoc_legendre_p_data(T n, T x)
{
static tr1::mt19937 r;
static boost::mt19937 r;
int l = real_cast<int>(floor(n));
tr1::uniform_int<> ui((std::max)(-l, -40), (std::min)(l, 40));
boost::uniform_int<> ui((std::max)(-l, -40), (std::min)(l, 40));
int m = ui(r);
T r1 = legendre_p(l, m, x);
return boost::math::make_tuple(n, m, x, r1);
}
int cpp_main(int argc, char*argv [])
int main(int argc, char*argv [])
{
using namespace boost::math::tools;
boost::math::ntl::RR::SetOutputPrecision(50);
boost::math::ntl::RR::SetPrecision(1000);
parameter_info<boost::math::ntl::RR> arg1, arg2;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1, arg2;
test_data<mp_t> data;
bool cont;
std::string line;
@@ -63,7 +58,7 @@ int cpp_main(int argc, char*argv [])
arg1.type |= dummy_param;
arg2.type |= dummy_param;
data.insert(&legendre_p_data<boost::math::ntl::RR>, arg1, arg2);
data.insert(&legendre_p_data<mp_t>, arg1, arg2);
std::cout << "Any more data [y/n]?";
std::getline(std::cin, line);
@@ -81,7 +76,7 @@ int cpp_main(int argc, char*argv [])
arg1.type |= dummy_param;
arg2.type |= dummy_param;
data.insert(&assoc_legendre_p_data<boost::math::ntl::RR>, arg1, arg2);
data.insert(&assoc_legendre_p_data<mp_t>, arg1, arg2);
std::cout << "Any more data [y/n]?";
std::getline(std::cin, line);
@@ -97,6 +92,7 @@ int cpp_main(int argc, char*argv [])
if(line == "")
line = "legendre_p.ipp";
std::ofstream ofs(line.c_str());
ofs << std::scientific << std::setprecision(40);
line.erase(line.find('.'));
write_code(ofs, data, line.c_str());

14
tools/log1p_expm1_data.cpp
View File

@@ -3,19 +3,18 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include <boost/math/special_functions/log1p.hpp>
#include <boost/math/special_functions/expm1.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
using namespace std;
struct data_generator
{
boost::math::tuple<boost::math::ntl::RR, boost::math::ntl::RR> operator()(boost::math::ntl::RR z)
boost::math::tuple<mp_t, mp_t> operator()(mp_t z)
{
return boost::math::make_tuple(boost::math::log1p(z), boost::math::expm1(z));
}
@@ -23,11 +22,8 @@ struct data_generator
int main(int argc, char* argv[])
{
boost::math::ntl::RR::SetPrecision(1000);
boost::math::ntl::RR::SetOutputPrecision(40);
parameter_info<boost::math::ntl::RR> arg1;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1;
test_data<mp_t> data;
std::cout << "Welcome.\n"
"This program will generate spot tests for the log1p and expm1 functions:\n\n";
@@ -52,7 +48,7 @@ int main(int argc, char* argv[])
if(line == "")
line = "log1p_expm1_data.ipp";
std::ofstream ofs(line.c_str());
ofs << std::scientific;
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, "log1p_expm1_data");
return 0;

71
tools/mp_t.hpp Normal file
View File

@@ -0,0 +1,71 @@
// (C) Copyright John Maddock 2014.
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#ifndef BOOST_MATH_TOOLS_MP_T
#define BOOST_MATH_TOOLS_MP_T
#ifndef BOOST_MATH_PRECISION
#define BOOST_MATH_PRECISION 1000
#endif
#if defined(BOOST_MATH_USE_MPFR)
#include <boost/multiprecision/mpfr.hpp>
typedef boost::multiprecision::number<boost::multiprecision::mpfr_float_backend<BOOST_MATH_PRECISION *301L / 1000L> > mp_t;
#else
#include <boost/multiprecision/cpp_bin_float.hpp>
typedef boost::multiprecision::number<boost::multiprecision::cpp_bin_float<BOOST_MATH_PRECISION, boost::multiprecision::digit_base_2> > mp_t;
#endif
inline mp_t ConvPrec(mp_t arg, int digits)
{
int e1, e2;
mp_t t = frexp(arg, &e1);
t = frexp(floor(ldexp(t, digits + 1)), &e2);
return ldexp(t, e1);
}
inline mp_t relative_error(mp_t a, mp_t b)
{
mp_t min_val = boost::math::tools::min_value<mp_t>();
mp_t max_val = boost::math::tools::max_value<mp_t>();
if((a != 0) && (b != 0))
{
// mp_tODO: use isfinite:
if(fabs(b) >= max_val)
{
if(fabs(a) >= max_val)
return 0; // one infinity is as good as another!
}
// If the result is denormalised, treat all denorms as equivalent:
if((a < min_val) && (a > 0))
a = min_val;
else if((a > -min_val) && (a < 0))
a = -min_val;
if((b < min_val) && (b > 0))
b = min_val;
else if((b > -min_val) && (b < 0))
b = -min_val;
return (std::max)(fabs((a - b) / a), fabs((a - b) / b));
}
// Handle special case where one or both are zero:
if(min_val == 0)
return fabs(a - b);
if(fabs(a) < min_val)
a = min_val;
if(fabs(b) < min_val)
b = min_val;
return (std::max)(fabs((a - b) / a), fabs((a - b) / b));
}
#endif // BOOST_MATH_TOOLS_MP_T

55
tools/rational_tests.cpp
View File

@@ -3,24 +3,23 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/tr1/random.hpp>
#include <boost/random.hpp>
#include <boost/math/tools/rational.hpp>
#include <iostream>
#include <fstream>
#include "mp_t.hpp"
int main()
{
using namespace boost::math;
using namespace boost::math::tools;
boost::math::ntl::RR::SetPrecision(500);
boost::math::ntl::RR::SetOutputPrecision(40);
std::cout << std::scientific << std::setprecision(40);
std::tr1::mt19937 rnd;
std::tr1::variate_generator<
std::tr1::mt19937,
std::tr1::uniform_int<> > gen(rnd, std::tr1::uniform_int<>(1, 12));
boost::mt19937 rnd;
boost::variate_generator<
boost::mt19937,
boost::uniform_int<> > gen(rnd, boost::uniform_int<>(1, 12));
for(unsigned i = 1; i < 12; ++i)
{
@@ -51,10 +50,10 @@ int main()
}
std::cout << " };\n";
boost::math::ntl::RR r1 = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(0.125), i);
boost::math::ntl::RR r2 = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(0.25), i);
boost::math::ntl::RR r3 = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(0.75), i);
boost::math::ntl::RR r4 = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(1) - boost::math::ntl::RR(1) / 64, i);
mp_t r1 = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(0.125), i);
mp_t r2 = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(0.25), i);
mp_t r3 = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(0.75), i);
mp_t r4 = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(1) - mp_t(1) / 64, i);
std::cout <<
" BOOST_CHECK_CLOSE(\n"
@@ -119,10 +118,10 @@ int main()
" static_cast<T>(" << r4 << "L),\n"
" tolerance);\n\n";
r1 = boost::math::tools::evaluate_even_polynomial(&coef[0], boost::math::ntl::RR(0.125), i);
r2 = boost::math::tools::evaluate_even_polynomial(&coef[0], boost::math::ntl::RR(0.25), i);
r3 = boost::math::tools::evaluate_even_polynomial(&coef[0], boost::math::ntl::RR(0.75), i);
r4 = boost::math::tools::evaluate_even_polynomial(&coef[0], boost::math::ntl::RR(1) - boost::math::ntl::RR(1) / 64, i);
r1 = boost::math::tools::evaluate_even_polynomial(&coef[0], mp_t(0.125), i);
r2 = boost::math::tools::evaluate_even_polynomial(&coef[0], mp_t(0.25), i);
r3 = boost::math::tools::evaluate_even_polynomial(&coef[0], mp_t(0.75), i);
r4 = boost::math::tools::evaluate_even_polynomial(&coef[0], mp_t(1) - mp_t(1) / 64, i);
std::cout <<
" BOOST_CHECK_CLOSE(\n"
@@ -189,10 +188,10 @@ int main()
if(i > 1)
{
r1 = boost::math::tools::evaluate_odd_polynomial(&coef[0], boost::math::ntl::RR(0.125), i);
r2 = boost::math::tools::evaluate_odd_polynomial(&coef[0], boost::math::ntl::RR(0.25), i);
r3 = boost::math::tools::evaluate_odd_polynomial(&coef[0], boost::math::ntl::RR(0.75), i);
r4 = boost::math::tools::evaluate_odd_polynomial(&coef[0], boost::math::ntl::RR(1) - boost::math::ntl::RR(1) / 64, i);
r1 = boost::math::tools::evaluate_odd_polynomial(&coef[0], mp_t(0.125), i);
r2 = boost::math::tools::evaluate_odd_polynomial(&coef[0], mp_t(0.25), i);
r3 = boost::math::tools::evaluate_odd_polynomial(&coef[0], mp_t(0.75), i);
r4 = boost::math::tools::evaluate_odd_polynomial(&coef[0], mp_t(1) - mp_t(1) / 64, i);
std::cout <<
" BOOST_CHECK_CLOSE(\n"
@@ -258,10 +257,10 @@ int main()
" tolerance);\n\n";
}
r1 = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(0.125), i);
r2 = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(0.25), i);
r3 = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(0.75), i);
r4 = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(1) - boost::math::ntl::RR(1) / 64, i);
r1 = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(0.125), i);
r2 = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(0.25), i);
r3 = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(0.75), i);
r4 = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(1) - mp_t(1) / 64, i);
coef.clear();
for(unsigned j = 0; j < i; ++j)
@@ -290,10 +289,10 @@ int main()
}
std::cout << " };\n";
boost::math::ntl::RR r1d = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(0.125), i);
boost::math::ntl::RR r2d = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(0.25), i);
boost::math::ntl::RR r3d = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(0.75), i);
boost::math::ntl::RR r4d = boost::math::tools::evaluate_polynomial(&coef[0], boost::math::ntl::RR(1) - boost::math::ntl::RR(1) / 64, i);
mp_t r1d = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(0.125), i);
mp_t r2d = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(0.25), i);
mp_t r3d = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(0.75), i);
mp_t r4d = boost::math::tools::evaluate_polynomial(&coef[0], mp_t(1) - mp_t(1) / 64, i);
std::cout <<
" BOOST_CHECK_CLOSE(\n"

25
tools/spherical_harmonic_data.cpp
View File

@@ -3,13 +3,11 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/math/tools/test_data.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include <boost/math/special_functions/spherical_harmonic.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include <boost/tr1/random.hpp>
#include <boost/random.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
using namespace boost::math;
@@ -28,13 +26,13 @@ float truncate_to_float(float const * pf)
template<class T>
boost::math::tuple<T, T, T, T, T, T> spherical_harmonic_data(T i)
{
static tr1::mt19937 r;
static boost::mt19937 r;
int n = real_cast<int>(floor(i));
std::tr1::uniform_int<> ui(0, (std::min)(n, 40));
boost::uniform_int<> ui(0, (std::min)(n, 40));
int m = ui(r);
std::tr1::uniform_real<float> ur(-2*constants::pi<float>(), 2*constants::pi<float>());
boost::uniform_real<float> ur(-2*constants::pi<float>(), 2*constants::pi<float>());
float _theta = ur(r);
float _phi = ur(r);
T theta = truncate_to_float(&_theta);
@@ -45,15 +43,12 @@ boost::math::tuple<T, T, T, T, T, T> spherical_harmonic_data(T i)
return boost::math::make_tuple(n, m, theta, phi, r1, r2);
}
int cpp_main(int argc, char*argv [])
int main(int argc, char*argv [])
{
using namespace boost::math::tools;
boost::math::ntl::RR::SetOutputPrecision(50);
boost::math::ntl::RR::SetPrecision(1000);
parameter_info<boost::math::ntl::RR> arg1, arg2, arg3;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1, arg2, arg3;
test_data<mp_t> data;
bool cont;
std::string line;
@@ -68,7 +63,7 @@ int cpp_main(int argc, char*argv [])
arg2.type |= dummy_param;
arg3 = arg2;
data.insert(&spherical_harmonic_data<boost::math::ntl::RR>, arg1);
data.insert(&spherical_harmonic_data<mp_t>, arg1);
std::cout << "Any more data [y/n]?";
std::getline(std::cin, line);
@@ -83,7 +78,7 @@ int cpp_main(int argc, char*argv [])
line = "spherical_harmonic.ipp";
std::ofstream ofs(line.c_str());
line.erase(line.find('.'));
ofs << std::scientific;
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, line.c_str());
return 0;

28
tools/tgamma_ratio_data.cpp
View File

@@ -3,44 +3,38 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/tools/test.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
using namespace std;
boost::math::tuple<boost::math::ntl::RR, boost::math::ntl::RR>
tgamma_ratio(const boost::math::ntl::RR& a, const boost::math::ntl::RR& delta)
boost::math::tuple<mp_t, mp_t>
tgamma_ratio(const mp_t& a, const mp_t& delta)
{
if(delta > a)
throw std::domain_error("");
boost::math::ntl::RR tg = boost::math::tgamma(a);
boost::math::ntl::RR r1 = tg / boost::math::tgamma(a + delta);
boost::math::ntl::RR r2 = tg / boost::math::tgamma(a - delta);
mp_t tg = boost::math::tgamma(a);
mp_t r1 = tg / boost::math::tgamma(a + delta);
mp_t r2 = tg / boost::math::tgamma(a - delta);
if((r1 > (std::numeric_limits<float>::max)()) || (r2 > (std::numeric_limits<float>::max)()))
throw std::domain_error("");
return boost::math::make_tuple(r1, r2);
}
boost::math::ntl::RR tgamma_ratio2(const boost::math::ntl::RR& a, const boost::math::ntl::RR& b)
mp_t tgamma_ratio2(const mp_t& a, const mp_t& b)
{
return boost::math::tgamma(a) / boost::math::tgamma(b);
}
int cpp_main(int argc, char*argv [])
int main(int argc, char*argv [])
{
boost::math::ntl::RR::SetPrecision(1000);
boost::math::ntl::RR::SetOutputPrecision(40);
parameter_info<boost::math::ntl::RR> arg1, arg2;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1, arg2;
test_data<mp_t> data;
bool cont;
std::string line;
@@ -88,7 +82,7 @@ int cpp_main(int argc, char*argv [])
if(line == "")
line = "tgamma_ratio_data.ipp";
std::ofstream ofs(line.c_str());
ofs << std::scientific;
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, "tgamma_ratio_data");
return 0;

23
tools/zeta_data.cpp
View File

@@ -3,24 +3,21 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/bindings/rr.hpp>
#include <boost/test/included/prg_exec_monitor.hpp>
#include <boost/math/special_functions/zeta.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/tools/test.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include "mp_t.hpp"
using namespace boost::math::tools;
using namespace std;
struct zeta_data_generator
{
boost::math::ntl::RR operator()(boost::math::ntl::RR z)
mp_t operator()(mp_t z)
{
std::cout << z << " ";
boost::math::ntl::RR result = boost::math::zeta(z);
mp_t result = boost::math::zeta(z);
std::cout << result << std::endl;
return result;
}
@@ -28,23 +25,20 @@ struct zeta_data_generator
struct zeta_data_generator2
{
boost::math::tuple<boost::math::ntl::RR, boost::math::ntl::RR> operator()(boost::math::ntl::RR z)
boost::math::tuple<mp_t, mp_t> operator()(mp_t z)
{
std::cout << -z << " ";
boost::math::ntl::RR result = boost::math::zeta(-z);
mp_t result = boost::math::zeta(-z);
std::cout << result << std::endl;
return boost::math::make_tuple(-z, result);
}
};
int cpp_main(int argc, char*argv [])
int main(int argc, char*argv [])
{
boost::math::ntl::RR::SetPrecision(500);
boost::math::ntl::RR::SetOutputPrecision(40);
parameter_info<boost::math::ntl::RR> arg1;
test_data<boost::math::ntl::RR> data;
parameter_info<mp_t> arg1;
test_data<mp_t> data;
bool cont;
std::string line;
@@ -70,6 +64,7 @@ int cpp_main(int argc, char*argv [])
if(line == "")
line = "zeta_data.ipp";
std::ofstream ofs(line.c_str());
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, "zeta_data");
return 0;