libs/math/test/complex_test.cpp

   1 //  (C) Copyright John Maddock 2005.
   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 #define BOOST_TEST_MAIN
   7 #include <boost/test/unit_test.hpp>
   8 #include <boost/test/tools/floating_point_comparison.hpp>
   9 #include <boost/type_traits/is_same.hpp>
  10 #include <boost/type_traits/is_floating_point.hpp>
  11 #include <boost/mpl/if.hpp>
  12 #include <boost/static_assert.hpp>
  13 #include <boost/math/complex.hpp>
  14
  15 #include <iostream>
  16 #include <iomanip>
  17 #include <cmath>
  18 #include <typeinfo>
  19
  20 #ifdef BOOST_NO_STDC_NAMESPACE
  21 namespace std{ using ::sqrt; using ::tan; using ::tanh; }
  22 #endif
  23
  24 #ifndef VERBOSE
  25 #undef BOOST_TEST_MESSAGE
  26 #define BOOST_TEST_MESSAGE(x)
  27 #endif
  28
  29 //
  30 // check_complex:
  31 // Verifies that expected value "a" and found value "b" have a relative error
  32 // less than "max_error" epsilons.  Note that relative error is calculated for
  33 // the complex number as a whole; this means that the error in the real or
  34 // imaginary parts alone can be much higher than max_error when the real and
  35 // imaginary parts are of very different magnitudes.  This is important, because
  36 // the Hull et al analysis of the acos and asin algorithms requires that very small
  37 // real/imaginary components can be safely ignored if they are negligible compared
  38 // to the other component.
  39 //
  40 template <class T>
  41 bool check_complex(const std::complex<T>& a, const std::complex<T>& b, int max_error)
  42 {
  43    //
  44    // a is the expected value, b is what was actually found,
  45    // compute | (a-b)/b | and compare with max_error which is the
  46    // multiple of E to permit:
  47    //
  48    bool result = true;
  49    static const std::complex<T> zero(0);
  50    static const T eps = std::pow(static_cast<T>(std::numeric_limits<T>::radix), static_cast<T>(1 - std::numeric_limits<T>::digits));
  51    if(a == zero)
  52    {
  53       if(b != zero)
  54       {
  55          if(boost::math::fabs(b) > eps)
  56          {
  57             result = false;
  58             BOOST_ERROR("Expected {0,0} but got: " << b);
  59          }
  60          else
  61          {
  62             BOOST_TEST_MESSAGE("Expected {0,0} but got: " << b);
  63          }
  64       }
  65       return result;
  66    }
  67    else if(b == zero)
  68    {
  69       if(boost::math::fabs(a) > eps)
  70       {
  71          BOOST_ERROR("Found {0,0} but expected: " << a);
  72          return false;;
  73       }
  74       else
  75       {
  76          BOOST_TEST_MESSAGE("Found {0,0} but expected: " << a);
  77       }
  78    }
  79
  80    if((boost::math::isnan)(a.real()))
  81    {
  82       BOOST_ERROR("Found non-finite value for real part: " << a);
  83    }
  84    if((boost::math::isnan)(a.imag()))
  85    {
  86       BOOST_ERROR("Found non-finite value for inaginary part: " << a);
  87    }
  88
  89    T rel = boost::math::fabs((b-a)/b) / eps;
  90    if( rel > max_error)
  91    {
  92       result = false;
  93       BOOST_ERROR("Error in result exceeded permitted limit of " << max_error << " (actual relative error was " << rel << "e).  Found " << b << " expected " << a);
  94    }
  95    return result;
  96 }
  97
  98 //
  99 // test_inverse_trig:
 100 // This is nothing more than a sanity check, computes trig(atrig(z))
 101 // and compare the result to z.  Note that:
 102 //
 103 // atrig(trig(z)) != z
 104 //
 105 // for certain z because the inverse trig functions are multi-valued, this
 106 // essentially rules this out as a testing method.  On the other hand:
 107 //
 108 // trig(atrig(z))
 109 //
 110 // can vary compare to z by an arbitrarily large amount.  For one thing we
 111 // have no control over the implementation of the trig functions, for another
 112 // even if both functions were accurate to 1ulp (as accurate as transcendental
 113 // number can get, thanks to the "table makers dilemma"), the errors can still
 114 // be arbitrarily large - often the inverse trig functions will map a very large
 115 // part of the complex domain into a small output domain, so you can never get
 116 // back exactly where you started from.  Consequently these tests are no more than
 117 // sanity checks (just verifies that signs are correct and so on).
 118 //
 119 template <class T>
 120 void test_inverse_trig(T)
 121 {
 122    using namespace std;
 123
 124    static const T interval = static_cast<T>(2.0L/128.0L);
 125
 126    T x, y;
 127
 128    std::cout << std::setprecision(std::numeric_limits<T>::digits10+2);
 129
 130    for(x = -1; x <= 1; x += interval)
 131    {
 132       for(y = -1; y <= 1; y += interval)
 133       {
 134          // acos:
 135          std::complex<T> val(x, y), inter, result;
 136          inter = boost::math::acos(val);
 137          result = cos(inter);
 138          if(!check_complex(val, result, 50))
 139          {
 140             std::cout << "Error in testing inverse complex cos for type " << typeid(T).name() << std::endl;
 141             std::cout << "   val=             " << val << std::endl;
 142             std::cout << "   acos(val) =      " << inter << std::endl;
 143             std::cout << "   cos(acos(val)) = " << result << std::endl;
 144          }
 145          // asin:
 146          inter = boost::math::asin(val);
 147          result = sin(inter);
 148          if(!check_complex(val, result, 5))
 149          {
 150             std::cout << "Error in testing inverse complex sin for type " << typeid(T).name() << std::endl;
 151             std::cout << "   val=             " << val << std::endl;
 152             std::cout << "   asin(val) =      " << inter << std::endl;
 153             std::cout << "   sin(asin(val)) = " << result << std::endl;
 154          }
 155       }
 156    }
 157
 158    static const T interval2 = static_cast<T>(3.0L/256.0L);
 159    for(x = -3; x <= 3; x += interval2)
 160    {
 161       for(y = -3; y <= 3; y += interval2)
 162       {
 163          // asinh:
 164          std::complex<T> val(x, y), inter, result;
 165          inter = boost::math::asinh(val);
 166          result = sinh(inter);
 167          if(!check_complex(val, result, 5))
 168          {
 169             std::cout << "Error in testing inverse complex sinh for type " << typeid(T).name() << std::endl;
 170             std::cout << "   val=               " << val << std::endl;
 171             std::cout << "   asinh(val) =       " << inter << std::endl;
 172             std::cout << "   sinh(asinh(val)) = " << result << std::endl;
 173          }
 174          // acosh:
 175          if(!((y == 0) && (x <= 1))) // can't test along the branch cut
 176          {
 177             inter = boost::math::acosh(val);
 178             result = cosh(inter);
 179             if(!check_complex(val, result, 60))
 180             {
 181                std::cout << "Error in testing inverse complex cosh for type " << typeid(T).name() << std::endl;
 182                std::cout << "   val=               " << val << std::endl;
 183                std::cout << "   acosh(val) =       " << inter << std::endl;
 184                std::cout << "   cosh(acosh(val)) = " << result << std::endl;
 185             }
 186          }
 187          //
 188          // There is a problem in testing atan and atanh:
 189          // The inverse functions map a large input range to a much
 190          // smaller output range, so at the extremes too rather different
 191          // inputs may map to the same output value once rounded to N places.
 192          // Consequently tan(atan(z)) can suffer from arbitrarily large errors
 193          // even if individually they each have a small error bound.  On the other
 194          // hand we can't test atan(tan(z)) either because atan is multi-valued, so
 195          // round-tripping in this direction isn't always possible.
 196          // The following heuristic is designed to make the best of a bad job,
 197          // using atan(tan(z)) where possible and tan(atan(z)) when it's not.
 198          //
 199          static const int tanh_error = 20;
 200          if((0 != x) && (0 != y) && ((std::fabs(y) < 1) || (std::fabs(x) < 1)))
 201          {
 202             // atanh:
 203             val = boost::math::atanh(val);
 204             inter = tanh(val);
 205             result = boost::math::atanh(inter);
 206             if(!check_complex(val, result, tanh_error))
 207             {
 208                std::cout << "Error in testing inverse complex tanh for type " << typeid(T).name() << std::endl;
 209                std::cout << "   val=               " << val << std::endl;
 210                std::cout << "   tanh(val) =        " << inter << std::endl;
 211                std::cout << "   atanh(tanh(val)) = " << result << std::endl;
 212             }
 213             // atan:
 214             if(!((x == 0) && (std::fabs(y) == 1))) // we can't test infinities here
 215             {
 216                val = std::complex<T>(x, y);
 217                val = boost::math::atan(val);
 218                inter = tan(val);
 219                result = boost::math::atan(inter);
 220                if(!check_complex(val, result, tanh_error))
 221                {
 222                   std::cout << "Error in testing inverse complex tan for type " << typeid(T).name() << std::endl;
 223                   std::cout << "   val=               " << val << std::endl;
 224                   std::cout << "   tan(val) =         " << inter << std::endl;
 225                   std::cout << "   atan(tan(val)) =   " << result << std::endl;
 226                }
 227             }
 228          }
 229          else
 230          {
 231             // atanh:
 232             inter = boost::math::atanh(val);
 233             result = tanh(inter);
 234             if(!check_complex(val, result, tanh_error))
 235             {
 236                std::cout << "Error in testing inverse complex atanh for type " << typeid(T).name() << std::endl;
 237                std::cout << "   val=                 " << val << std::endl;
 238                std::cout << "   atanh(val) =         " << inter << std::endl;
 239                std::cout << "   tanh(atanh(val)) =   " << result << std::endl;
 240             }
 241             // atan:
 242             if(!((x == 0) && (std::fabs(y) == 1))) // we can't test infinities here
 243             {
 244                inter = boost::math::atan(val);
 245                result = tan(inter);
 246                if(!check_complex(val, result, tanh_error))
 247                {
 248                   std::cout << "Error in testing inverse complex atan for type " << typeid(T).name() << std::endl;
 249                   std::cout << "   val=                 " << val << std::endl;
 250                   std::cout << "   atan(val) =          " << inter << std::endl;
 251                   std::cout << "   tan(atan(val)) =     " << result << std::endl;
 252                }
 253             }
 254          }
 255       }
 256    }
 257 }
 258
 259 //
 260 // check_spots:
 261 // Various spot values, mostly the C99 special cases (infinites and NAN's).
 262 // TODO: add spot checks for the Wolfram spot values.
 263 //
 264 template <class T>
 265 void check_spots(const T&)
 266 {
 267    typedef std::complex<T> ct;
 268    ct result;
 269    static const T two = 2.0;
 270    T eps = std::pow(two, T(1-std::numeric_limits<T>::digits)); // numeric_limits<>::epsilon way too small to be useful on Darwin.
 271    static const T zero = 0;
 272    static const T mzero = -zero;
 273    static const T one = 1;
 274    static const T pi = boost::math::constants::pi<T>();
 275    static const T half_pi = boost::math::constants::half_pi<T>();
 276    static const T quarter_pi = half_pi / 2;
 277    static const T three_quarter_pi = boost::math::constants::three_quarters_pi<T>();
 278    T infinity = std::numeric_limits<T>::infinity();
 279    bool test_infinity = std::numeric_limits<T>::has_infinity;
 280    T nan = 0;
 281    bool test_nan = false;
 282 #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x564))
 283    // numeric_limits reports that a quiet NaN is present
 284    // but an attempt to access it will terminate the program!!!!
 285    if(std::numeric_limits<T>::has_quiet_NaN)
 286       nan = std::numeric_limits<T>::quiet_NaN();
 287    if((boost::math::isnan)(nan))
 288       test_nan = true;
 289 #endif
 290 #if defined(__DECCXX) && !defined(_IEEE_FP)
 291    // Tru64 cxx traps infinities unless the -ieee option is used:
 292    test_infinity = false;
 293 #endif
 294
 295    //
 296    // C99 spot tests for acos:
 297    //
 298    result = boost::math::acos(ct(zero));
 299    check_complex(ct(half_pi), result, 2);
 300
 301    result = boost::math::acos(ct(mzero));
 302    check_complex(ct(half_pi), result, 2);
 303
 304    result = boost::math::acos(ct(zero, mzero));
 305    check_complex(ct(half_pi), result, 2);
 306
 307    result = boost::math::acos(ct(mzero, mzero));
 308    check_complex(ct(half_pi), result, 2);
 309
 310    if(test_nan)
 311    {
 312       result = boost::math::acos(ct(zero,nan));
 313       BOOST_CHECK_CLOSE(result.real(), half_pi, eps*200);
 314       BOOST_CHECK((boost::math::isnan)(result.imag()));
 315
 316       result = boost::math::acos(ct(mzero,nan));
 317       BOOST_CHECK_CLOSE(result.real(), half_pi, eps*200);
 318       BOOST_CHECK((boost::math::isnan)(result.imag()));
 319    }
 320    if(test_infinity)
 321    {
 322       result = boost::math::acos(ct(zero, infinity));
 323       BOOST_CHECK_CLOSE(result.real(), half_pi, eps*200);
 324       BOOST_CHECK(result.imag() == -infinity);
 325
 326       result = boost::math::acos(ct(zero, -infinity));
 327       BOOST_CHECK_CLOSE(result.real(), half_pi, eps*200);
 328       BOOST_CHECK(result.imag() == infinity);
 329    }
 330
 331    if(test_nan)
 332    {
 333       result = boost::math::acos(ct(one, nan));
 334       BOOST_CHECK((boost::math::isnan)(result.real()));
 335       BOOST_CHECK((boost::math::isnan)(result.imag()));
 336    }
 337    if(test_infinity)
 338    {
 339       result = boost::math::acos(ct(-infinity, one));
 340       BOOST_CHECK_CLOSE(result.real(), pi, eps*200);
 341       BOOST_CHECK(result.imag() == -infinity);
 342
 343       result = boost::math::acos(ct(infinity, one));
 344       BOOST_CHECK(result.real() == 0);
 345       BOOST_CHECK(result.imag() == -infinity);
 346
 347       result = boost::math::acos(ct(-infinity, -one));
 348       BOOST_CHECK_CLOSE(result.real(), pi, eps*200);
 349       BOOST_CHECK(result.imag() == infinity);
 350
 351       result = boost::math::acos(ct(infinity, -one));
 352       BOOST_CHECK(result.real() == 0);
 353       BOOST_CHECK(result.imag() == infinity);
 354
 355       result = boost::math::acos(ct(-infinity, infinity));
 356       BOOST_CHECK_CLOSE(result.real(), three_quarter_pi, eps*200);
 357       BOOST_CHECK(result.imag() == -infinity);
 358
 359       result = boost::math::acos(ct(infinity, infinity));
 360       BOOST_CHECK_CLOSE(result.real(), quarter_pi, eps*200);
 361       BOOST_CHECK(result.imag() == -infinity);
 362
 363       result = boost::math::acos(ct(-infinity, -infinity));
 364       BOOST_CHECK_CLOSE(result.real(), three_quarter_pi, eps*200);
 365       BOOST_CHECK(result.imag() == infinity);
 366
 367       result = boost::math::acos(ct(infinity, -infinity));
 368       BOOST_CHECK_CLOSE(result.real(), quarter_pi, eps*200);
 369       BOOST_CHECK(result.imag() == infinity);
 370    }
 371    if(test_nan)
 372    {
 373       result = boost::math::acos(ct(infinity, nan));
 374       BOOST_CHECK((boost::math::isnan)(result.real()));
 375       BOOST_CHECK(std::fabs(result.imag()) == infinity);
 376
 377       result = boost::math::acos(ct(-infinity, nan));
 378       BOOST_CHECK((boost::math::isnan)(result.real()));
 379       BOOST_CHECK(std::fabs(result.imag()) == infinity);
 380
 381       result = boost::math::acos(ct(nan, zero));
 382       BOOST_CHECK((boost::math::isnan)(result.real()));
 383       BOOST_CHECK((boost::math::isnan)(result.imag()));
 384
 385       result = boost::math::acos(ct(nan, -zero));
 386       BOOST_CHECK((boost::math::isnan)(result.real()));
 387       BOOST_CHECK((boost::math::isnan)(result.imag()));
 388
 389       result = boost::math::acos(ct(nan, one));
 390       BOOST_CHECK((boost::math::isnan)(result.real()));
 391       BOOST_CHECK((boost::math::isnan)(result.imag()));
 392
 393       result = boost::math::acos(ct(nan, -one));
 394       BOOST_CHECK((boost::math::isnan)(result.real()));
 395       BOOST_CHECK((boost::math::isnan)(result.imag()));
 396
 397       result = boost::math::acos(ct(nan, nan));
 398       BOOST_CHECK((boost::math::isnan)(result.real()));
 399       BOOST_CHECK((boost::math::isnan)(result.imag()));
 400
 401       result = boost::math::acos(ct(nan, infinity));
 402       BOOST_CHECK((boost::math::isnan)(result.real()));
 403       BOOST_CHECK(result.imag() == -infinity);
 404
 405       result = boost::math::acos(ct(nan, -infinity));
 406       BOOST_CHECK((boost::math::isnan)(result.real()));
 407       BOOST_CHECK(result.imag() == infinity);
 408    }
 409    if(boost::math::signbit(mzero))
 410    {
 411       result = boost::math::acos(ct(-1.25f, zero));
 412       BOOST_CHECK(result.real() > 0);
 413       BOOST_CHECK(result.imag() < 0);
 414       result = boost::math::asin(ct(-1.75f, mzero));
 415       BOOST_CHECK(result.real() < 0);
 416       BOOST_CHECK(result.imag() < 0);
 417       result = boost::math::atan(ct(mzero, -1.75f));
 418       BOOST_CHECK(result.real() < 0);
 419       BOOST_CHECK(result.imag() < 0);
 420
 421       result = boost::math::acos(ct(zero, zero));
 422       BOOST_CHECK(result.real() > 0);
 423       BOOST_CHECK(result.imag() == 0);
 424       BOOST_CHECK((boost::math::signbit)(result.imag()));
 425       result = boost::math::acos(ct(zero, mzero));
 426       BOOST_CHECK(result.real() > 0);
 427       BOOST_CHECK(result.imag() == 0);
 428       BOOST_CHECK(0 == (boost::math::signbit)(result.imag()));
 429       result = boost::math::acos(ct(mzero, zero));
 430       BOOST_CHECK(result.real() > 0);
 431       BOOST_CHECK(result.imag() == 0);
 432       BOOST_CHECK((boost::math::signbit)(result.imag()));
 433       result = boost::math::acos(ct(mzero, mzero));
 434       BOOST_CHECK(result.real() > 0);
 435       BOOST_CHECK(result.imag() == 0);
 436       BOOST_CHECK(0 == (boost::math::signbit)(result.imag()));
 437    }
 438
 439    //
 440    // C99 spot tests for acosh:
 441    //
 442    result = boost::math::acosh(ct(zero, zero));
 443    BOOST_CHECK(result.real() == 0);
 444    BOOST_CHECK_CLOSE(result.imag(), half_pi, eps*200);
 445
 446    result = boost::math::acosh(ct(zero, mzero));
 447    BOOST_CHECK(result.real() == 0);
 448    BOOST_CHECK_CLOSE(result.imag(), -half_pi, eps*200);
 449
 450    result = boost::math::acosh(ct(mzero, zero));
 451    BOOST_CHECK(result.real() == 0);
 452    BOOST_CHECK_CLOSE(result.imag(), half_pi, eps*200);
 453
 454    result = boost::math::acosh(ct(mzero, mzero));
 455    BOOST_CHECK(result.real() == 0);
 456    BOOST_CHECK_CLOSE(result.imag(), -half_pi, eps*200);
 457
 458    if(test_infinity)
 459    {
 460       result = boost::math::acosh(ct(one, infinity));
 461       BOOST_CHECK(result.real() == infinity);
 462       BOOST_CHECK_CLOSE(result.imag(), half_pi, eps*200);
 463
 464       result = boost::math::acosh(ct(one, -infinity));
 465       BOOST_CHECK(result.real() == infinity);
 466       BOOST_CHECK_CLOSE(result.imag(), -half_pi, eps*200);
 467    }
 468
 469    if(test_nan)
 470    {
 471       result = boost::math::acosh(ct(one, nan));
 472       BOOST_CHECK((boost::math::isnan)(result.real()));
 473       BOOST_CHECK((boost::math::isnan)(result.imag()));
 474    }
 475    if(test_infinity)
 476    {
 477       result = boost::math::acosh(ct(-infinity, one));
 478       BOOST_CHECK(result.real() == infinity);
 479       BOOST_CHECK_CLOSE(result.imag(), pi, eps*200);
 480
 481       result = boost::math::acosh(ct(infinity, one));
 482       BOOST_CHECK(result.real() == infinity);
 483       BOOST_CHECK(result.imag() == 0);
 484
 485       result = boost::math::acosh(ct(-infinity, -one));
 486       BOOST_CHECK(result.real() == infinity);
 487       BOOST_CHECK_CLOSE(result.imag(), -pi, eps*200);
 488
 489       result = boost::math::acosh(ct(infinity, -one));
 490       BOOST_CHECK(result.real() == infinity);
 491       BOOST_CHECK(result.imag() == 0);
 492
 493       result = boost::math::acosh(ct(-infinity, infinity));
 494       BOOST_CHECK(result.real() == infinity);
 495       BOOST_CHECK_CLOSE(result.imag(), three_quarter_pi, eps*200);
 496
 497       result = boost::math::acosh(ct(infinity, infinity));
 498       BOOST_CHECK(result.real() == infinity);
 499       BOOST_CHECK_CLOSE(result.imag(), quarter_pi, eps*200);
 500
 501       result = boost::math::acosh(ct(-infinity, -infinity));
 502       BOOST_CHECK(result.real() == infinity);
 503       BOOST_CHECK_CLOSE(result.imag(), -three_quarter_pi, eps*200);
 504
 505       result = boost::math::acosh(ct(infinity, -infinity));
 506       BOOST_CHECK(result.real() == infinity);
 507       BOOST_CHECK_CLOSE(result.imag(), -quarter_pi, eps*200);
 508    }
 509
 510    if(test_nan)
 511    {
 512       result = boost::math::acosh(ct(infinity, nan));
 513       BOOST_CHECK(result.real() == infinity);
 514       BOOST_CHECK((boost::math::isnan)(result.imag()));
 515
 516       result = boost::math::acosh(ct(-infinity, nan));
 517       BOOST_CHECK(result.real() == infinity);
 518       BOOST_CHECK((boost::math::isnan)(result.imag()));
 519
 520       result = boost::math::acosh(ct(nan, one));
 521       BOOST_CHECK((boost::math::isnan)(result.real()));
 522       BOOST_CHECK((boost::math::isnan)(result.imag()));
 523
 524       result = boost::math::acosh(ct(nan, infinity));
 525       BOOST_CHECK(result.real() == infinity);
 526       BOOST_CHECK((boost::math::isnan)(result.imag()));
 527
 528       result = boost::math::acosh(ct(nan, -one));
 529       BOOST_CHECK((boost::math::isnan)(result.real()));
 530       BOOST_CHECK((boost::math::isnan)(result.imag()));
 531
 532       result = boost::math::acosh(ct(nan, -infinity));
 533       BOOST_CHECK(result.real() == infinity);
 534       BOOST_CHECK((boost::math::isnan)(result.imag()));
 535
 536       result = boost::math::acosh(ct(nan, nan));
 537       BOOST_CHECK((boost::math::isnan)(result.real()));
 538       BOOST_CHECK((boost::math::isnan)(result.imag()));
 539    }
 540    if(boost::math::signbit(mzero))
 541    {
 542       result = boost::math::acosh(ct(-2.5f, zero));
 543       BOOST_CHECK(result.real() > 0);
 544       BOOST_CHECK(result.imag() > 0);
 545    }
 546    //
 547    // C99 spot checks for asinh:
 548    //
 549    result = boost::math::asinh(ct(zero, zero));
 550    BOOST_CHECK(result.real() == 0);
 551    BOOST_CHECK(result.imag() == 0);
 552
 553    result = boost::math::asinh(ct(mzero, zero));
 554    BOOST_CHECK(result.real() == 0);
 555    BOOST_CHECK(result.imag() == 0);
 556
 557    result = boost::math::asinh(ct(zero, mzero));
 558    BOOST_CHECK(result.real() == 0);
 559    BOOST_CHECK(result.imag() == 0);
 560
 561    result = boost::math::asinh(ct(mzero, mzero));
 562    BOOST_CHECK(result.real() == 0);
 563    BOOST_CHECK(result.imag() == 0);
 564
 565    if(test_infinity)
 566    {
 567       result = boost::math::asinh(ct(one, infinity));
 568       BOOST_CHECK(result.real() == infinity);
 569       BOOST_CHECK_CLOSE(result.imag(), half_pi, eps*200);
 570
 571       result = boost::math::asinh(ct(one, -infinity));
 572       BOOST_CHECK(result.real() == infinity);
 573       BOOST_CHECK_CLOSE(result.imag(), -half_pi, eps*200);
 574
 575       result = boost::math::asinh(ct(-one, -infinity));
 576       BOOST_CHECK(result.real() == -infinity);
 577       BOOST_CHECK_CLOSE(result.imag(), -half_pi, eps*200);
 578
 579       result = boost::math::asinh(ct(-one, infinity));
 580       BOOST_CHECK(result.real() == -infinity);
 581       BOOST_CHECK_CLOSE(result.imag(), half_pi, eps*200);
 582    }
 583
 584    if(test_nan)
 585    {
 586       result = boost::math::asinh(ct(one, nan));
 587       BOOST_CHECK((boost::math::isnan)(result.real()));
 588       BOOST_CHECK((boost::math::isnan)(result.imag()));
 589
 590       result = boost::math::asinh(ct(-one, nan));
 591       BOOST_CHECK((boost::math::isnan)(result.real()));
 592       BOOST_CHECK((boost::math::isnan)(result.imag()));
 593
 594       result = boost::math::asinh(ct(zero, nan));
 595       BOOST_CHECK((boost::math::isnan)(result.real()));
 596       BOOST_CHECK((boost::math::isnan)(result.imag()));
 597    }
 598
 599    if(test_infinity)
 600    {
 601       result = boost::math::asinh(ct(infinity, one));
 602       BOOST_CHECK(result.real() == infinity);
 603       BOOST_CHECK(result.imag() == 0);
 604
 605       result = boost::math::asinh(ct(infinity, -one));
 606       BOOST_CHECK(result.real() == infinity);
 607       BOOST_CHECK(result.imag() == 0);
 608
 609       result = boost::math::asinh(ct(-infinity, -one));
 610       BOOST_CHECK(result.real() == -infinity);
 611       BOOST_CHECK(result.imag() == 0);
 612
 613       result = boost::math::asinh(ct(-infinity, one));
 614       BOOST_CHECK(result.real() == -infinity);
 615       BOOST_CHECK(result.imag() == 0);
 616
 617       result = boost::math::asinh(ct(infinity, infinity));
 618       BOOST_CHECK(result.real() == infinity);
 619       BOOST_CHECK_CLOSE(result.imag(), quarter_pi, eps*200);
 620
 621       result = boost::math::asinh(ct(infinity, -infinity));
 622       BOOST_CHECK(result.real() == infinity);
 623       BOOST_CHECK_CLOSE(result.imag(), -quarter_pi, eps*200);
 624
 625       result = boost::math::asinh(ct(-infinity, -infinity));
 626       BOOST_CHECK(result.real() == -infinity);
 627       BOOST_CHECK_CLOSE(result.imag(), -quarter_pi, eps*200);
 628
 629       result = boost::math::asinh(ct(-infinity, infinity));
 630       BOOST_CHECK(result.real() == -infinity);
 631       BOOST_CHECK_CLOSE(result.imag(), quarter_pi, eps*200);
 632    }
 633
 634    if(test_nan)
 635    {
 636       result = boost::math::asinh(ct(infinity, nan));
 637       BOOST_CHECK(result.real() == infinity);
 638       BOOST_CHECK((boost::math::isnan)(result.imag()));
 639
 640       result = boost::math::asinh(ct(-infinity, nan));
 641       BOOST_CHECK(result.real() == -infinity);
 642       BOOST_CHECK((boost::math::isnan)(result.imag()));
 643
 644       result = boost::math::asinh(ct(nan, zero));
 645       BOOST_CHECK((boost::math::isnan)(result.real()));
 646       BOOST_CHECK(result.imag() == 0);
 647
 648       result = boost::math::asinh(ct(nan,  mzero));
 649       BOOST_CHECK((boost::math::isnan)(result.real()));
 650       BOOST_CHECK(result.imag() == 0);
 651
 652       result = boost::math::asinh(ct(nan, one));
 653       BOOST_CHECK((boost::math::isnan)(result.real()));
 654       BOOST_CHECK((boost::math::isnan)(result.imag()));
 655
 656       result = boost::math::asinh(ct(nan,  -one));
 657       BOOST_CHECK((boost::math::isnan)(result.real()));
 658       BOOST_CHECK((boost::math::isnan)(result.imag()));
 659
 660       result = boost::math::asinh(ct(nan,  nan));
 661       BOOST_CHECK((boost::math::isnan)(result.real()));
 662       BOOST_CHECK((boost::math::isnan)(result.imag()));
 663
 664       result = boost::math::asinh(ct(nan, infinity));
 665       BOOST_CHECK(std::fabs(result.real()) == infinity);
 666       BOOST_CHECK((boost::math::isnan)(result.imag()));
 667
 668       result = boost::math::asinh(ct(nan,  -infinity));
 669       BOOST_CHECK(std::fabs(result.real()) == infinity);
 670       BOOST_CHECK((boost::math::isnan)(result.imag()));
 671    }
 672    if(boost::math::signbit(mzero))
 673    {
 674       result = boost::math::asinh(ct(zero, 1.5f));
 675       BOOST_CHECK(result.real() > 0);
 676       BOOST_CHECK(result.imag() > 0);
 677    }
 678
 679    //
 680    // C99 special cases for atanh:
 681    //
 682    result = boost::math::atanh(ct(zero, zero));
 683    BOOST_CHECK(result.real() == zero);
 684    BOOST_CHECK(result.imag() == zero);
 685
 686    result = boost::math::atanh(ct(mzero, zero));
 687    BOOST_CHECK(result.real() == zero);
 688    BOOST_CHECK(result.imag() == zero);
 689
 690    result = boost::math::atanh(ct(zero, mzero));
 691    BOOST_CHECK(result.real() == zero);
 692    BOOST_CHECK(result.imag() == zero);
 693
 694    result = boost::math::atanh(ct(mzero, mzero));
 695    BOOST_CHECK(result.real() == zero);
 696    BOOST_CHECK(result.imag() == zero);
 697
 698    if(test_nan)
 699    {
 700       result = boost::math::atanh(ct(zero, nan));
 701       BOOST_CHECK(result.real() == zero);
 702       BOOST_CHECK((boost::math::isnan)(result.imag()));
 703
 704       result = boost::math::atanh(ct(-zero, nan));
 705       BOOST_CHECK(result.real() == zero);
 706       BOOST_CHECK((boost::math::isnan)(result.imag()));
 707    }
 708
 709    if(test_infinity)
 710    {
 711       result = boost::math::atanh(ct(one, zero));
 712       BOOST_CHECK_EQUAL(result.real(), infinity);
 713       BOOST_CHECK_EQUAL(result.imag(), zero);
 714
 715       result = boost::math::atanh(ct(-one, zero));
 716       BOOST_CHECK_EQUAL(result.real(), -infinity);
 717       BOOST_CHECK_EQUAL(result.imag(), zero);
 718
 719       result = boost::math::atanh(ct(-one, -zero));
 720       BOOST_CHECK_EQUAL(result.real(), -infinity);
 721       BOOST_CHECK_EQUAL(result.imag(), zero);
 722
 723       result = boost::math::atanh(ct(one, -zero));
 724       BOOST_CHECK_EQUAL(result.real(), infinity);
 725       BOOST_CHECK_EQUAL(result.imag(), zero);
 726
 727       result = boost::math::atanh(ct(pi, infinity));
 728       BOOST_CHECK_EQUAL(result.real(), zero);
 729       BOOST_CHECK_CLOSE(result.imag(), half_pi, eps*200);
 730
 731       result = boost::math::atanh(ct(pi, -infinity));
 732       BOOST_CHECK_EQUAL(result.real(), zero);
 733       BOOST_CHECK_CLOSE(result.imag(), -half_pi, eps*200);
 734
 735       result = boost::math::atanh(ct(-pi, -infinity));
 736       BOOST_CHECK_EQUAL(result.real(), zero);
 737       BOOST_CHECK_CLOSE(result.imag(), -half_pi, eps*200);
 738
 739       result = boost::math::atanh(ct(-pi, infinity));
 740       BOOST_CHECK_EQUAL(result.real(), zero);
 741       BOOST_CHECK_CLOSE(result.imag(), half_pi, eps*200);
 742    }
 743    if(test_nan)
 744    {
 745       result = boost::math::atanh(ct(pi, nan));
 746       BOOST_CHECK((boost::math::isnan)(result.real()));
 747       BOOST_CHECK((boost::math::isnan)(result.imag()));
 748
 749       result = boost::math::atanh(ct(-pi, nan));
 750       BOOST_CHECK((boost::math::isnan)(result.real()));
 751       BOOST_CHECK((boost::math::isnan)(result.imag()));
 752    }
 753
 754    if(test_infinity)
 755    {
 756       result = boost::math::atanh(ct(infinity, pi));
 757       BOOST_CHECK(result.real() == zero);
 758       BOOST_CHECK_CLOSE(result.imag(), half_pi, eps*200);
 759
 760       result = boost::math::atanh(ct(infinity, -pi));
 761       BOOST_CHECK_EQUAL(result.real(), zero);
 762       BOOST_CHECK_CLOSE(result.imag(), -half_pi, eps*200);
 763
 764       result = boost::math::atanh(ct(-infinity, -pi));
 765       BOOST_CHECK_EQUAL(result.real(), zero);
 766       BOOST_CHECK_CLOSE(result.imag(), -half_pi, eps*200);
 767
 768       result = boost::math::atanh(ct(-infinity, pi));
 769       BOOST_CHECK_EQUAL(result.real(), zero);
 770       BOOST_CHECK_CLOSE(result.imag(), half_pi, eps*200);
 771
 772       result = boost::math::atanh(ct(infinity, infinity));
 773       BOOST_CHECK_EQUAL(result.real(), zero);
 774       BOOST_CHECK_CLOSE(result.imag(), half_pi, eps*200);
 775
 776       result = boost::math::atanh(ct(infinity, -infinity));
 777       BOOST_CHECK_EQUAL(result.real(), zero);
 778       BOOST_CHECK_CLOSE(result.imag(), -half_pi, eps*200);
 779
 780       result = boost::math::atanh(ct(-infinity, -infinity));
 781       BOOST_CHECK_EQUAL(result.real(), zero);
 782       BOOST_CHECK_CLOSE(result.imag(), -half_pi, eps*200);
 783
 784       result = boost::math::atanh(ct(-infinity, infinity));
 785       BOOST_CHECK_EQUAL(result.real(), zero);
 786       BOOST_CHECK_CLOSE(result.imag(), half_pi, eps*200);
 787    }
 788
 789    if(test_nan)
 790    {
 791       result = boost::math::atanh(ct(infinity, nan));
 792       BOOST_CHECK(result.real() == 0);
 793       BOOST_CHECK((boost::math::isnan)(result.imag()));
 794
 795       result = boost::math::atanh(ct(-infinity, nan));
 796       BOOST_CHECK(result.real() == 0);
 797       BOOST_CHECK((boost::math::isnan)(result.imag()));
 798
 799       result = boost::math::atanh(ct(nan, pi));
 800       BOOST_CHECK((boost::math::isnan)(result.real()));
 801       BOOST_CHECK((boost::math::isnan)(result.imag()));
 802
 803       result = boost::math::atanh(ct(nan, -pi));
 804       BOOST_CHECK((boost::math::isnan)(result.real()));
 805       BOOST_CHECK((boost::math::isnan)(result.imag()));
 806
 807       result = boost::math::atanh(ct(nan, infinity));
 808       BOOST_CHECK(result.real() == 0);
 809       BOOST_CHECK_CLOSE(result.imag(), half_pi, eps*200);
 810
 811       result = boost::math::atanh(ct(nan, -infinity));
 812       BOOST_CHECK(result.real() == 0);
 813       BOOST_CHECK_CLOSE(result.imag(), -half_pi, eps*200);
 814
 815       result = boost::math::atanh(ct(nan, nan));
 816       BOOST_CHECK((boost::math::isnan)(result.real()));
 817       BOOST_CHECK((boost::math::isnan)(result.imag()));
 818
 819    }
 820    if(boost::math::signbit(mzero))
 821    {
 822       result = boost::math::atanh(ct(-2.0f, mzero));
 823       BOOST_CHECK(result.real() < 0);
 824       BOOST_CHECK(result.imag() < 0);
 825    }
 826 }
 827
 828 //
 829 // test_boundaries:
 830 // This is an accuracy test, sets the real and imaginary components
 831 // of the input argument to various "boundary conditions" that exist
 832 // inside the implementation.  Then computes the result at double precision
 833 // and again at float precision.  The double precision result will be
 834 // computed using the "regular" code, where as the float precision versions
 835 // will calculate the result using the "exceptional value" handlers, so
 836 // we end up comparing the values calculated by two different methods.
 837 //
 838 const float boundaries[] = {
 839    0,
 840    1,
 841    2,
 842    (std::numeric_limits<float>::max)(),
 843    (std::numeric_limits<float>::min)(),
 844    std::numeric_limits<float>::epsilon(),
 845    std::sqrt((std::numeric_limits<float>::max)()) / 8,
 846    static_cast<float>(4) * std::sqrt((std::numeric_limits<float>::min)()),
 847    0.6417F,
 848    1.5F,
 849    std::sqrt((std::numeric_limits<float>::max)()) / 2,
 850    std::sqrt((std::numeric_limits<float>::min)()),
 851    1.0F / 0.3F,
 852 };
 853
 854 void do_test_boundaries(float x, float y)
 855 {
 856    std::complex<float> r1 = boost::math::asin(std::complex<float>(x, y));
 857    std::complex<double> dr = boost::math::asin(std::complex<double>(x, y));
 858    std::complex<float> r2(static_cast<float>(dr.real()), static_cast<float>(dr.imag()));
 859    check_complex(r2, r1, 5);
 860    r1 = boost::math::acos(std::complex<float>(x, y));
 861    dr = boost::math::acos(std::complex<double>(x, y));
 862    r2 = std::complex<float>(std::complex<double>(dr.real(), dr.imag()));
 863    check_complex(r2, r1, 5);
 864    r1 = boost::math::atanh(std::complex<float>(x, y));
 865    dr = boost::math::atanh(std::complex<double>(x, y));
 866    r2 = std::complex<float>(std::complex<double>(dr.real(), dr.imag()));
 867    check_complex(r2, r1, 5);
 868 }
 869
 870 void test_boundaries(float x, float y)
 871 {
 872    do_test_boundaries(x, y);
 873    do_test_boundaries(-x, y);
 874    do_test_boundaries(-x, -y);
 875    do_test_boundaries(x, -y);
 876 }
 877
 878 void test_boundaries(float x)
 879 {
 880    for(unsigned i = 0; i < sizeof(boundaries)/sizeof(float); ++i)
 881    {
 882       test_boundaries(x, boundaries[i]);
 883       test_boundaries(x, boundaries[i] + std::numeric_limits<float>::epsilon()*boundaries[i]);
 884       test_boundaries(x, boundaries[i] - std::numeric_limits<float>::epsilon()*boundaries[i]);
 885    }
 886 }
 887
 888 void test_boundaries()
 889 {
 890    for(unsigned i = 0; i < sizeof(boundaries)/sizeof(float); ++i)
 891    {
 892       test_boundaries(boundaries[i]);
 893       test_boundaries(boundaries[i] + std::numeric_limits<float>::epsilon()*boundaries[i]);
 894       test_boundaries(boundaries[i] - std::numeric_limits<float>::epsilon()*boundaries[i]);//here
 895    }
 896 }
 897
 898
 899 BOOST_AUTO_TEST_CASE( test_main )
 900 {
 901    std::cout << "Running complex trig sanity checks for type float." << std::endl;
 902    test_inverse_trig(float(0));
 903    std::cout << "Running complex trig sanity checks for type double." << std::endl;
 904    test_inverse_trig(double(0));
 905    //test_inverse_trig((long double)(0));
 906
 907    std::cout << "Running complex trig spot checks for type float." << std::endl;
 908    check_spots(float(0));
 909    std::cout << "Running complex trig spot checks for type double." << std::endl;
 910    check_spots(double(0));
 911 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 912    std::cout << "Running complex trig spot checks for type long double." << std::endl;
 913    check_spots((long double)(0));
 914 #endif
 915
 916    std::cout << "Running complex trig boundary and accuracy tests." << std::endl;
 917    test_boundaries();
 918 }
 919
 920
 921