libs/math/tools/hyp_1f1_big_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 #define BOOST_MATH_MAX_SERIES_ITERATION_POLICY 10000000
   7
   8 #include <boost/math/special_functions/hypergeometric_1f1.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 #define BOOST_MATH_USE_MPFR
  16 #include "mp_t.hpp"
  17
  18 #include <boost/multiprecision/mpfi.hpp>
  19
  20 using namespace boost::math::tools;
  21 using namespace boost::math;
  22 using namespace std;
  23 using namespace boost::multiprecision;
  24
  25 typedef mpfi_float_1000 mpfi_type;
  26
  27 mp_t hypergeometric_1f1_generic_series(mp_t a_, mp_t b_, mp_t z_)
  28 {
  29    using namespace boost::math::tools;
  30    using namespace boost::math;
  31    using namespace std;
  32    using namespace boost::multiprecision;
  33
  34    mpfi_type a(a_), b(b_), z(z_), sum(0), term(1), diff, term0(0);
  35    unsigned n = 0;
  36    bool cont = true;
  37
  38    unsigned max_n;
  39    if (b < 0)
  40       max_n = itrunc(-b) + 10000;
  41    else
  42       max_n = 10000000;
  43
  44    mpfi_type overflow_limit("1.189731495357231765e+4900");  // a bit less than LDBL_MAX for extended long doubles.
  45
  46    do
  47    {
  48       sum += term;
  49       term *= (((a + n) / ((b + n) * (n + 1))) * z);
  50       ++n;
  51       diff = fabs(term / sum);
  52       if (n > max_n)
  53       {
  54          std::cout << "Aborting series evaluation due to too many iterations...\n";
  55          throw evaluation_error("");
  56       }
  57       if (fabs(upper(sum)) > overflow_limit)
  58       {
  59          std::cout << "Aborting series evaluation due to over large sum...\n";
  60          throw evaluation_error("");
  61       }
  62       cont = (fabs(upper(diff)) > 1e-40) || (b + n < 0) || (fabs(term0) < fabs(term));
  63       term0 = term;
  64       //std::cout << upper(term) << " " << upper(sum) << " " << upper(diff) << " " << cont << std::endl;
  65    } while (cont);
  66
  67    mp_t r = mp_t(width(sum) / median(sum));
  68    if (fabs(r) > 1e-40)
  69    {
  70       std::cout << "Aborting to to error in result of " << r << std::endl;
  71       throw evaluation_error("");
  72    }
  73    std::cout << "Found error in sum was " << r << std::endl;
  74
  75    return mp_t(median(sum));
  76 }
  77
  78
  79 struct hypergeometric_1f1_gen
  80 {
  81    mp_t operator()(mp_t a1, mp_t a2, mp_t z)
  82    {
  83       mp_t result;
  84       try {
  85          result = hypergeometric_1f1_generic_series(a1, a2, z);
  86          std::cout << a1 << " " << a2 << " " << z << " " << result << std::endl;
  87       }
  88       catch (...)
  89       {
  90          throw std::domain_error("");
  91       }
  92       if (fabs(result) > (std::numeric_limits<double>::max)())
  93       {
  94          std::cout << "Rejecting over large value\n";
  95          throw std::domain_error("");
  96       }
  97       return result;
  98    }
  99 };
 100
 101
 102 int main(int, char* [])
 103 {
 104    parameter_info<mp_t> arg1, arg2, arg3;
 105    test_data<mp_t> data;
 106
 107    std::cout << "Welcome.\n"
 108       "This program will generate spot tests for 1F1 (Yeh!!):\n";
 109
 110    std::string line;
 111    //bool cont;
 112
 113    std::vector<mp_t> v;
 114    random_ns::mt19937 rnd;
 115    random_ns::uniform_real_distribution<float> ur_a(0, 1);
 116
 117    mp_t p = ur_a(rnd);
 118    p *= 1e6;
 119    v.push_back(p);
 120    v.push_back(-p);
 121    p = ur_a(rnd);
 122    p *= 1e5;
 123    v.push_back(p);
 124    v.push_back(-p);
 125    p = ur_a(rnd);
 126    p *= 1e4;
 127    v.push_back(p);
 128    v.push_back(-p);
 129    p = ur_a(rnd);
 130    p *= 1e3;
 131    v.push_back(p);
 132    v.push_back(-p);
 133    p = ur_a(rnd);
 134    p *= 1e2;
 135    v.push_back(p);
 136    v.push_back(-p);
 137    p = ur_a(rnd);
 138    p *= 1e-5;
 139    v.push_back(p);
 140    v.push_back(-p);
 141    p = ur_a(rnd);
 142    p *= 1e-12;
 143    v.push_back(p);
 144    v.push_back(-p);
 145    p = ur_a(rnd);
 146    p *= 1e-30;
 147    v.push_back(p);
 148    v.push_back(-p);
 149
 150    for (unsigned i = 0; i < v.size(); ++i)
 151    {
 152       for (unsigned j = 0; j < v.size(); ++j)
 153       {
 154          for (unsigned k = 0; k < v.size(); ++k)
 155          {
 156             std::cout << i << " " << j << " " << k << std::endl;
 157             std::cout << v[i] << " " << (v[j] * 3) / 2 << " " << (v[j] * 5) / 4 << std::endl;
 158             arg1 = make_single_param(v[i]);
 159             arg2 = make_single_param(mp_t((v[j] * 3) / 2));
 160             arg3 = make_single_param(mp_t((v[k] * 5) / 4));
 161             data.insert(hypergeometric_1f1_gen(), arg1, arg2, arg3);
 162          }
 163       }
 164    }
 165
 166
 167    std::cout << "Enter name of test data file [default=hypergeometric_1f1.ipp]";
 168    std::getline(std::cin, line);
 169    boost::algorithm::trim(line);
 170    if(line == "")
 171       line = "hypergeometric_1f1.ipp";
 172    std::ofstream ofs(line.c_str());
 173    ofs << std::scientific << std::setprecision(40);
 174    write_code(ofs, data, line.c_str());
 175
 176    return 0;
 177 }
 178
 179