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

@@ -578,28 +578,28 @@ T rising_factorial_ratio(T a, T b, int k)
 template <class T>
 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<T>::value >= 100);
 };
 template <>
 struct Pn_size<float>
 {
    BOOST_STATIC_CONSTANT(unsigned, value = 15); // ~8-15 digit accuracy
-   BOOST_STATIC_ASSERT(::boost::math::max_factorial<float>::value >= 15);
+   BOOST_STATIC_ASSERT(::boost::math::max_factorial<float>::value >= 30);
 };
 template <>
 struct Pn_size<double>
 {
    BOOST_STATIC_CONSTANT(unsigned, value = 30); // 16-20 digit accuracy
-   BOOST_STATIC_ASSERT(::boost::math::max_factorial<double>::value >= 30);
+   BOOST_STATIC_ASSERT(::boost::math::max_factorial<double>::value >= 60);
 };
 template <>
 struct Pn_size<long double>
 {
    BOOST_STATIC_CONSTANT(unsigned, value = 50); // ~35-50 digit accuracy
-   BOOST_STATIC_ASSERT(::boost::math::max_factorial<long double>::value >= 50);
+   BOOST_STATIC_ASSERT(::boost::math::max_factorial<long double>::value >= 100);
 };
 
 template <class T, class L>

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 <boost/math/tools/roots.hpp>
+
+namespace NTL{
+
+   //
+   // Inverse trig functions:
+   //
+   struct asin_root
+   {
+      asin_root(NTL::RR const& target) : t(target){}
+
+      std::tr1::tuple<NTL::RR, NTL::RR> 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<NTL::RR>()/2),
+         NTL::RR(boost::math::constants::pi<NTL::RR>()/2),
+         boost::math::tools::digits<NTL::RR>());
+   }
+
+} // namespace NTL
+
+
 #endif // BOOST_MATH_TOOLS_NTL_HPP
 
 

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 <boost/math/tools/ntl.hpp>
+#include <boost/test/included/test_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 "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<float>(r);
+   return force_truncate(&f);
+}
+
+struct gamma_inverse_generator_a
+{
+   std::tr1::tuple<NTL::RR, NTL::RR> 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<NTL::RR> arg1, arg2;
+   test_data<NTL::RR> 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>(NTL::RR(0), -4, 24);
+   arg2 = make_random_param<NTL::RR>(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;
+}
+

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 <boost/math/tools/ntl.hpp>
+#include <boost/test/included/test_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 "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<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);
+
+struct ibeta_inv_data_generator
+{
+   std::tr1::tuple<NTL::RR, NTL::RR, NTL::RR, NTL::RR, NTL::RR> operator()(NTL::RR ap, NTL::RR bp, NTL::RR x_)
+   {
+      float a = truncate_to_float(real_cast<float>(gen() * pow(NTL::RR(10), ap)));      
+      float b = truncate_to_float(real_cast<float>(gen() * pow(NTL::RR(10), bp))); 
+      float x = truncate_to_float(real_cast<float>(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<NTL::RR> arg1, arg2, arg3;
+   test_data<NTL::RR> 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