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)
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
15 #include <boost/multiprecision/mpfr.hpp>
16 #include <boost/math/special_functions/hypergeometric_1F1.hpp>
17 #include <boost/math/special_functions/hypergeometric_pFq.hpp>
18 #include <boost/math/special_functions/relative_difference.hpp>
20 #include <boost/random.hpp>
23 #include <boost/iostreams/tee.hpp>
24 #include <boost/iostreams/stream.hpp>
26 typedef double test_type
;
27 using boost::multiprecision::mpfr_float
;
31 // We convert assertions into exceptions, so we can log them and continue:
33 void assertion_failed(char const * expr
, char const *, char const * file
, long line
)
35 std::ostringstream oss
;
36 oss
<< file
<< ":" << line
<< " Assertion failed: " << expr
;
37 throw std::runtime_error(oss
.str());
42 void print_value(double x
, std::ostream
& os
= std::cout
)
45 double m
= std::frexp(x
, &e
);
46 m
= std::ldexp(m
, 54);
48 std::int64_t val
= (std::int64_t)m
;
49 BOOST_MATH_ASSERT(std::ldexp((double)val
, e
) == x
);
50 os
<< "std::ldexp((double)" << val
<< ", " << e
<< ")";
54 void print_row(double a
, double b
, double z
, mpfr_float result
, std::ostream
& os
= std::cout
)
62 os
<< ", SC_(" << std::setprecision(45) << result
<< ") }}" << std::endl
;
67 error_data(double a
, double b
, double z
, std::intmax_t e
)
68 : a(a
), b(b
), z(z
), error(e
) {}
71 bool operator<(const error_data
& other
)const
73 return error
< other
.error
;
80 test_type max_a
, max_b
, max_z
, min_a
, min_b
, min_z
;
82 unsigned number_of_samples
;
84 std::ofstream log_stream
, incalculable_stream
, unevaluated_stream
, bins_stream
;
87 std::cout
<< "Enter range for a: ";
88 std::cin
>> min_a
>> max_a
;
89 std::cout
<< "Enter range for b: ";
90 std::cin
>> min_b
>> max_b
;
91 std::cout
<< "Enter range for z: ";
92 std::cin
>> min_z
>> max_z
;
93 std::cout
<< "Enter number of samples: ";
94 std::cin
>> number_of_samples
;
95 std::cout
<< "Enter basename for log files: ";
98 typedef boost::iostreams::tee_device
<std::ostream
, std::ostream
> tee_sink
;
99 typedef boost::iostreams::stream
<tee_sink
> tee_stream
;
101 log_stream
.open((basename
+ ".log").c_str());
102 tee_stream
tee_log(tee_sink(std::cout
, log_stream
));
103 incalculable_stream
.open((basename
+ "_incalculable.log").c_str());
104 unevaluated_stream
.open((basename
+ "_unevaluated.log").c_str());
105 bins_stream
.open((basename
+ "_bins.csv").c_str());
106 tee_stream
tee_bins(tee_sink(std::cout
, bins_stream
));
108 boost::random::mt19937
gen(std::time(0));
109 boost::random::uniform_real_distribution
<test_type
> a_dist(min_a
, max_a
);
110 boost::random::uniform_real_distribution
<test_type
> b_dist(min_b
, max_b
);
111 boost::random::uniform_real_distribution
<test_type
> z_dist(min_z
, max_z
);
113 std::multiset
<error_data
> errors
;
114 std::map
<std::pair
<int, int>, int> bins
;
116 unsigned incalculable
= 0;
117 unsigned evaluation_errors
= 0;
118 test_type max_error
= 0;
122 test_type a
= a_dist(gen
);
123 test_type b
= b_dist(gen
);
124 test_type z
= z_dist(gen
);
125 test_type found
, expected
;
126 mpfr_float mp_expected
;
129 mp_expected
= boost::math::hypergeometric_pFq_precision({ mpfr_float(a
) }, { mpfr_float(b
) }, mpfr_float(z
), 25, 200.0);
130 expected
= (test_type
)mp_expected
;
132 catch (const std::exception
&)
134 // Unable to compute reference value:
136 tee_log
<< "Unable to compute reference value in reasonable time: " << std::endl
;
137 print_row(a
, b
, z
, mpfr_float(0), tee_log
);
138 incalculable_stream
<< std::setprecision(6) << std::scientific
<< a
<< "," << b
<< "," << z
<< "\n";
143 found
= boost::math::hypergeometric_1F1(a
, b
, z
);
145 catch (const std::exception
&)
149 log_stream
<< "Unexpected exception calculating value: " << std::endl
;
150 print_row(a
, b
, z
, mp_expected
, log_stream
);
151 unevaluated_stream
<< std::setprecision(6) << std::scientific
<< a
<< "," << b
<< "," << z
<< "\n";
154 test_type err
= boost::math::epsilon_difference(found
, expected
);
157 tee_log
<< "New maximum error is: " << err
<< std::endl
;
158 print_row(a
, b
, z
, mp_expected
, tee_log
);
162 errors
.insert(error_data(a
, b
, z
, boost::math::lltrunc(err
)));
166 errors
.insert(error_data(a
, b
, z
, INT_MAX
));
169 if (number_of_samples
% 500 == 0)
170 std::cout
<< number_of_samples
<< " samples to go" << std::endl
;
171 } while (number_of_samples
);
173 tee_log
<< "Max error found was: " << max_error
<< std::endl
;
175 unsigned current_bin
= 0;
177 unsigned old_lim
= 0;
179 while (errors
.size())
183 //std::cout << "Enter upper limit for bin " << current_bin << ": ";
185 auto p
= errors
.upper_bound(error_data(0, 0, 0, lim
));
186 int bin_count
= std::distance(errors
.begin(), p
);
189 std::ofstream
os((basename
+ "_errors_" + std::to_string(current_bin
+ 1) + ".csv").c_str());
190 os
<< "a,b,z,error\n";
191 bins
[std::make_pair(old_lim
, lim
)] = bin_count
;
192 for (auto pos
= errors
.begin(); pos
!= p
; ++pos
)
194 os
<< pos
->a
<< "," << pos
->b
<< "," << pos
->z
<< "," << pos
->error
<< "\n";
196 errors
.erase(errors
.begin(), p
);
201 tee_bins
<< "Results:\n\n";
202 tee_bins
<< "#bin,Range,2^N,Count\n";
204 for (auto p
= bins
.begin(); p
!= bins
.end(); ++p
, ++hash
)
206 tee_bins
<< hash
<< "," << p
->first
.first
<< "-" << p
->first
.second
<< "," << hash
+1 << "," << p
->second
<< std::endl
;
208 if (evaluation_errors
)
210 tee_bins
<< ",Failed,," << evaluation_errors
<< std::endl
;
214 tee_bins
<< ",Incalculable,," << incalculable
<< std::endl
;
217 catch (const std::exception
& e
)
219 std::cout
<< "Terminating with unhandled exception: " << e
.what() << std::endl
;