ceph/src/boost/libs/math/test/test_roots.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 <pch.hpp>
   7
   8 #define BOOST_TEST_MAIN
   9 #include <boost/test/unit_test.hpp>
  10 #include <boost/test/floating_point_comparison.hpp>
  11 #include <boost/test/results_collector.hpp>
  12 #include <boost/math/special_functions/beta.hpp>
  13 #include <boost/math/tools/roots.hpp>
  14 #include <boost/test/results_collector.hpp>
  15 #include <boost/test/unit_test.hpp>
  16 #include <boost/array.hpp>
  17 #include "table_type.hpp"
  18 #include <iostream>
  19 #include <iomanip>
  20
  21 #define BOOST_CHECK_CLOSE_EX(a, b, prec, i) \
  22    {\
  23       unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
  24       BOOST_CHECK_CLOSE(a, b, prec); \
  25       if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\
  26       {\
  27          std::cerr << "Failure was at row " << i << std::endl;\
  28          std::cerr << std::setprecision(35); \
  29          std::cerr << "{ " << data[i][0] << " , " << data[i][1] << " , " << data[i][2];\
  30          std::cerr << " , " << data[i][3] << " , " << data[i][4] << " , " << data[i][5] << " } " << std::endl;\
  31       }\
  32    }
  33
  34
  35 //
  36 // Implement various versions of inverse of the incomplete beta
  37 // using different root finding algorithms, and deliberately "bad"
  38 // starting conditions: that way we get all the pathological cases
  39 // we could ever wish for!!!
  40 //
  41
  42 template <class T, class Policy>
  43 struct ibeta_roots_1   // for first order algorithms
  44 {
  45    ibeta_roots_1(T _a, T _b, T t, bool inv = false)
  46       : a(_a), b(_b), target(t), invert(inv) {}
  47
  48    T operator()(const T& x)
  49    {
  50       return boost::math::detail::ibeta_imp(a, b, x, Policy(), invert, true) - target;
  51    }
  52 private:
  53    T a, b, target;
  54    bool invert;
  55 };
  56
  57 template <class T, class Policy>
  58 struct ibeta_roots_2   // for second order algorithms
  59 {
  60    ibeta_roots_2(T _a, T _b, T t, bool inv = false)
  61       : a(_a), b(_b), target(t), invert(inv) {}
  62
  63    boost::math::tuple<T, T> operator()(const T& x)
  64    {
  65       typedef typename boost::math::lanczos::lanczos<T, Policy>::type L;
  66       T f = boost::math::detail::ibeta_imp(a, b, x, Policy(), invert, true) - target;
  67       T f1 = invert ?
  68          -boost::math::detail::ibeta_power_terms(b, a, 1 - x, x, L(), true, Policy())
  69                : boost::math::detail::ibeta_power_terms(a, b, x, 1 - x, L(), true, Policy());
  70       T y = 1 - x;
  71       if(y == 0)
  72          y = boost::math::tools::min_value<T>() * 8;
  73       f1 /= y * x;
  74
  75       // make sure we don't have a zero derivative:
  76       if(f1 == 0)
  77          f1 = (invert ? -1 : 1) * boost::math::tools::min_value<T>() * 64;
  78
  79       return boost::math::make_tuple(f, f1);
  80    }
  81 private:
  82    T a, b, target;
  83    bool invert;
  84 };
  85
  86 template <class T, class Policy>
  87 struct ibeta_roots_3   // for third order algorithms
  88 {
  89    ibeta_roots_3(T _a, T _b, T t, bool inv = false)
  90       : a(_a), b(_b), target(t), invert(inv) {}
  91
  92    boost::math::tuple<T, T, T> operator()(const T& x)
  93    {
  94       typedef typename boost::math::lanczos::lanczos<T, Policy>::type L;
  95       T f = boost::math::detail::ibeta_imp(a, b, x, Policy(), invert, true) - target;
  96       T f1 = invert ?
  97                -boost::math::detail::ibeta_power_terms(b, a, 1 - x, x, L(), true, Policy())
  98                : boost::math::detail::ibeta_power_terms(a, b, x, 1 - x, L(), true, Policy());
  99       T y = 1 - x;
 100       if(y == 0)
 101          y = boost::math::tools::min_value<T>() * 8;
 102       f1 /= y * x;
 103       T f2 = f1 * (-y * a + (b - 2) * x + 1) / (y * x);
 104       if(invert)
 105          f2 = -f2;
 106
 107       // make sure we don't have a zero derivative:
 108       if(f1 == 0)
 109          f1 = (invert ? -1 : 1) * boost::math::tools::min_value<T>() * 64;
 110
 111       return boost::math::make_tuple(f, f1, f2);
 112    }
 113 private:
 114    T a, b, target;
 115    bool invert;
 116 };
 117
 118 double inverse_ibeta_bisect(double a, double b, double z)
 119 {
 120    typedef boost::math::policies::policy<> pol;
 121    bool invert = false;
 122    int bits = std::numeric_limits<double>::digits;
 123
 124    //
 125    // special cases, we need to have these because there may be other
 126    // possible answers:
 127    //
 128    if(z == 1) return 1;
 129    if(z == 0) return 0;
 130
 131    //
 132    // We need a good estimate of the error in the incomplete beta function
 133    // so that we don't set the desired precision too high.  Assume that 3-bits
 134    // are lost each time the arguments increase by a factor of 10:
 135    //
 136    using namespace std;
 137    int bits_lost = static_cast<int>(ceil(log10((std::max)(a, b)) * 3));
 138    if(bits_lost < 0)
 139       bits_lost = 3;
 140    else
 141       bits_lost += 3;
 142    int precision = bits - bits_lost;
 143
 144    double min = 0;
 145    double max = 1;
 146    boost::math::tools::eps_tolerance<double> tol(precision);
 147    return boost::math::tools::bisect(ibeta_roots_1<double, pol>(a, b, z, invert), min, max, tol).first;
 148 }
 149
 150 double inverse_ibeta_newton(double a, double b, double z)
 151 {
 152    double guess = 0.5;
 153    bool invert = false;
 154    int bits = std::numeric_limits<double>::digits;
 155
 156    //
 157    // special cases, we need to have these because there may be other
 158    // possible answers:
 159    //
 160    if(z == 1) return 1;
 161    if(z == 0) return 0;
 162
 163    //
 164    // We need a good estimate of the error in the incomplete beta function
 165    // so that we don't set the desired precision too high.  Assume that 3-bits
 166    // are lost each time the arguments increase by a factor of 10:
 167    //
 168    using namespace std;
 169    int bits_lost = static_cast<int>(ceil(log10((std::max)(a, b)) * 3));
 170    if(bits_lost < 0)
 171       bits_lost = 3;
 172    else
 173       bits_lost += 3;
 174    int precision = bits - bits_lost;
 175
 176    double min = 0;
 177    double max = 1;
 178    return boost::math::tools::newton_raphson_iterate(ibeta_roots_2<double, boost::math::policies::policy<> >(a, b, z, invert), guess, min, max, precision);
 179 }
 180
 181 double inverse_ibeta_halley(double a, double b, double z)
 182 {
 183    double guess = 0.5;
 184    bool invert = false;
 185    int bits = std::numeric_limits<double>::digits;
 186
 187    //
 188    // special cases, we need to have these because there may be other
 189    // possible answers:
 190    //
 191    if(z == 1) return 1;
 192    if(z == 0) return 0;
 193
 194    //
 195    // We need a good estimate of the error in the incomplete beta function
 196    // so that we don't set the desired precision too high.  Assume that 3-bits
 197    // are lost each time the arguments increase by a factor of 10:
 198    //
 199    using namespace std;
 200    int bits_lost = static_cast<int>(ceil(log10((std::max)(a, b)) * 3));
 201    if(bits_lost < 0)
 202       bits_lost = 3;
 203    else
 204       bits_lost += 3;
 205    int precision = bits - bits_lost;
 206
 207    double min = 0;
 208    double max = 1;
 209    return boost::math::tools::halley_iterate(ibeta_roots_3<double, boost::math::policies::policy<> >(a, b, z, invert), guess, min, max, precision);
 210 }
 211
 212 double inverse_ibeta_schroder(double a, double b, double z)
 213 {
 214    double guess = 0.5;
 215    bool invert = false;
 216    int bits = std::numeric_limits<double>::digits;
 217
 218    //
 219    // special cases, we need to have these because there may be other
 220    // possible answers:
 221    //
 222    if(z == 1) return 1;
 223    if(z == 0) return 0;
 224
 225    //
 226    // We need a good estimate of the error in the incomplete beta function
 227    // so that we don't set the desired precision too high.  Assume that 3-bits
 228    // are lost each time the arguments increase by a factor of 10:
 229    //
 230    using namespace std;
 231    int bits_lost = static_cast<int>(ceil(log10((std::max)(a, b)) * 3));
 232    if(bits_lost < 0)
 233       bits_lost = 3;
 234    else
 235       bits_lost += 3;
 236    int precision = bits - bits_lost;
 237
 238    double min = 0;
 239    double max = 1;
 240    return boost::math::tools::schroder_iterate(ibeta_roots_3<double, boost::math::policies::policy<> >(a, b, z, invert), guess, min, max, precision);
 241 }
 242
 243
 244 template <class Real, class T>
 245 void test_inverses(const T& data)
 246 {
 247    using namespace std;
 248    typedef Real                   value_type;
 249
 250    value_type precision = static_cast<value_type>(ldexp(1.0, 1-boost::math::policies::digits<value_type, boost::math::policies::policy<> >()/2)) * 150;
 251    if(boost::math::policies::digits<value_type, boost::math::policies::policy<> >() < 50)
 252       precision = 1;   // 1% or two decimal digits, all we can hope for when the input is truncated
 253
 254    for(unsigned i = 0; i < data.size(); ++i)
 255    {
 256       //
 257       // These inverse tests are thrown off if the output of the
 258       // incomplete beta is too close to 1: basically there is insuffient
 259       // information left in the value we're using as input to the inverse
 260       // to be able to get back to the original value.
 261       //
 262       if(data[i][5] == 0)
 263       {
 264          BOOST_CHECK_EQUAL(inverse_ibeta_halley(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(0));
 265          BOOST_CHECK_EQUAL(inverse_ibeta_schroder(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(0));
 266          BOOST_CHECK_EQUAL(inverse_ibeta_newton(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(0));
 267          BOOST_CHECK_EQUAL(inverse_ibeta_bisect(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(0));
 268       }
 269       else if((1 - data[i][5] > 0.001)
 270          && (fabs(data[i][5]) > 2 * boost::math::tools::min_value<value_type>())
 271          && (fabs(data[i][5]) > 2 * boost::math::tools::min_value<double>()))
 272       {
 273          value_type inv = inverse_ibeta_halley(Real(data[i][0]), Real(data[i][1]), Real(data[i][5]));
 274          BOOST_CHECK_CLOSE_EX(Real(data[i][2]), inv, precision, i);
 275          inv = inverse_ibeta_schroder(Real(data[i][0]), Real(data[i][1]), Real(data[i][5]));
 276          BOOST_CHECK_CLOSE_EX(Real(data[i][2]), inv, precision, i);
 277          inv = inverse_ibeta_newton(Real(data[i][0]), Real(data[i][1]), Real(data[i][5]));
 278          BOOST_CHECK_CLOSE_EX(Real(data[i][2]), inv, precision, i);
 279          inv = inverse_ibeta_bisect(Real(data[i][0]), Real(data[i][1]), Real(data[i][5]));
 280          BOOST_CHECK_CLOSE_EX(Real(data[i][2]), inv, precision, i);
 281       }
 282       else if(1 == data[i][5])
 283       {
 284          BOOST_CHECK_EQUAL(inverse_ibeta_halley(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(1));
 285          BOOST_CHECK_EQUAL(inverse_ibeta_schroder(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(1));
 286          BOOST_CHECK_EQUAL(inverse_ibeta_newton(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(1));
 287          BOOST_CHECK_EQUAL(inverse_ibeta_bisect(Real(data[i][0]), Real(data[i][1]), Real(data[i][5])), value_type(1));
 288       }
 289
 290    }
 291 }
 292
 293 #ifndef SC_
 294 #define SC_(x) static_cast<typename table_type<T>::type>(BOOST_JOIN(x, L))
 295 #endif
 296
 297 template <class T>
 298 void test_beta(T, const char* /* name */)
 299 {
 300    //
 301    // The actual test data is rather verbose, so it's in a separate file
 302    //
 303    // The contents are as follows, each row of data contains
 304    // five items, input value a, input value b, integration limits x, beta(a, b, x) and ibeta(a, b, x):
 305    //
 306 #  include "ibeta_small_data.ipp"
 307
 308    test_inverses<T>(ibeta_small_data);
 309
 310 #  include "ibeta_data.ipp"
 311
 312    test_inverses<T>(ibeta_data);
 313
 314 #  include "ibeta_large_data.ipp"
 315
 316    test_inverses<T>(ibeta_large_data);
 317 }
 318
 319 BOOST_AUTO_TEST_CASE( test_main )
 320 {
 321    test_beta(0.1, "double");
 322
 323 }