libs/multiprecision/test/test_mpfi.cpp

   1 ///////////////////////////////////////////////////////////////
   2 //  Copyright 2012 John Maddock. Distributed under the Boost
   3 //  Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
   5
   6 #ifdef _MSC_VER
   7 #define _SCL_SECURE_NO_WARNINGS
   8 #endif
   9
  10 #include "test.hpp"
  11 #include <boost/multiprecision/mpfi.hpp>
  12 #include <boost/multiprecision/random.hpp>
  13
  14 using namespace boost::multiprecision;
  15 using namespace boost::random;
  16
  17 void test_exp()
  18 {
  19    std::cout << "Testing exp\n";
  20
  21    mpfr_float_50 val = 1.25;
  22
  23    for (unsigned i = 0; i < 2000; ++i)
  24    {
  25       mpfr_float_100 a(val);
  26       mpfr_float_100 b  = exp(a);
  27       mpfi_float_50  in = val;
  28       in                = exp(in);
  29       BOOST_CHECK((boost::math::isfinite)(in));
  30       BOOST_CHECK(lower(in) <= b);
  31       BOOST_CHECK(upper(in) >= b);
  32       b  = log(a);
  33       in = val;
  34       in = log(in);
  35       BOOST_CHECK((boost::math::isfinite)(in));
  36       BOOST_CHECK(lower(in) <= b);
  37       BOOST_CHECK(upper(in) >= b);
  38       val *= 1.01;
  39    }
  40    val = 1;
  41    for (unsigned i = 0; i < 2000; ++i)
  42    {
  43       mpfr_float_100 a(val);
  44       mpfr_float_100 b  = exp(a);
  45       mpfi_float_50  in = val;
  46       in                = exp(in);
  47       BOOST_CHECK((boost::math::isfinite)(in));
  48       BOOST_CHECK(lower(in) <= b);
  49       BOOST_CHECK(upper(in) >= b);
  50       b  = log(a);
  51       in = val;
  52       in = log(in);
  53       BOOST_CHECK((boost::math::isfinite)(in));
  54       BOOST_CHECK(lower(in) <= b);
  55       BOOST_CHECK(upper(in) >= b);
  56       val /= 1.01;
  57    }
  58 }
  59
  60 void test_pow()
  61 {
  62    std::cout << "Testing pow function\n";
  63
  64    mt19937                                  gen;
  65    uniform_real_distribution<mpfr_float_50> dist1(0, 400);
  66
  67    for (unsigned i = 0; i < 5000; ++i)
  68    {
  69       mpfr_float_50 base, p;
  70       base = dist1(gen);
  71       p    = dist1(gen);
  72       mpfr_float_100 a, b, r;
  73       a = base;
  74       b = p;
  75       r = pow(a, b);
  76       mpfi_float_50 ai, bi, ri;
  77       ai = base;
  78       bi = p;
  79       ri = pow(ai, bi);
  80       BOOST_CHECK((boost::math::isfinite)(ri));
  81       BOOST_CHECK(lower(ri) <= r);
  82       BOOST_CHECK(upper(ri) >= r);
  83    }
  84 }
  85
  86 void test_trig()
  87 {
  88    std::cout << "Testing trig functions\n";
  89
  90    mt19937                                  gen;
  91    uniform_real_distribution<mpfr_float_50> dist1(-1.57079632679, 1.57079632679);
  92    uniform_real_distribution<mpfr_float_50> dist2(-1, 1);
  93
  94    for (unsigned i = 0; i < 5000; ++i)
  95    {
  96       mpfr_float_50 val;
  97       val               = dist1(gen);
  98       mpfr_float_100 a  = val;
  99       mpfr_float_100 b  = sin(a);
 100       mpfi_float_50  a2 = val;
 101       mpfi_float_50  b2 = sin(a2);
 102       BOOST_CHECK((boost::math::isfinite)(b2));
 103       BOOST_CHECK(lower(b2) <= b);
 104       BOOST_CHECK(upper(b2) >= b);
 105       b  = cos(a);
 106       b2 = cos(a2);
 107       BOOST_CHECK((boost::math::isfinite)(b2));
 108       BOOST_CHECK(lower(b2) <= b);
 109       BOOST_CHECK(upper(b2) >= b);
 110       b  = tan(a);
 111       b2 = tan(a2);
 112       BOOST_CHECK((boost::math::isfinite)(b2));
 113       BOOST_CHECK(lower(b2) <= b);
 114       BOOST_CHECK(upper(b2) >= b);
 115    }
 116    for (unsigned i = 0; i < 5000; ++i)
 117    {
 118       mpfr_float_50 val;
 119       val               = dist2(gen);
 120       mpfr_float_100 a  = val;
 121       mpfr_float_100 b  = asin(a);
 122       mpfi_float_50  a2 = val;
 123       mpfi_float_50  b2 = asin(a2);
 124       BOOST_CHECK((boost::math::isfinite)(b2));
 125       BOOST_CHECK(lower(b2) <= b);
 126       BOOST_CHECK(upper(b2) >= b);
 127       b  = acos(a);
 128       b2 = acos(a2);
 129       BOOST_CHECK((boost::math::isfinite)(b2));
 130       BOOST_CHECK(lower(b2) <= b);
 131       BOOST_CHECK(upper(b2) >= b);
 132       b  = atan(a);
 133       b2 = atan(a2);
 134       BOOST_CHECK((boost::math::isfinite)(b2));
 135       BOOST_CHECK(lower(b2) <= b);
 136       BOOST_CHECK(upper(b2) >= b);
 137    }
 138 }
 139
 140 void test_hyp()
 141 {
 142    std::cout << "Testing hyperbolic trig functions\n";
 143
 144    mt19937                                  gen;
 145    uniform_real_distribution<mpfr_float_50> dist1(-10, 10);
 146    uniform_real_distribution<mpfr_float_50> dist2(-1, 1);
 147
 148    for (unsigned i = 0; i < 5000; ++i)
 149    {
 150       mpfr_float_50 val;
 151       val               = dist1(gen);
 152       mpfr_float_100 a  = val;
 153       mpfr_float_100 b  = sinh(a);
 154       mpfi_float_50  a2 = val;
 155       mpfi_float_50  b2 = sinh(a2);
 156       BOOST_CHECK((boost::math::isfinite)(b2));
 157       BOOST_CHECK(lower(b2) <= b);
 158       BOOST_CHECK(upper(b2) >= b);
 159       b  = cosh(a);
 160       b2 = cosh(a2);
 161       BOOST_CHECK((boost::math::isfinite)(b2));
 162       BOOST_CHECK(lower(b2) <= b);
 163       BOOST_CHECK(upper(b2) >= b);
 164       b  = tanh(a);
 165       b2 = tanh(a2);
 166       BOOST_CHECK((boost::math::isfinite)(b2));
 167       BOOST_CHECK(lower(b2) <= b);
 168       BOOST_CHECK(upper(b2) >= b);
 169    }
 170 }
 171
 172 void test_intervals()
 173 {
 174    mpfi_float_50 a(1, 2);
 175    mpfi_float_50 b(1.5, 2.5);
 176    BOOST_CHECK_EQUAL(lower(a), 1);
 177    BOOST_CHECK_EQUAL(upper(a), 2);
 178    BOOST_CHECK_EQUAL(median(a), 1.5);
 179    BOOST_CHECK_EQUAL(width(a), 1);
 180    mpfi_float_50 r = intersect(a, b);
 181    BOOST_CHECK_EQUAL(lower(r), 1.5);
 182    BOOST_CHECK_EQUAL(upper(r), 2);
 183    r = hull(a, b);
 184    BOOST_CHECK_EQUAL(lower(r), 1);
 185    BOOST_CHECK_EQUAL(upper(r), 2.5);
 186    BOOST_CHECK(overlap(a, b));
 187    BOOST_CHECK(in(mpfr_float_50(1.5), a));
 188    BOOST_CHECK(in(mpfr_float_50(1), a));
 189    BOOST_CHECK(in(mpfr_float_50(2), a));
 190    BOOST_CHECK(!zero_in(a));
 191    b = mpfi_float_50(1.5, 1.75);
 192    BOOST_CHECK(subset(b, a));
 193    BOOST_CHECK(proper_subset(b, a));
 194    BOOST_CHECK(!empty(a));
 195    BOOST_CHECK(!singleton(a));
 196    b = mpfi_float_50(5, 6);
 197    r = intersect(a, b);
 198    BOOST_CHECK(empty(r));
 199 }
 200
 201 #ifdef TEST_SPECIAL
 202 #include "math/table_type.hpp"
 203 #include <boost/math/special_functions.hpp>
 204
 205 #define T mpfi_float_50
 206
 207 typedef number<mpfi_float_backend<25> > mpfi_float_25;
 208
 209 void test_log1p_expm1()
 210 {
 211 #include "../../math/test/log1p_expm1_data.ipp"
 212
 213    std::cout << std::setprecision(std::numeric_limits<mpfi_float_25>::max_digits10);
 214    std::cout << "Testing log1p and expm1\n";
 215
 216    for (unsigned i = 0; i < log1p_expm1_data.size(); ++i)
 217    {
 218       mpfi_float_25 in(log1p_expm1_data[i][0]);
 219       mpfi_float_25 out = boost::math::log1p(in);
 220       mpfi_float_25 expected(log1p_expm1_data[i][1]);
 221       if (!subset(expected, out))
 222       {
 223          std::cout << in << std::endl;
 224          std::cout << out << std::endl;
 225          std::cout << expected << std::endl;
 226          BOOST_CHECK(lower(out) <= lower(expected));
 227          BOOST_CHECK(upper(out) >= upper(expected));
 228       }
 229       out      = boost::math::expm1(in);
 230       expected = mpfi_float_25(log1p_expm1_data[i][2]);
 231       if (!subset(expected, out))
 232       {
 233          std::cout << in << std::endl;
 234          std::cout << out << std::endl;
 235          std::cout << expected << std::endl;
 236          BOOST_CHECK(lower(out) <= lower(expected));
 237          BOOST_CHECK(upper(out) >= upper(expected));
 238       }
 239    }
 240 }
 241
 242 void test_bessel()
 243 {
 244 #include "../../math/test/bessel_i_int_data.ipp"
 245 #include "../../math/test/bessel_i_data.ipp"
 246
 247    std::cout << std::setprecision(std::numeric_limits<mpfi_float_25>::max_digits10);
 248    std::cout << "Testing Bessel Functions\n";
 249
 250    for (unsigned i = 0; i < bessel_i_int_data.size(); ++i)
 251    {
 252       int           v = boost::lexical_cast<int>(static_cast<const char*>(bessel_i_int_data[i][0]));
 253       mpfi_float_25 in(bessel_i_int_data[i][1]);
 254       mpfi_float_25 out = boost::math::cyl_bessel_i(v, in);
 255       mpfi_float_25 expected(bessel_i_int_data[i][2]);
 256       if (!subset(expected, out))
 257       {
 258          std::cout << in << std::endl;
 259          std::cout << out << std::endl;
 260          std::cout << expected << std::endl;
 261          BOOST_CHECK(lower(out) <= lower(expected));
 262          BOOST_CHECK(upper(out) >= upper(expected));
 263       }
 264    }
 265 }
 266
 267 #endif
 268
 269 int main()
 270 {
 271 #ifdef TEST_SPECIAL
 272    test_log1p_expm1();
 273    test_bessel();
 274 #endif
 275    test_intervals();
 276    test_exp();
 277    test_pow();
 278    test_trig();
 279    test_hyp();
 280    return boost::report_errors();
 281 }