libs/multiprecision/test/test_rat_float_interconv.cpp

   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;
  36    stopwatch()
  37    {
  38       m_start = Clock::now();
  39    }
  40    duration elapsed()
  41    {
  42       return Clock::now() - m_start;
  43    }
  44    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>
  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>
  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>
 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>
 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>
 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
 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
 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 }