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