1 // (C) Copyright John Maddock 2008.
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 #include <pch.hpp>
7
8 #include <boost/math/concepts/real_concept.hpp>
9 #include <boost/math/tools/test.hpp>
10 #define BOOST_TEST_MAIN
11 #include <boost/test/unit_test.hpp>
12 #include <boost/test/tools/floating_point_comparison.hpp>
13 #include <boost/math/special_functions/next.hpp>
14 #include <boost/math/special_functions/ulp.hpp>
15 #include <boost/multiprecision/cpp_bin_float.hpp>
16 #include <iostream>
17 #include <iomanip>
18
19 #ifdef BOOST_MSVC
20 #pragma warning(disable:4127)
21 #endif
22
23 #if !defined(_CRAYC) && !defined(__CUDACC__) && (!defined(__GNUC__) || (__GNUC__ > 3) || ((__GNUC__ == 3) && (__GNUC_MINOR__ > 3)))
24 #if (defined(_M_IX86_FP) && (_M_IX86_FP >= 2)) || defined(__SSE2__) || defined(TEST_SSE2)
25 #include <float.h>
26 #include "xmmintrin.h"
27 #define TEST_SSE2
28 #endif
29 #endif
30
31
32 template <class T>
test_value(const T & val,const char * name)33 void test_value(const T& val, const char* name)
34 {
35 using namespace boost::math;
36 T upper = tools::max_value<T>();
37 T lower = -upper;
38
39 std::cout << "Testing type " << name << " with initial value " << val << std::endl;
40
41 BOOST_CHECK_EQUAL(float_distance(float_next(val), val), -1);
42 BOOST_CHECK(float_next(val) > val);
43 BOOST_CHECK_EQUAL(float_distance(float_prior(val), val), 1);
44 BOOST_CHECK(float_prior(val) < val);
45 BOOST_CHECK_EQUAL(float_distance((boost::math::nextafter)(val, upper), val), -1);
46 BOOST_CHECK((boost::math::nextafter)(val, upper) > val);
47 BOOST_CHECK_EQUAL(float_distance((boost::math::nextafter)(val, lower), val), 1);
48 BOOST_CHECK((boost::math::nextafter)(val, lower) < val);
49 BOOST_CHECK_EQUAL(float_distance(float_next(float_next(val)), val), -2);
50 BOOST_CHECK_EQUAL(float_distance(float_prior(float_prior(val)), val), 2);
51 BOOST_CHECK_EQUAL(float_distance(float_prior(float_prior(val)), float_next(float_next(val))), 4);
52 BOOST_CHECK_EQUAL(float_distance(float_prior(float_next(val)), val), 0);
53 BOOST_CHECK_EQUAL(float_distance(float_next(float_prior(val)), val), 0);
54 BOOST_CHECK_EQUAL(float_prior(float_next(val)), val);
55 BOOST_CHECK_EQUAL(float_next(float_prior(val)), val);
56
57 BOOST_CHECK_EQUAL(float_distance(float_advance(val, 4), val), -4);
58 BOOST_CHECK_EQUAL(float_distance(float_advance(val, -4), val), 4);
59 if(std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == std::denorm_present))
60 {
61 BOOST_CHECK_EQUAL(float_distance(float_advance(float_next(float_next(val)), 4), float_next(float_next(val))), -4);
62 BOOST_CHECK_EQUAL(float_distance(float_advance(float_next(float_next(val)), -4), float_next(float_next(val))), 4);
63 }
64 if(val > 0)
65 {
66 T n = val + ulp(val);
67 T fn = float_next(val);
68 if(n > fn)
69 {
70 BOOST_CHECK_LE(ulp(val), boost::math::tools::min_value<T>());
71 }
72 else
73 {
74 BOOST_CHECK_EQUAL(fn, n);
75 }
76 }
77 else if(val == 0)
78 {
79 BOOST_CHECK_GE(boost::math::tools::min_value<T>(), ulp(val));
80 }
81 else
82 {
83 T n = val - ulp(val);
84 T fp = float_prior(val);
85 if(n < fp)
86 {
87 BOOST_CHECK_LE(ulp(val), boost::math::tools::min_value<T>());
88 }
89 else
90 {
91 BOOST_CHECK_EQUAL(fp, n);
92 }
93 }
94 }
95
96 template <class T>
test_values(const T & val,const char * name)97 void test_values(const T& val, const char* name)
98 {
99 static const T a = static_cast<T>(1.3456724e22);
100 static const T b = static_cast<T>(1.3456724e-22);
101 static const T z = 0;
102 static const T one = 1;
103 static const T two = 2;
104
105 std::cout << "Testing type " << name << std::endl;
106
107 T den = (std::numeric_limits<T>::min)() / 4;
108 if(den != 0)
109 {
110 std::cout << "Denormals are active\n";
111 }
112 else
113 {
114 std::cout << "Denormals are flushed to zero.\n";
115 }
116
117 test_value(a, name);
118 test_value(-a, name);
119 test_value(b, name);
120 test_value(-b, name);
121 test_value(boost::math::tools::epsilon<T>(), name);
122 test_value(-boost::math::tools::epsilon<T>(), name);
123 test_value(boost::math::tools::min_value<T>(), name);
124 test_value(-boost::math::tools::min_value<T>(), name);
125 if (std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == std::denorm_present) && ((std::numeric_limits<T>::min)() / 2 != 0))
126 {
127 test_value(z, name);
128 test_value(-z, name);
129 }
130 test_value(one, name);
131 test_value(-one, name);
132 test_value(two, name);
133 test_value(-two, name);
134 #if defined(TEST_SSE2)
135 if((_mm_getcsr() & (_MM_FLUSH_ZERO_ON | 0x40)) == 0)
136 {
137 #endif
138 if(std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == std::denorm_present) && ((std::numeric_limits<T>::min)() / 2 != 0))
139 {
140 test_value(std::numeric_limits<T>::denorm_min(), name);
141 test_value(-std::numeric_limits<T>::denorm_min(), name);
142 test_value(2 * std::numeric_limits<T>::denorm_min(), name);
143 test_value(-2 * std::numeric_limits<T>::denorm_min(), name);
144 }
145 #if defined(TEST_SSE2)
146 }
147 #endif
148 static const int primes[] = {
149 11, 13, 17, 19, 23, 29,
150 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
151 73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
152 127, 131, 137, 139, 149, 151, 157, 163, 167, 173,
153 179, 181, 191, 193, 197, 199, 211, 223, 227, 229,
154 233, 239, 241, 251, 257, 263, 269, 271, 277, 281,
155 283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
156 353, 359, 367, 373, 379, 383, 389, 397, 401, 409,
157 419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
158 };
159
160 for(unsigned i = 0; i < sizeof(primes)/sizeof(primes[0]); ++i)
161 {
162 T v1 = val;
163 T v2 = val;
164 for(int j = 0; j < primes[i]; ++j)
165 {
166 v1 = boost::math::float_next(v1);
167 v2 = boost::math::float_prior(v2);
168 }
169 BOOST_CHECK_EQUAL(boost::math::float_distance(v1, val), -primes[i]);
170 BOOST_CHECK_EQUAL(boost::math::float_distance(v2, val), primes[i]);
171 BOOST_CHECK_EQUAL(boost::math::float_advance(val, primes[i]), v1);
172 BOOST_CHECK_EQUAL(boost::math::float_advance(val, -primes[i]), v2);
173 }
174 if(std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_infinity))
175 {
176 BOOST_CHECK_EQUAL(boost::math::float_prior(std::numeric_limits<T>::infinity()), (std::numeric_limits<T>::max)());
177 BOOST_CHECK_EQUAL(boost::math::float_next(-std::numeric_limits<T>::infinity()), -(std::numeric_limits<T>::max)());
178 BOOST_MATH_CHECK_THROW(boost::math::float_prior(-std::numeric_limits<T>::infinity()), std::domain_error);
179 BOOST_MATH_CHECK_THROW(boost::math::float_next(std::numeric_limits<T>::infinity()), std::domain_error);
180 if(boost::math::policies:: BOOST_MATH_OVERFLOW_ERROR_POLICY == boost::math::policies::throw_on_error)
181 {
182 BOOST_MATH_CHECK_THROW(boost::math::float_prior(-(std::numeric_limits<T>::max)()), std::overflow_error);
183 BOOST_MATH_CHECK_THROW(boost::math::float_next((std::numeric_limits<T>::max)()), std::overflow_error);
184 }
185 else
186 {
187 BOOST_CHECK_EQUAL(boost::math::float_prior(-(std::numeric_limits<T>::max)()), -std::numeric_limits<T>::infinity());
188 BOOST_CHECK_EQUAL(boost::math::float_next((std::numeric_limits<T>::max)()), std::numeric_limits<T>::infinity());
189 }
190 }
191 //
192 // We need to test float_distance over multiple orders of magnitude,
193 // the only way to get an accurate true result is to count the representations
194 // between the two end points, but we can only really do this for type float:
195 //
196 if (std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::digits < 30) && (std::numeric_limits<T>::radix == 2))
197 {
198 T left, right, dist, fresult;
199 boost::uintmax_t result;
200
201 left = static_cast<T>(0.1);
202 right = left * static_cast<T>(4.2);
203 dist = boost::math::float_distance(left, right);
204 // We have to use a wider integer type for the accurate count, since there
205 // aren't enough bits in T to get a true result if the values differ
206 // by more than a factor of 2:
207 result = 0;
208 for (; left != right; ++result, left = boost::math::float_next(left));
209 fresult = static_cast<T>(result);
210 BOOST_CHECK_EQUAL(fresult, dist);
211
212 left = static_cast<T>(-0.1);
213 right = left * static_cast<T>(4.2);
214 dist = boost::math::float_distance(right, left);
215 result = 0;
216 for (; left != right; ++result, left = boost::math::float_prior(left));
217 fresult = static_cast<T>(result);
218 BOOST_CHECK_EQUAL(fresult, dist);
219
220 left = static_cast<T>(-1.1) * (std::numeric_limits<T>::min)();
221 right = static_cast<T>(-4.1) * left;
222 dist = boost::math::float_distance(left, right);
223 result = 0;
224 for (; left != right; ++result, left = boost::math::float_next(left));
225 fresult = static_cast<T>(result);
226 BOOST_CHECK_EQUAL(fresult, dist);
227 }
228 }
229
BOOST_AUTO_TEST_CASE(test_main)230 BOOST_AUTO_TEST_CASE( test_main )
231 {
232 test_values(1.0f, "float");
233 test_values(1.0, "double");
234 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
235 test_values(1.0L, "long double");
236 test_values(boost::math::concepts::real_concept(0), "real_concept");
237 #endif
238
239 //
240 // Test some multiprecision types:
241 //
242 test_values(boost::multiprecision::cpp_bin_float_quad(0), "cpp_bin_float_quad");
243 // This is way to slow to test routinely:
244 //test_values(boost::multiprecision::cpp_bin_float_single(0), "cpp_bin_float_single");
245 test_values(boost::multiprecision::cpp_bin_float_50(0), "cpp_bin_float_50");
246
247 #if defined(TEST_SSE2)
248
249 int mmx_flags = _mm_getcsr(); // We'll restore these later.
250
251 #ifdef _WIN32
252 // These tests fail pretty badly on Linux x64, especially with Intel-12.1
253 _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
254 std::cout << "Testing again with Flush-To-Zero set" << std::endl;
255 std::cout << "SSE2 control word is: " << std::hex << _mm_getcsr() << std::endl;
256 test_values(1.0f, "float");
257 test_values(1.0, "double");
258 _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_OFF);
259 #endif
260 BOOST_ASSERT((_mm_getcsr() & 0x40) == 0);
261 _mm_setcsr(_mm_getcsr() | 0x40);
262 std::cout << "Testing again with Denormals-Are-Zero set" << std::endl;
263 std::cout << "SSE2 control word is: " << std::hex << _mm_getcsr() << std::endl;
264 test_values(1.0f, "float");
265 test_values(1.0, "double");
266
267 // Restore the MMX flags:
268 _mm_setcsr(mmx_flags);
269 #endif
270
271 }
272
273
274