ceph/src/boost/libs/math/tools/hyp_1f1_reg_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 #define BOOST_MATH_USE_MPFR
   8 #include "mp_t.hpp"
   9 #include <boost/math/special_functions/hypergeometric_1F1.hpp>
  10 #include <boost/math/special_functions/hypergeometric_pFq.hpp>
  11 #include <boost/math/constants/constants.hpp>
  12 #include <boost/lexical_cast.hpp>
  13 #include <fstream>
  14 #include <map>
  15 #include <boost/math/tools/test_data.hpp>
  16 #include <boost/random.hpp>
  17
  18 #include <boost/multiprecision/mpfr.hpp>
  19
  20 using namespace boost::math::tools;
  21 using namespace boost::math;
  22 using namespace std;
  23 using namespace boost::multiprecision;
  24
  25
  26 struct hypergeometric_1f1_gen
  27 {
  28    mp_t operator()(mp_t a1, mp_t a2, mp_t z)
  29    {
  30       mp_t result;
  31       try {
  32          result = (mp_t)boost::math::hypergeometric_pFq_precision({ mpfr_float(a1) }, { mpfr_float(mpfr_float(a2)) }, mpfr_float(z), 50, 25.0) / boost::multiprecision::tgamma(a2);
  33          std::cout << a1 << " " << a2 << " " << z << " " << result << std::endl;
  34       }
  35       catch (...)
  36       {
  37          throw std::domain_error("");
  38       }
  39       if (fabs(result) > (std::numeric_limits<double>::max)())
  40       {
  41          std::cout << "Rejecting over large value\n";
  42          throw std::domain_error("");
  43       }
  44       return result;
  45    }
  46 };
  47
  48
  49 int main(int, char* [])
  50 {
  51    parameter_info<mp_t> arg1, arg2, arg3;
  52    test_data<mp_t> data;
  53
  54    std::cout << "Welcome.\n"
  55       "This program will generate spot tests for 1F1 (Yeh!!):\n";
  56
  57    std::string line;
  58    //bool cont;
  59
  60    std::vector<mp_t> v;
  61    random_ns::mt19937 rnd;
  62    random_ns::uniform_real_distribution<float> ur_a(0, 1);
  63
  64    mp_t p = ur_a(rnd);
  65    p *= 1e6;
  66    v.push_back(p);
  67    v.push_back(-p);
  68    p = ur_a(rnd);
  69    p *= 1e5;
  70    v.push_back(p);
  71    v.push_back(-p);
  72    p = ur_a(rnd);
  73    p *= 1e4;
  74    v.push_back(p);
  75    v.push_back(-p);
  76    p = ur_a(rnd);
  77    p *= 1e3;
  78    v.push_back(p);
  79    v.push_back(-p);
  80    p = ur_a(rnd);
  81    p *= 1e2;
  82    v.push_back(p);
  83    v.push_back(-p);
  84    p = ur_a(rnd);
  85    p *= 1e-5;
  86    v.push_back(p);
  87    v.push_back(-p);
  88    p = ur_a(rnd);
  89    p *= 1e-12;
  90    v.push_back(p);
  91    v.push_back(-p);
  92    p = ur_a(rnd);
  93    p *= 1e-30;
  94    v.push_back(p);
  95    v.push_back(-p);
  96
  97    for (unsigned i = 0; i < v.size(); ++i)
  98    {
  99       for (unsigned j = 0; j < v.size(); ++j)
 100       {
 101          for (unsigned k = 0; k < v.size(); ++k)
 102          {
 103             std::cout << i << " " << j << " " << k << std::endl;
 104             std::cout << v[i] << " " << (v[j] * 3) / 2 << " " << (v[j] * 5) / 4 << std::endl;
 105             arg1 = make_single_param(v[i]);
 106             arg2 = make_single_param(mp_t((v[j] * 3) / 2));
 107             arg3 = make_single_param(mp_t((v[k] * 5) / 4));
 108             data.insert(hypergeometric_1f1_gen(), arg1, arg2, arg3);
 109          }
 110       }
 111    }
 112
 113
 114    std::cout << "Enter name of test data file [default=hypergeometric_1f1.ipp]";
 115    std::getline(std::cin, line);
 116    boost::algorithm::trim(line);
 117    if(line == "")
 118       line = "hypergeometric_1f1.ipp";
 119    std::ofstream ofs(line.c_str());
 120    ofs << std::scientific << std::setprecision(40);
 121    write_code(ofs, data, line.c_str());
 122
 123    return 0;
 124 }
 125
 126