1 // Copyright John Maddock 2013.
2
3 // Use, modification and distribution are subject to the
4 // Boost Software License, Version 1.0.
5 // (See accompanying file LICENSE_1_0.txt
6 // or copy at http://www.boost.org/LICENSE_1_0.txt)
7
8 #ifdef _MSC_VER
9 #define _SCL_SECURE_NO_WARNINGS
10 #endif
11
12 #if defined(TEST1) || defined(TEST2) || defined(TEST3) || defined(TEST4)
13 #include <boost/multiprecision/cpp_bin_float.hpp>
14 #include <boost/multiprecision/cpp_int.hpp>
15 #else
16 #include <boost/multiprecision/mpfr.hpp>
17 #endif
18 #include <boost/math/special_functions/next.hpp>
19
20 #include <boost/random/mersenne_twister.hpp>
21 #include <boost/random/uniform_int.hpp>
22 #include <boost/chrono.hpp>
23 #include "test.hpp"
24 #include <boost/array.hpp>
25 #include <iostream>
26 #include <iomanip>
27
28 #ifdef BOOST_MSVC
29 #pragma warning(disable : 4127)
30 #endif
31
32 template <class Clock>
33 struct stopwatch
34 {
35 typedef typename Clock::duration duration;
stopwatchstopwatch36 stopwatch()
37 {
38 m_start = Clock::now();
39 }
elapsedstopwatch40 duration elapsed()
41 {
42 return Clock::now() - m_start;
43 }
resetstopwatch44 void reset()
45 {
46 m_start = Clock::now();
47 }
48
49 private:
50 typename Clock::time_point m_start;
51 };
52
53 template <class T>
54 struct exponent_type
55 {
56 typedef int type;
57 };
58 template <class T, boost::multiprecision::expression_template_option ET>
59 struct exponent_type<boost::multiprecision::number<T, ET> >
60 {
61 typedef typename T::exponent_type type;
62 };
63
64 template <class T>
generate_random_float()65 T generate_random_float()
66 {
67 BOOST_MATH_STD_USING
68 typedef typename exponent_type<T>::type e_type;
69 static boost::random::mt19937 gen;
70 T val = gen();
71 T prev_val = -1;
72 while (val != prev_val)
73 {
74 val *= (gen.max)();
75 prev_val = val;
76 val += gen();
77 }
78 e_type e;
79 val = frexp(val, &e);
80
81 static const int max_exponent_value = (std::min)(static_cast<int>(std::numeric_limits<T>::max_exponent - std::numeric_limits<T>::digits - 20), 2000);
82 static boost::random::uniform_int_distribution<e_type> ui(0, max_exponent_value);
83 return ldexp(val, ui(gen));
84 }
85
86 template <class Float, class Rat>
do_round_trip(const Float & val)87 void do_round_trip(const Float& val)
88 {
89 #ifndef BOOST_MP_NO_CXX11_EXPLICIT_CONVERSION_OPERATORS
90 BOOST_MATH_STD_USING
91 Rat rat(val);
92 Float new_f(rat);
93 BOOST_CHECK_EQUAL(val, new_f);
94 //
95 // Try adding or subtracting an insignificant amount
96 // (0.25ulp) from rat and check that it rounds to the same value:
97 //
98 typename exponent_type<Float>::type e;
99 Float t = frexp(val, &e);
100 (void)t; // warning suppression
101 e -= std::numeric_limits<Float>::digits + 2;
102 BOOST_ASSERT(val == (val + ldexp(Float(1), e)));
103 Rat delta, rounded;
104 typedef typename boost::multiprecision::component_type<Rat>::type i_type;
105 i_type i(1);
106 i <<= (e < 0 ? -e : e);
107 if (e > 0)
108 delta.assign(i);
109 else
110 delta = Rat(i_type(1), i);
111 rounded = rat + delta;
112 new_f = static_cast<Float>(rounded);
113 BOOST_CHECK_EQUAL(val, new_f);
114 rounded = rat - delta;
115 new_f = static_cast<Float>(rounded);
116 BOOST_CHECK_EQUAL(val, new_f);
117
118 delta /= 2;
119 rounded = rat + delta;
120 new_f = static_cast<Float>(rounded);
121 BOOST_CHECK_EQUAL(val, new_f);
122 rounded = rat - delta;
123 new_f = static_cast<Float>(rounded);
124 BOOST_CHECK_EQUAL(val, new_f);
125
126 delta /= 2;
127 rounded = rat + delta;
128 new_f = static_cast<Float>(rounded);
129 BOOST_CHECK_EQUAL(val, new_f);
130 rounded = rat - delta;
131 new_f = static_cast<Float>(rounded);
132 BOOST_CHECK_EQUAL(val, new_f);
133 #endif
134 }
135
136 template <class Float, class Rat>
test_round_trip()137 void test_round_trip()
138 {
139 #ifndef BOOST_MP_NO_CXX11_EXPLICIT_CONVERSION_OPERATORS
140 std::cout << "Testing types " << typeid(Float).name() << " <<==>> " << typeid(Rat).name() << std::endl;
141 std::cout << "digits = " << std::numeric_limits<Float>::digits << std::endl;
142 std::cout << "digits10 = " << std::numeric_limits<Float>::digits10 << std::endl;
143 #ifndef BOOST_NO_CXX11_NUMERIC_LIMITS
144 std::cout << "max_digits10 = " << std::numeric_limits<Float>::max_digits10 << std::endl;
145 #endif
146
147 stopwatch<boost::chrono::high_resolution_clock> w;
148
149 int count = 0;
150
151 #ifndef CI_SUPPRESS_KNOWN_ISSUES
152 while (boost::chrono::duration_cast<boost::chrono::duration<double> >(w.elapsed()).count() < 200)
153 #else
154 while (boost::chrono::duration_cast<boost::chrono::duration<double> >(w.elapsed()).count() < 50)
155 #endif
156 {
157 Float val = generate_random_float<Float>();
158 do_round_trip<Float, Rat>(val);
159 do_round_trip<Float, Rat>(Float(-val));
160 do_round_trip<Float, Rat>(Float(1 / val));
161 do_round_trip<Float, Rat>(Float(-1 / val));
162 count += 4;
163 if (boost::detail::test_errors() > 100)
164 break;
165 }
166
167 std::cout << "Execution time = " << boost::chrono::duration_cast<boost::chrono::duration<double> >(w.elapsed()).count() << "s" << std::endl;
168 std::cout << "Total values tested: " << count << std::endl;
169 #endif
170 }
171
172 template <class Int>
generate_random_int()173 Int generate_random_int()
174 {
175 static boost::random::mt19937 gen;
176 static boost::random::uniform_int_distribution<boost::random::mt19937::result_type> d(1, 20);
177
178 int lim;
179 Int cppi(0);
180
181 lim = d(gen);
182
183 for (int i = 0; i < lim; ++i)
184 {
185 cppi *= (gen.max)();
186 cppi += gen();
187 }
188 return cppi;
189 }
190
191 template <class Float, class Rat>
test_random_rationals()192 void test_random_rationals()
193 {
194 #ifndef BOOST_MP_NO_CXX11_EXPLICIT_CONVERSION_OPERATORS
195 std::cout << "Testing types " << typeid(Float).name() << " <<==>> " << typeid(Rat).name() << std::endl;
196 std::cout << "digits = " << std::numeric_limits<Float>::digits << std::endl;
197 std::cout << "digits10 = " << std::numeric_limits<Float>::digits10 << std::endl;
198 #ifndef BOOST_NO_CXX11_NUMERIC_LIMITS
199 std::cout << "max_digits10 = " << std::numeric_limits<Float>::max_digits10 << std::endl;
200 #endif
201
202 typedef typename boost::multiprecision::component_type<Rat>::type i_type;
203 stopwatch<boost::chrono::high_resolution_clock> w;
204
205 int count = 0;
206
207 #ifndef CI_SUPPRESS_KNOWN_ISSUES
208 while (boost::chrono::duration_cast<boost::chrono::duration<double> >(w.elapsed()).count() < 200)
209 #else
210 while (boost::chrono::duration_cast<boost::chrono::duration<double> >(w.elapsed()).count() < 50)
211 #endif
212 {
213 Rat rat(generate_random_int<i_type>(), generate_random_int<i_type>());
214 Float f(rat);
215 Rat new_rat(f); // rounded value
216 int c = new_rat.compare(rat);
217 if (c < 0)
218 {
219 // If f was rounded down, next float up must be above the original value:
220 f = boost::math::float_next(f);
221 new_rat.assign(f);
222 BOOST_CHECK(new_rat >= rat);
223 }
224 else if (c > 0)
225 {
226 // If f was rounded up, next float down must be below the original value:
227 f = boost::math::float_prior(f);
228 new_rat.assign(f);
229 BOOST_CHECK(new_rat <= rat);
230 }
231 else
232 {
233 // Values were equal... nothing to test.
234 }
235 if (boost::detail::test_errors() > 100)
236 break;
237 }
238
239 std::cout << "Execution time = " << boost::chrono::duration_cast<boost::chrono::duration<double> >(w.elapsed()).count() << "s" << std::endl;
240 std::cout << "Total values tested: " << count << std::endl;
241 #endif
242 }
243
244 #if defined(TEST2)
245
double_spot_tests()246 void double_spot_tests()
247 {
248 boost::multiprecision::cpp_rational rat = 1;
249 boost::multiprecision::cpp_rational twiddle(boost::multiprecision::cpp_int(1), boost::multiprecision::cpp_int(boost::multiprecision::cpp_int(1) << 54));
250 rat += boost::multiprecision::cpp_rational(boost::multiprecision::cpp_int(1), boost::multiprecision::cpp_int(boost::multiprecision::cpp_int(1) << 50));
251
252 double d = rat.convert_to<double>();
253
254 rat += twiddle;
255 BOOST_CHECK_EQUAL(d, rat.convert_to<double>());
256 rat += twiddle;
257 // tie: round to even rounds down
258 BOOST_CHECK_EQUAL(d, rat.convert_to<double>());
259 rat += twiddle;
260 BOOST_CHECK_NE(d, rat.convert_to<double>());
261 rat -= twiddle;
262 BOOST_CHECK_EQUAL(d, rat.convert_to<double>());
263 rat += boost::multiprecision::cpp_rational(boost::multiprecision::cpp_int(1), boost::multiprecision::cpp_int(boost::multiprecision::cpp_int(1) << 52));
264 // tie, but last bit is now a 1 so we round up:
265 BOOST_CHECK_NE(d, rat.convert_to<double>());
266 }
267
268 #endif
269
main()270 int main()
271 {
272 using namespace boost::multiprecision;
273 #if defined(TEST1) && !defined(BOOST_MSVC)
274 test_round_trip<number<cpp_bin_float<113, digit_base_2, void, boost::int16_t> >, cpp_rational>();
275 #elif defined(TEST2)
276 double_spot_tests();
277 test_round_trip<double, cpp_rational>();
278 #elif defined(TEST3) && !defined(BOOST_MSVC)
279 test_random_rationals<number<cpp_bin_float<113, digit_base_2, void, boost::int16_t> >, cpp_rational>();
280 #elif defined(TEST4)
281 test_random_rationals<double, cpp_rational>();
282 #elif defined(TEST5)
283 // This does not work: gmp does not correctly round integer to float or
284 // rational to float conversions:
285 test_round_trip<double, mpq_rational>();
286 #elif defined(TEST6)
287 test_round_trip<mpfr_float_100, mpq_rational>();
288 #elif defined(TEST7)
289 test_random_rationals<mpfr_float_100, mpq_rational>();
290 #elif defined(TEST8)
291 test_random_rationals<double, mpq_rational>();
292 #endif
293 return boost::report_errors();
294 }
295