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)
6 #define BOOST_MATH_MAX_SERIES_ITERATION_POLICY 10000000
8 #define BOOST_MATH_USE_MPFR
10 #include <boost/math/special_functions/hypergeometric_1f1.hpp>
11 #include <boost/math/constants/constants.hpp>
12 #include <boost/lexical_cast.hpp>
15 #include <boost/math/tools/test_data.hpp>
16 #include <boost/random.hpp>
18 #include <boost/multiprecision/mpfi.hpp>
20 using namespace boost::math::tools
;
21 using namespace boost::math
;
23 using namespace boost::multiprecision
;
25 typedef mpfi_float_1000 mpfi_type
;
27 mp_t
hypergeometric_1f1_generic_series(mp_t a_
, mp_t b_
, mp_t z_
)
29 using namespace boost::math::tools
;
30 using namespace boost::math
;
32 using namespace boost::multiprecision
;
34 mpfi_type
a(a_
), b(b_
), z(z_
), sum(0), term(1), diff
, term0(0);
40 max_n
= itrunc(-b
) + 10000;
44 mpfi_type
overflow_limit("1.189731495357231765e+4900"); // a bit less than LDBL_MAX for extended long doubles.
49 term
*= (((a
+ n
) / ((b
+ n
) * (n
+ 1))) * z
);
51 diff
= fabs(term
/ sum
);
54 std::cout
<< "Aborting series evaluation due to too many iterations...\n";
55 throw evaluation_error("");
57 if (fabs(upper(sum
)) > overflow_limit
)
59 std::cout
<< "Aborting series evaluation due to over large sum...\n";
60 throw evaluation_error("");
62 cont
= (fabs(upper(diff
)) > 1e-40) || (b
+ n
< 0) || (fabs(term0
) < fabs(term
));
64 //std::cout << upper(term) << " " << upper(sum) << " " << upper(diff) << " " << cont << std::endl;
67 mp_t r
= mp_t(width(sum
) / median(sum
));
70 std::cout
<< "Aborting to to error in result of " << r
<< std::endl
;
71 throw evaluation_error("");
73 std::cout
<< "Found error in sum was " << r
<< std::endl
;
75 return mp_t(median(sum
));
79 struct hypergeometric_1f1_gen
81 mp_t
operator()(mp_t a1
, mp_t a2
, mp_t z
)
85 result
= hypergeometric_1f1_generic_series(a1
, a2
, z
);
86 std::cout
<< a1
<< " " << a2
<< " " << z
<< " " << result
<< std::endl
;
90 throw std::domain_error("");
92 if (fabs(result
) > (std::numeric_limits
<double>::max
)())
94 std::cout
<< "Rejecting over large value\n";
95 throw std::domain_error("");
102 int main(int, char* [])
104 parameter_info
<mp_t
> arg1
, arg2
, arg3
;
105 test_data
<mp_t
> data
;
107 std::cout
<< "Welcome.\n"
108 "This program will generate spot tests for 1F1 (Yeh!!):\n";
114 random_ns::mt19937 rnd
;
115 random_ns::uniform_real_distribution
<float> ur_a(0, 1);
150 for (unsigned i
= 0; i
< v
.size(); ++i
)
152 for (unsigned j
= 0; j
< v
.size(); ++j
)
154 for (unsigned k
= 0; k
< v
.size(); ++k
)
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
);
167 std::cout
<< "Enter name of test data file [default=hypergeometric_1f1.ipp]";
168 std::getline(std::cin
, line
);
169 boost::algorithm::trim(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());