From 4ff84d85ab188a4d2ddc30d51066897609e99f2c Mon Sep 17 00:00:00 2001 From: John Maddock Date: Mon, 14 Aug 2006 17:29:59 +0000 Subject: [PATCH] Added sample programs that generate test data for incomplete beta and gamma inverses. Added asin to ntl.hpp: needed for incomplete beta inverse. Fixed limits on table size in the incomplete beta. [SVN r3146] --- include/boost/math/special_functions/beta.hpp | 10 +- include/boost/math/tools/ntl.hpp | 40 ++++++++ tools/gamma_P_inva_data.cpp | 75 +++++++++++++++ tools/ibeta_inv_data.cpp | 93 +++++++++++++++++++ 4 files changed, 213 insertions(+), 5 deletions(-) create mode 100644 tools/gamma_P_inva_data.cpp create mode 100644 tools/ibeta_inv_data.cpp diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index c6a5acbee..1be4217fb 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -578,28 +578,28 @@ T rising_factorial_ratio(T a, T b, int k) template struct Pn_size { - // This is likely to be enough for ~100 digit accuracy + // This is likely to be enough for ~35-50 digit accuracy // but it's hard to quantify exactly: - BOOST_STATIC_CONSTANT(unsigned, value = 100); + BOOST_STATIC_CONSTANT(unsigned, value = 50); BOOST_STATIC_ASSERT(::boost::math::max_factorial::value >= 100); }; template <> struct Pn_size { BOOST_STATIC_CONSTANT(unsigned, value = 15); // ~8-15 digit accuracy - BOOST_STATIC_ASSERT(::boost::math::max_factorial::value >= 15); + BOOST_STATIC_ASSERT(::boost::math::max_factorial::value >= 30); }; template <> struct Pn_size { BOOST_STATIC_CONSTANT(unsigned, value = 30); // 16-20 digit accuracy - BOOST_STATIC_ASSERT(::boost::math::max_factorial::value >= 30); + BOOST_STATIC_ASSERT(::boost::math::max_factorial::value >= 60); }; template <> struct Pn_size { BOOST_STATIC_CONSTANT(unsigned, value = 50); // ~35-50 digit accuracy - BOOST_STATIC_ASSERT(::boost::math::max_factorial::value >= 50); + BOOST_STATIC_ASSERT(::boost::math::max_factorial::value >= 100); }; template diff --git a/include/boost/math/tools/ntl.hpp b/include/boost/math/tools/ntl.hpp index 6c935d301..c58d94712 100644 --- a/include/boost/math/tools/ntl.hpp +++ b/include/boost/math/tools/ntl.hpp @@ -55,6 +55,7 @@ namespace NTL r.e += exp; return r; } + } // namespace NTL namespace boost{ namespace math{ namespace tools{ @@ -267,6 +268,45 @@ void setprecision(std::ostream& /*os*/, NTL::quad_float, int p) } // namespace tools } // namespace math } // namespaceboost + +// +// The following *must* all occur after the definitions above +// since we rely on them being already defined: +// +#include + +namespace NTL{ + + // + // Inverse trig functions: + // + struct asin_root + { + asin_root(NTL::RR const& target) : t(target){} + + std::tr1::tuple operator()(NTL::RR const& p) + { + NTL::RR f0 = sin(p) - t; + NTL::RR f1 = cos(p); + return std::tr1::make_tuple(f0, f1); + } + private: + NTL::RR t; + }; + + NTL::RR asin(NTL::RR z) + { + return boost::math::tools::newton_raphson_iterate( + asin_root(z), + NTL::RR(0), + NTL::RR(-boost::math::constants::pi()/2), + NTL::RR(boost::math::constants::pi()/2), + boost::math::tools::digits()); + } + +} // namespace NTL + + #endif // BOOST_MATH_TOOLS_NTL_HPP diff --git a/tools/gamma_P_inva_data.cpp b/tools/gamma_P_inva_data.cpp new file mode 100644 index 000000000..a7a378f90 --- /dev/null +++ b/tools/gamma_P_inva_data.cpp @@ -0,0 +1,75 @@ +// (C) Copyright John Maddock 2006. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include +#include +#include +#include + +#include +#include "ntl_rr_lanczos.hpp" + +using namespace boost::math::tools; + +// +// Force trunctation to float precision of input values: +// we must ensure that the input values are exactly representable +// in whatever type we are testing, or the output values will all +// be thrown off: +// +float external_f; +float force_truncate(const float* f) +{ + external_f = *f; + return external_f; +} + +float truncate_to_float(NTL::RR r) +{ + float f = boost::math::tools::real_cast(r); + return force_truncate(&f); +} + +struct gamma_inverse_generator_a +{ + std::tr1::tuple operator()(const NTL::RR x, const NTL::RR p) + { + NTL::RR x1 = boost::math::gamma_P_inva(x, p); + NTL::RR x2 = boost::math::gamma_Q_inva(x, p); + std::cout << "Inverse for " << x << " " << p << std::endl; + return std::tr1::make_tuple(x1, x2); + } +}; + + +int test_main(int argc, char*argv []) +{ + NTL::RR::SetPrecision(1000); + NTL::RR::SetOutputPrecision(100); + + bool cont; + std::string line; + + parameter_info arg1, arg2; + test_data 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(NTL::RR(0), -4, 24); + arg2 = make_random_param(NTL::RR(0), NTL::RR(1), 15); + data.insert(gamma_inverse_generator_a(), arg1, arg2); + + line = "igamma_inva_data.ipp"; + std::ofstream ofs(line.c_str()); + write_code(ofs, data, "igamma_inva_data"); + + return 0; +} + diff --git a/tools/ibeta_inv_data.cpp b/tools/ibeta_inv_data.cpp new file mode 100644 index 000000000..a37945a8a --- /dev/null +++ b/tools/ibeta_inv_data.cpp @@ -0,0 +1,93 @@ +// (C) Copyright John Maddock 2006. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#if 1 +#include +#include +#include +#include +#include +#include +#include + +#include +#include "ntl_rr_lanczos.hpp" + +using namespace boost::math::tools; + +// +// Force trunctation to float precision of input values: +// we must ensure that the input values are exactly representable +// in whatever type we are testing, or the output values will all +// be thrown off: +// +float external_f; +float force_truncate(const float* f) +{ + external_f = *f; + return external_f; +} + +float truncate_to_float(NTL::RR r) +{ + float f = boost::math::tools::real_cast(r); + return force_truncate(&f); +} + +std::tr1::mt19937 rnd; +std::tr1::uniform_real ur_a(1.0F, 5.0F); +std::tr1::variate_generator > gen(rnd, ur_a); +std::tr1::uniform_real ur_a2(0.0F, 100.0F); +std::tr1::variate_generator > gen2(rnd, ur_a2); + +struct ibeta_inv_data_generator +{ + std::tr1::tuple operator()(NTL::RR ap, NTL::RR bp, NTL::RR x_) + { + float a = truncate_to_float(real_cast(gen() * pow(NTL::RR(10), ap))); + float b = truncate_to_float(real_cast(gen() * pow(NTL::RR(10), bp))); + float x = truncate_to_float(real_cast(x_)); + std::cout << a << " " << b << " " << x << std::flush; + NTL::RR inv = boost::math::ibeta_inv(NTL::RR(a), NTL::RR(b), NTL::RR(x)); + std::cout << " " << inv << std::flush; + NTL::RR invc = boost::math::ibetac_inv(NTL::RR(a), NTL::RR(b), NTL::RR(x)); + std::cout << " " << invc << std::endl; + return std::tr1::make_tuple(a, b, x, inv, invc); + } +}; + +int test_main(int argc, char*argv []) +{ + NTL::RR::SetPrecision(1000); + NTL::RR::SetOutputPrecision(100); + + bool cont; + std::string line; + + parameter_info arg1, arg2, arg3; + test_data 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(NTL::RR(-5), NTL::RR(6), 11); + arg2 = make_periodic_param(NTL::RR(-5), NTL::RR(6), 11); + arg3 = make_random_param(NTL::RR(0.0001), NTL::RR(1), 10); + + arg1.type |= dummy_param; + arg2.type |= dummy_param; + arg3.type |= dummy_param; + + data.insert(ibeta_inv_data_generator(), arg1, arg2, arg3); + + line = "ibeta_inv_data.ipp"; + std::ofstream ofs(line.c_str()); + write_code(ofs, data, "ibeta_inv_data"); + + return 0; +} + +#endif