]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
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 | ||
f67539c2 | 53 | std::cout << "Enter range for parameter a: "; |
92f5a8d4 | 54 | std::cin >> a_start >> a_end; |
f67539c2 | 55 | std::cout << "Enter range for parameter b: "; |
92f5a8d4 TL |
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 |