ceph/src/boost/libs/math/tools/ibeta_data.cpp

   1 //  (C) Copyright John Maddock 2006.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #include "mp_t.hpp"
   7 #include <boost/math/special_functions/gamma.hpp>
   8 #include <boost/math/special_functions/beta.hpp>
   9 #include <boost/math/constants/constants.hpp>
  10 #include <boost/lexical_cast.hpp>
  11 #include <fstream>
  12 #include <map>
  13 #include <boost/math/tools/test_data.hpp>
  14 #include <boost/random.hpp>
  15
  16 using namespace boost::math::tools;
  17 using namespace boost::math;
  18 using namespace std;
  19
  20 template <class T>
  21 struct ibeta_fraction1_t
  22 {
  23    typedef std::pair<T, T> result_type;
  24
  25    ibeta_fraction1_t(T a_, T b_, T x_) : a(a_), b(b_), x(x_), k(1) {}
  26
  27    result_type operator()()
  28    {
  29       T aN;
  30       if(k & 1)
  31       {
  32          int m = (k - 1) / 2;
  33          aN = -(a + m) * (a + b + m) * x;
  34          aN /= a + 2*m;
  35          aN /= a + 2*m + 1;
  36       }
  37       else
  38       {
  39          int m = k / 2;
  40          aN = m * (b - m) *x;
  41          aN /= a + 2*m - 1;
  42          aN /= a + 2*m;
  43       }
  44       ++k;
  45       return std::make_pair(aN, T(1));
  46    }
  47
  48 private:
  49    T a, b, x;
  50    int k;
  51 };
  52
  53 //
  54 // This function caches previous calls to beta
  55 // just so we can speed things up a bit:
  56 //
  57 template <class T>
  58 T get_beta(T a, T b)
  59 {
  60    static std::map<std::pair<T, T>, T> m;
  61
  62    if(a < b)
  63       std::swap(a, b);
  64
  65    std::pair<T, T> p(a, b);
  66    typename std::map<std::pair<T, T>, T>::const_iterator i = m.find(p);
  67    if(i != m.end())
  68       return i->second;
  69
  70    T r = beta(a, b);
  71    p.first = a;
  72    p.second = b;
  73    m[p] = r;
  74
  75    return r;
  76 }
  77
  78 //
  79 // compute the continued fraction:
  80 //
  81 template <class T>
  82 T get_ibeta_fraction1(T a, T b, T x)
  83 {
  84    ibeta_fraction1_t<T> f(a, b, x);
  85    T fract = boost::math::tools::continued_fraction_a(f, boost::math::policies::digits<T, boost::math::policies::policy<> >());
  86    T denom = (a * (fract + 1));
  87    T num = pow(x, a) * pow(1 - x, b);
  88    if(num == 0)
  89       return 0;
  90    else if(denom == 0)
  91       return -1;
  92    return num / denom;
  93 }
  94 //
  95 // calculate the incomplete beta from the fraction:
  96 //
  97 template <class T>
  98 std::pair<T,T> ibeta_fraction1(T a, T b, T x)
  99 {
 100    T bet = get_beta(a, b);
 101    if(x > ((a+1)/(a+b+2)))
 102    {
 103       T fract = get_ibeta_fraction1(b, a, 1-x);
 104       if(fract/bet > 0.75)
 105       {
 106          fract = get_ibeta_fraction1(a, b, x);
 107          return std::make_pair(fract, bet - fract);
 108       }
 109       return std::make_pair(bet - fract, fract);
 110    }
 111    T fract = get_ibeta_fraction1(a, b, x);
 112    if(fract/bet > 0.75)
 113    {
 114       fract = get_ibeta_fraction1(b, a, 1-x);
 115       return std::make_pair(bet - fract, fract);
 116    }
 117    return std::make_pair(fract, bet - fract);
 118
 119 }
 120 //
 121 // calculate the regularised incomplete beta from the fraction:
 122 //
 123 template <class T>
 124 std::pair<T,T> ibeta_fraction1_regular(T a, T b, T x)
 125 {
 126    T bet = get_beta(a, b);
 127    if(x > ((a+1)/(a+b+2)))
 128    {
 129       T fract = get_ibeta_fraction1(b, a, 1-x);
 130       if(fract == 0)
 131          bet = 1;  // normalise so we don't get 0/0
 132       else if(bet == 0)
 133          return std::make_pair(T(-1), T(-1)); // Yikes!!
 134       if(fract / bet > 0.75)
 135       {
 136          fract = get_ibeta_fraction1(a, b, x);
 137          return std::make_pair(fract / bet, 1 - (fract / bet));
 138       }
 139       return std::make_pair(1 - (fract / bet), fract / bet);
 140    }
 141    T fract = get_ibeta_fraction1(a, b, x);
 142    if(fract / bet > 0.75)
 143    {
 144       fract = get_ibeta_fraction1(b, a, 1-x);
 145       return std::make_pair(1 - (fract / bet), fract / bet);
 146    }
 147    return std::make_pair(fract / bet, 1 - (fract / bet));
 148 }
 149
 150 //
 151 // we absolutely must truncate the input values to float
 152 // precision: we have to be certain that the input values
 153 // can be represented exactly in whatever width floating
 154 // point type we are testing, otherwise the output will
 155 // necessarily be off.
 156 //
 157 float external_f;
 158 float force_truncate(const float* f)
 159 {
 160    external_f = *f;
 161    return external_f;
 162 }
 163
 164 float truncate_to_float(mp_t r)
 165 {
 166    float f = boost::math::tools::real_cast<float>(r);
 167    return force_truncate(&f);
 168 }
 169
 170 boost::mt19937 rnd;
 171 boost::uniform_real<float> ur_a(1.0F, 5.0F);
 172 boost::variate_generator<boost::mt19937, boost::uniform_real<float> > gen(rnd, ur_a);
 173 boost::uniform_real<float> ur_a2(0.0F, 100.0F);
 174 boost::variate_generator<boost::mt19937, boost::uniform_real<float> > gen2(rnd, ur_a2);
 175
 176 struct beta_data_generator
 177 {
 178    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_)
 179    {
 180       float a = truncate_to_float(real_cast<float>(gen() * pow(mp_t(10), ap)));
 181       float b = truncate_to_float(real_cast<float>(gen() * pow(mp_t(10), bp)));
 182       float x = truncate_to_float(real_cast<float>(x_));
 183       std::cout << a << " " << b << " " << x << std::endl;
 184       std::pair<mp_t, mp_t> ib_full = ibeta_fraction1(mp_t(a), mp_t(b), mp_t(x));
 185       std::pair<mp_t, mp_t> ib_reg = ibeta_fraction1_regular(mp_t(a), mp_t(b), mp_t(x));
 186       return boost::math::make_tuple(a, b, x, ib_full.first, ib_full.second, ib_reg.first, ib_reg.second);
 187    }
 188 };
 189
 190 // medium sized values:
 191 struct beta_data_generator_medium
 192 {
 193    boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t, mp_t, mp_t> operator()(mp_t x_)
 194    {
 195       mp_t a = gen2();
 196       mp_t b = gen2();
 197       mp_t x = x_;
 198       a = ConvPrec(a, 22);
 199       b = ConvPrec(b, 22);
 200       x = ConvPrec(x, 22);
 201       std::cout << a << " " << b << " " << x << std::endl;
 202       //mp_t exp_beta = boost::math::beta(a, b, x);
 203       std::pair<mp_t, mp_t> ib_full = ibeta_fraction1(mp_t(a), mp_t(b), mp_t(x));
 204       /*exp_beta = boost::math::tools::relative_error(ib_full.first, exp_beta);
 205       if(exp_beta > 1e-40)
 206       {
 207          std::cout << exp_beta << std::endl;
 208       }*/
 209       std::pair<mp_t, mp_t> ib_reg = ibeta_fraction1_regular(mp_t(a), mp_t(b), mp_t(x));
 210       return boost::math::make_tuple(a, b, x, ib_full.first, ib_full.second, ib_reg.first, ib_reg.second);
 211    }
 212 };
 213
 214 struct beta_data_generator_small
 215 {
 216    boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t, mp_t, mp_t> operator()(mp_t x_)
 217    {
 218       float a = truncate_to_float(gen2()/10);
 219       float b = truncate_to_float(gen2()/10);
 220       float x = truncate_to_float(real_cast<float>(x_));
 221       std::cout << a << " " << b << " " << x << std::endl;
 222       std::pair<mp_t, mp_t> ib_full = ibeta_fraction1(mp_t(a), mp_t(b), mp_t(x));
 223       std::pair<mp_t, mp_t> ib_reg = ibeta_fraction1_regular(mp_t(a), mp_t(b), mp_t(x));
 224       return boost::math::make_tuple(a, b, x, ib_full.first, ib_full.second, ib_reg.first, ib_reg.second);
 225    }
 226 };
 227
 228 struct beta_data_generator_int
 229 {
 230    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_)
 231    {
 232       float x = truncate_to_float(real_cast<float>(x_));
 233       std::cout << a << " " << b << " " << x << std::endl;
 234       std::pair<mp_t, mp_t> ib_full = ibeta_fraction1(a, b, mp_t(x));
 235       std::pair<mp_t, mp_t> ib_reg = ibeta_fraction1_regular(a, b, mp_t(x));
 236       return boost::math::make_tuple(a, b, x, ib_full.first, ib_full.second, ib_reg.first, ib_reg.second);
 237    }
 238 };
 239
 240 int main(int, char* [])
 241 {
 242    parameter_info<mp_t> arg1, arg2, arg3, arg4, arg5;
 243    test_data<mp_t> data;
 244
 245    std::cout << "Welcome.\n"
 246       "This program will generate spot tests for the incomplete beta functions:\n"
 247       "  beta(a, b, x) and ibeta(a, b, x)\n\n"
 248       "This is not an interactive program be prepared for a long wait!!!\n\n";
 249
 250    arg1 = make_periodic_param(mp_t(-5), mp_t(6), 11);
 251    arg2 = make_periodic_param(mp_t(-5), mp_t(6), 11);
 252    arg3 = make_random_param(mp_t(0.0001), mp_t(1), 10);
 253    arg4 = make_random_param(mp_t(0.0001), mp_t(1), 100 /*500*/);
 254    arg5 = make_periodic_param(mp_t(1), mp_t(41), 10);
 255
 256    arg1.type |= dummy_param;
 257    arg2.type |= dummy_param;
 258    arg3.type |= dummy_param;
 259    arg4.type |= dummy_param;
 260    arg5.type |= dummy_param;
 261
 262    // comment out all but one of the following when running
 263    // or this program will take forever to complete!
 264    //data.insert(beta_data_generator(), arg1, arg2, arg3);
 265    //data.insert(beta_data_generator_medium(), arg4);
 266    //data.insert(beta_data_generator_small(), arg4);
 267    data.insert(beta_data_generator_int(), arg5, arg5, arg3);
 268
 269    test_data<mp_t>::const_iterator i, j;
 270    i = data.begin();
 271    j = data.end();
 272    while(i != j)
 273    {
 274       mp_t v1 = beta((*i)[0], (*i)[1], (*i)[2]);
 275       mp_t v2 = relative_error(v1, (*i)[3]);
 276       std::string s = boost::lexical_cast<std::string>((*i)[3]);
 277       mp_t v3 = boost::lexical_cast<mp_t>(s);
 278       mp_t v4 = relative_error(v3, (*i)[3]);
 279       if(v2 > 1e-40)
 280       {
 281          std::cout << v2 << std::endl;
 282       }
 283       if(v4 > 1e-60)
 284       {
 285          std::cout << v4 << std::endl;
 286       }
 287       ++ i;
 288    }
 289
 290    std::cout << "Enter name of test data file [default=ibeta_data.ipp]";
 291    std::string line;
 292    std::getline(std::cin, line);
 293    boost::algorithm::trim(line);
 294    if(line == "")
 295       line = "ibeta_data.ipp";
 296    std::ofstream ofs(line.c_str());
 297    ofs << std::scientific << std::setprecision(40);
 298    write_code(ofs, data, "ibeta_data");
 299
 300    return 0;
 301 }
 302
 303