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

   1 //  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 #define BOOST_ENABLE_ASSERT_HANDLER
   7 #define BOOST_MATH_MAX_SERIES_ITERATION_POLICY INT_MAX
   8 // for consistent behaviour across compilers/platforms:
   9 #define BOOST_MATH_PROMOTE_DOUBLE_POLICY false
  10 // overflow to infinity is OK, we treat these as zero error as long as the sign is correct!
  11 #define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
  12
  13 #include <iostream>
  14 #include <ctime>
  15 #include <boost/multiprecision/mpfr.hpp>
  16 #include <boost/multiprecision/cpp_bin_float.hpp>
  17 #include <boost/math/special_functions/hypergeometric_1F1.hpp>
  18 #include <boost/math/special_functions/hypergeometric_pFq.hpp>
  19 #include <boost/math/special_functions/relative_difference.hpp>
  20
  21 #include <boost/random.hpp>
  22 #include <set>
  23 #include <fstream>
  24 #include <boost/iostreams/tee.hpp>
  25 #include <boost/iostreams/stream.hpp>
  26
  27 using boost::multiprecision::mpfr_float;
  28
  29 namespace boost {
  30    //
  31    // We convert assertions into exceptions, so we can log them and continue:
  32    //
  33    void assertion_failed(char const * expr, char const *, char const * file, long line)
  34    {
  35       std::ostringstream oss;
  36       oss << file << ":" << line << " Assertion failed: " << expr;
  37       throw std::runtime_error(oss.str());
  38    }
  39
  40 }
  41
  42 typedef boost::multiprecision::cpp_bin_float_quad test_type;
  43
  44 int main()
  45 {
  46    using std::floor;
  47    using std::ceil;
  48    try {
  49       test_type a_start, a_end;
  50       test_type b_start, b_end;
  51       test_type a_mult, b_mult;
  52
  53       std::cout << "Enter range for paramater a: ";
  54       std::cin >> a_start >> a_end;
  55       std::cout << "Enter range for paramater b: ";
  56       std::cin >> b_start >> b_end;
  57       std::cout << "Enter multiplier for a parameter: ";
  58       std::cin >> a_mult;
  59       std::cout << "Enter multiplier for b parameter: ";
  60       std::cin >> b_mult;
  61
  62       double error_limit = 200;
  63       double time_limit = 10.0;
  64
  65       for (test_type a = a_start; a < a_end; a_start < 0 ? a /= a_mult : a *= a_mult)
  66       {
  67          for (test_type b = b_start; b < b_end; b_start < 0 ? b /= b_mult : b *= b_mult)
  68          {
  69             test_type z_mult = 2;
  70             test_type last_good = 0;
  71             test_type bad = 0;
  72             try {
  73                for (test_type z = 1; z < 1e10; z *= z_mult, z_mult *= 2)
  74                {
  75                   // std::cout << "z = " << z << std::endl;
  76                   boost::uintmax_t max_iter = 1000;
  77                   test_type calc = boost::math::tools::function_ratio_from_forwards_recurrence(boost::math::detail::hypergeometric_1F1_recurrence_a_and_b_coefficients<test_type>(a, b, z), std::numeric_limits<test_type>::epsilon() * 2, max_iter);
  78                   test_type reference = (test_type)(boost::math::hypergeometric_pFq_precision({ mpfr_float(a) }, { mpfr_float(b) }, mpfr_float(z), 50, time_limit) / boost::math::hypergeometric_pFq_precision({ mpfr_float(a + 1) }, { mpfr_float(b + 1) }, mpfr_float(z), std::numeric_limits<test_type>::digits10 * 2, time_limit));
  79                   double err = (double)boost::math::epsilon_difference(reference, calc);
  80
  81                   if (err < error_limit)
  82                   {
  83                      last_good = z;
  84                      break;
  85                   }
  86                   else
  87                   {
  88                      bad = z;
  89                   }
  90                }
  91             }
  92             catch (const std::exception& e)
  93             {
  94                std::cout << "Unexpected exception: " << e.what() << std::endl;
  95                std::cout << "For a = " << a << " b = " << b << " z = " << bad * z_mult / 2 << std::endl;
  96             }
  97             test_type z_limit;
  98             if (0 == bad)
  99                z_limit = 1;  // Any z is large enough
 100             else if (0 == last_good)
 101                z_limit = std::numeric_limits<test_type > ::infinity();
 102             else
 103             {
 104                //
 105                // At this stage last_good and bad should bracket the edge of the domain, bisect to narrow things down:
 106                //
 107                z_limit = last_good == 0 ? 0 : boost::math::tools::bisect([&a, b, error_limit, time_limit](test_type z)
 108                {
 109                   boost::uintmax_t max_iter = 1000;
 110                   test_type calc = boost::math::tools::function_ratio_from_forwards_recurrence(boost::math::detail::hypergeometric_1F1_recurrence_a_and_b_coefficients<test_type>(a, b, z), std::numeric_limits<test_type>::epsilon() * 2, max_iter);
 111                   test_type reference = (test_type)(boost::math::hypergeometric_pFq_precision({ mpfr_float(a) }, { mpfr_float(b) }, mpfr_float(z), 50, time_limit + 20) / boost::math::hypergeometric_pFq_precision({ mpfr_float(a + 1) }, { mpfr_float(b + 1) }, mpfr_float(z), std::numeric_limits<test_type>::digits10 * 2, time_limit + 20));
 112                   test_type err = boost::math::epsilon_difference(reference, calc);
 113                   return err < error_limit ? 1 : -1;
 114                }, bad, last_good, boost::math::tools::equal_floor()).first;
 115                z_limit = floor(z_limit + 2);  // Give ourselves some headroom!
 116             }
 117             // std::cout << "z_limit = " << z_limit << std::endl;
 118             //
 119             // Now over again for backwards recurrence domain at the same points:
 120             //
 121             bad = z_limit > 1e10 ? 1e10 : z_limit;
 122             last_good = 0;
 123             z_mult = 1.1;
 124             for (test_type z = bad; z > 1; z /= z_mult, z_mult *= 2)
 125             {
 126                // std::cout << "z = " << z << std::endl;
 127                try {
 128                   boost::uintmax_t max_iter = 1000;
 129                   test_type calc = boost::math::tools::function_ratio_from_backwards_recurrence(boost::math::detail::hypergeometric_1F1_recurrence_a_and_b_coefficients<test_type>(a, b, z), std::numeric_limits<test_type>::epsilon() * 2, max_iter);
 130                   test_type reference = (test_type)(boost::math::hypergeometric_pFq_precision({ mpfr_float(a) }, { mpfr_float(b) }, mpfr_float(z), 50, time_limit) / boost::math::hypergeometric_pFq_precision({ mpfr_float(a - 1) }, { mpfr_float(b - 1) }, mpfr_float(z), std::numeric_limits<test_type>::digits10 * 2, time_limit));
 131                   test_type err = boost::math::epsilon_difference(reference, calc);
 132
 133                   if (err < error_limit)
 134                   {
 135                      last_good = z;
 136                      break;
 137                   }
 138                   else
 139                   {
 140                      bad = z;
 141                   }
 142                }
 143                catch (const std::exception& e)
 144                {
 145                   bad = z;
 146                   std::cout << "Unexpected exception: " << e.what() << std::endl;
 147                   std::cout << "For a = " << a << " b = " << b << " z = " << z << std::endl;
 148                }
 149             }
 150             test_type lower_z_limit;
 151             if (last_good < 1)
 152                lower_z_limit = 0;
 153             else if (last_good >= bad)
 154             {
 155                boost::uintmax_t max_iter = 1000;
 156                test_type z = bad;
 157                test_type calc = boost::math::tools::function_ratio_from_forwards_recurrence(boost::math::detail::hypergeometric_1F1_recurrence_a_and_b_coefficients<test_type>(a, b, z), std::numeric_limits<test_type>::epsilon() * 2, max_iter);
 158                test_type reference = (test_type)(boost::math::hypergeometric_pFq_precision({ mpfr_float(a) }, { mpfr_float(b) }, mpfr_float(z), 50, time_limit) / boost::math::hypergeometric_pFq_precision({ mpfr_float(a + 1) }, { mpfr_float(b + 1) }, mpfr_float(z), std::numeric_limits<test_type>::digits10 * 2, time_limit));
 159                test_type err = boost::math::epsilon_difference(reference, calc);
 160                if (err < error_limit)
 161                {
 162                   lower_z_limit = bad;   //  Both forwards and backwards iteration work!!!
 163                }
 164                else
 165                   throw std::runtime_error("Internal logic failed!");
 166             }
 167             else
 168             {
 169                //
 170                // At this stage last_good and bad should bracket the edge of the domain, bisect to narrow things down:
 171                //
 172                lower_z_limit = last_good == 0 ? 0 : boost::math::tools::bisect([&a, b, error_limit, time_limit](test_type z)
 173                {
 174                   boost::uintmax_t max_iter = 1000;
 175                   test_type calc = boost::math::tools::function_ratio_from_backwards_recurrence(boost::math::detail::hypergeometric_1F1_recurrence_a_and_b_coefficients<test_type>(a, b, z), std::numeric_limits<test_type>::epsilon() * 2, max_iter);
 176                   test_type reference = (test_type)(boost::math::hypergeometric_pFq_precision({ mpfr_float(a) }, { mpfr_float(b) }, mpfr_float(z), 50, time_limit + 20) / boost::math::hypergeometric_pFq_precision({ mpfr_float(a - 1) }, { mpfr_float(b - 1) }, mpfr_float(z), std::numeric_limits<test_type>::digits10 * 2, time_limit + 20));
 177                   test_type err = boost::math::epsilon_difference(reference, calc);
 178                   return err < error_limit ? 1 : -1;
 179                }, last_good, bad, boost::math::tools::equal_floor()).first;
 180                z_limit = ceil(z_limit - 2);  // Give ourselves some headroom!
 181             }
 182
 183             std::cout << std::setprecision(std::numeric_limits<test_type>::max_digits10) << "{ " << a << ", " << b << ", " << lower_z_limit << ", " << z_limit << "}," << std::endl;
 184          }
 185       }
 186    }
 187    catch (const std::exception& e)
 188    {
 189       std::cout << "Unexpected exception: " << e.what() << std::endl;
 190    }
 191    return 0;
 192 }
 193