test_next.cpp 10 KB

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