libs/math/test/test_next.cpp

   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>
  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>
  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 mulyiple 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
 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