text: //  (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/bindings/rr.hpp>
  #include <boost/test/included/test_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>
  
  using namespace boost::math::tools;
  using namespace boost::math;
  using namespace std;
  
  template <class T>
  struct ibeta_fraction1_t
  {
     typedef std::pair<T, T> result_type;
  
     ibeta_fraction1_t(T a_, T b_, T x_) : a(a_), b(b_), x(x_), k(1) {}
  
     result_type operator()()
     {
        T aN;
        if(k & 1)
        {
           int m = (k - 1) / 2;
           aN = -(a + m) * (a + b + m) * x;
           aN /= a + 2*m;
           aN /= a + 2*m + 1;
        }
        else
        {
           int m = k / 2;
           aN = m * (b - m) *x;
           aN /= a + 2*m - 1;
           aN /= a + 2*m;
        }
        ++k;
        return std::make_pair(aN, T(1));
     }
  
  private:
     T a, b, x;
     int k;
  };
  
  //
  // This function caches previous calls to beta
  // just so we can speed things up a bit:
  //
  template <class T>
  T get_beta(T a, T b)
  {
     static std::map<std::pair<T, T>, T> m;
  
     if(a < b)
        std::swap(a, b);
  
     std::pair<T, T> p(a, b);
     typename std::map<std::pair<T, T>, T>::const_iterator i = m.find(p);
     if(i != m.end())
        return i->second;
  
     T r = beta(a, b);
     p.first = a;
     p.second = b;
     m[p] = r;
  
     return r;
  }
  
  //
  // compute the continued fraction:
  //
  template <class T>
  T get_ibeta_fraction1(T a, T b, T x)
  {
     ibeta_fraction1_t<T> f(a, b, x);
     T fract = boost::math::tools::continued_fraction_a(f, boost::math::policies::digits<T, boost::math::policies::policy<> >());
     T denom = (a * (fract + 1));
     T num = pow(x, a) * pow(1 - x, b);
     if(num == 0)
        return 0;
     else if(denom == 0)
        return -1;
     return num / denom;
  }
  //
  // calculate the incomplete beta from the fraction:
  //
  template <class T>
  std::pair<T,T> ibeta_fraction1(T a, T b, T x)
  {
     T bet = get_beta(a, b);
     if(x > ((a+1)/(a+b+2)))
     {
        T fract = get_ibeta_fraction1(b, a, 1-x);
        if(fract/bet > 0.75)
        {
           fract = get_ibeta_fraction1(a, b, x);
           return std::make_pair(fract, bet - fract);
        }
        return std::make_pair(bet - fract, fract);
     }
     T fract = get_ibeta_fraction1(a, b, x);
     if(fract/bet > 0.75)
     {
        fract = get_ibeta_fraction1(b, a, 1-x);
        return std::make_pair(bet - fract, fract);
     }
     return std::make_pair(fract, bet - fract);
  
  }
  //
  // calculate the regularised incomplete beta from the fraction:
  //
  template <class T>
  std::pair<T,T> ibeta_fraction1_regular(T a, T b, T x)
  {
     T bet = get_beta(a, b);
     if(x > ((a+1)/(a+b+2)))
     {
        T fract = get_ibeta_fraction1(b, a, 1-x);
        if(fract == 0)
           bet = 1;  // normalise so we don't get 0/0
        else if(bet == 0)
           return std::make_pair(T(-1), T(-1)); // Yikes!!
        if(fract / bet > 0.75)
        {
           fract = get_ibeta_fraction1(a, b, x);
           return std::make_pair(fract / bet, 1 - (fract / bet));
        }
        return std::make_pair(1 - (fract / bet), fract / bet);
     }
     T fract = get_ibeta_fraction1(a, b, x);
     if(fract / bet > 0.75)
     {
        fract = get_ibeta_fraction1(b, a, 1-x);
        return std::make_pair(1 - (fract / bet), fract / bet);
     }
     return std::make_pair(fract / bet, 1 - (fract / bet));
  }
  
  //
  // we absolutely must trunctate the input values to float
  // precision: we have to be certain that the input values
  // can be represented exactly in whatever width floating
  // point type we are testing, otherwise the output will 
  // necessarily be off.
  //
  float external_f;
  float force_truncate(const float* f)
  {
     external_f = *f;
     return external_f;
  }
  
  float truncate_to_float(boost::math::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 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_)
     {
        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 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));
        return boost::math::make_tuple(a, b, x, ib_full.first, ib_full.second, ib_reg.first, ib_reg.second);
     }
  };
  
  // 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::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);
        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));
        /*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));
        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_)
     {
        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));
        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_)
     {
        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));
        return boost::math::make_tuple(a, b, x, ib_full.first, ib_full.second, ib_reg.first, ib_reg.second);
     }
  };
  
  int test_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;
  
     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.type |= dummy_param;
     arg2.type |= dummy_param;
     arg3.type |= dummy_param;
     arg4.type |= dummy_param;
     arg5.type |= dummy_param;
  
     // comment out all but one of the following when running
     // or this program will take forever to complete!
     //data.insert(beta_data_generator(), arg1, arg2, arg3);
     //data.insert(beta_data_generator_medium(), arg4);
     //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;
     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]);
        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]);
        if(v2 > 1e-40)
        {
           std::cout << v2 << std::endl;
        }
        if(v4 > 1e-60)
        {
           std::cout << v4 << std::endl;
        }
        ++ i;
     }
  
     std::cout << "Enter name of test data file [default=ibeta_data.ipp]";
     std::string line;
     std::getline(std::cin, line);
     boost::algorithm::trim(line);
     if(line == "")
        line = "ibeta_data.ipp";
     std::ofstream ofs(line.c_str());
     write_code(ofs, data, "ibeta_data");
     
     return 0;
  }
  
  
  