//  (C) Copyright John Maddock 2005.

//  Use, modification and distribution are subject to the

//  Boost Software License, Version 1.0. (See accompanying file

//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

/*

 * © Portions copyright (c) 2006-2007 Nokia Corporation.  All rights reserved.

*/

#include <boost/test/test_tools.hpp>

#include <boost/test/included/test_exec_monitor.hpp>

#include <boost/test/floating_point_comparison.hpp>

#include <boost/type_traits/is_same.hpp>

#include <boost/type_traits/is_floating_point.hpp>

#include <boost/mpl/if.hpp>

#include <boost/static_assert.hpp>

#include <boost/math/complex.hpp>

#include <iostream>

#include <iomanip>

#include <cmath>

#include <typeinfo>

#ifdef BOOST_NO_STDC_NAMESPACE 

namespace std{ using ::sqrt; using ::tan; using ::tanh; }

#endif

#ifndef VERBOSE

#undef BOOST_MESSAGE

#define BOOST_MESSAGE(x)

#endif

//

// check_complex:

// Verifies that expected value "a" and found value "b" have a relative error

// less than "max_error" epsilons.  Note that relative error is calculated for

// the complex number as a whole; this means that the error in the real or 

// imaginary parts alone can be much higher than max_error when the real and 

// imaginary parts are of very different magnitudes.  This is important, because

// the Hull et al analysis of the acos and asin algorithms requires that very small

// real/imaginary components can be safely ignored if they are negligible compared

// to the other component.

//

template <class T>

bool check_complex(const std::complex<T>& a, const std::complex<T>& b, int max_error)

{

   //

   // a is the expected value, b is what was actually found,

   // compute | (a-b)/b | and compare with max_error which is the 

   // multiple of E to permit:

   //

   bool result = true;

   static const std::complex<T> zero(0);

   static const T eps = std::pow(static_cast<T>(std::numeric_limits<T>::radix), 1 - std::numeric_limits<T>::digits);

   if(a == zero)

   {

      if(b != zero)

      {

         if(boost::math::fabs(b) > eps)

         {

            result = false;

            BOOST_ERROR("Expected {0,0} but got: " << b);

         }

         else

         {

            BOOST_MESSAGE("Expected {0,0} but got: " << b);

         }

      }

      return result;

   }

   else if(b == zero)

   {

      if(boost::math::fabs(a) > eps)

      {

         BOOST_ERROR("Found {0,0} but expected: " << a);

         return false;;

      }

      else

      {

         BOOST_MESSAGE("Found {0,0} but expected: " << a);

      }

   }

   T rel = boost::math::fabs((b-a)/b) / eps;

   if( rel > max_error)

   {

      result = false;

      BOOST_ERROR("Error in result exceeded permitted limit of " << max_error << " (actual relative error was " << rel << "e).  Found " << b << " expected " << a);

   }

   return result;

}

//

// test_inverse_trig:

// This is nothing more than a sanity check, computes trig(atrig(z)) 

// and compare the result to z.  Note that:

//

// atrig(trig(z)) != z

//

// for certain z because the inverse trig functions are multi-valued, this 

// essentially rules this out as a testing method.  On the other hand:

//

// trig(atrig(z))

//

// can vary compare to z by an arbitrarily large amount.  For one thing we 

// have no control over the implementation of the trig functions, for another

// even if both functions were accurate to 1ulp (as accurate as transcendental

// number can get, thanks to the "table makers dilemma"), the errors can still

// be arbitrarily large - often the inverse trig functions will map a very large

// part of the complex domain into a small output domain, so you can never get

// back exactly where you started from.  Consequently these tests are no more than

// sanity checks (just verifies that signs are correct and so on).

//

template <class T>

void test_inverse_trig(T)

{

   using namespace std;

   static const T interval = static_cast<T>(2.0L/128.0L);

   T x, y;

   std::cout << std::setprecision(std::numeric_limits<T>::digits10+2);

   for(x = -1; x <= 1; x += interval)

   {

      for(y = -1; y <= 1; y += interval)

      {

         // acos:

         std::complex<T> val(x, y), inter, result;

         inter = boost::math::acos(val);

         result = cos(inter);

         if(!check_complex(val, result, 50))

         {

            std::cout << "Error in testing inverse complex cos for type " << typeid(T).name() << std::endl;

            std::cout << "   val=             " << val << std::endl;

            std::cout << "   acos(val) =      " << inter << std::endl;

            std::cout << "   cos(acos(val)) = " << result << std::endl;

         }

         // asin:

         inter = boost::math::asin(val);

         result = sin(inter);

         if(!check_complex(val, result, 5))

         {

            std::cout << "Error in testing inverse complex sin for type " << typeid(T).name() << std::endl;

            std::cout << "   val=             " << val << std::endl;

            std::cout << "   asin(val) =      " << inter << std::endl;

            std::cout << "   sin(asin(val)) = " << result << std::endl;

         }

      }

   }

   static const T interval2 = static_cast<T>(3.0L/256.0L);

   for(x = -3; x <= 3; x += interval2)

   {

      for(y = -3; y <= 3; y += interval2)

      {

         // asinh:

         std::complex<T> val(x, y), inter, result;

         inter = boost::math::asinh(val);

         result = sinh(inter);

         if(!check_complex(val, result, 5))

         {

            std::cout << "Error in testing inverse complex sinh for type " << typeid(T).name() << std::endl;

            std::cout << "   val=               " << val << std::endl;

            std::cout << "   asinh(val) =       " << inter << std::endl;

            std::cout << "   sinh(asinh(val)) = " << result << std::endl;

         }

         // acosh:

         if(!((y == 0) && (x <= 1))) // can't test along the branch cut

         {

            inter = boost::math::acosh(val);

            result = cosh(inter);

            if(!check_complex(val, result, 60))

            {

               std::cout << "Error in testing inverse complex cosh for type " << typeid(T).name() << std::endl;

               std::cout << "   val=               " << val << std::endl;

               std::cout << "   acosh(val) =       " << inter << std::endl;

               std::cout << "   cosh(acosh(val)) = " << result << std::endl;

            }

         }

         //

         // There is a problem in testing atan and atanh:

         // The inverse functions map a large input range to a much

         // smaller output range, so at the extremes too rather different

         // inputs may map to the same output value once rounded to N places.

         // Consequently tan(atan(z)) can suffer from arbitrarily large errors

         // even if individually they each have a small error bound.  On the other

         // hand we can't test atan(tan(z)) either because atan is multi-valued, so

         // round-tripping in this direction isn't always possible.

         // The following heuristic is designed to make the best of a bad job,

         // using atan(tan(z)) where possible and tan(atan(z)) when it's not.

         //

         static const int tanh_error = 20;

         if((0 != x) && (0 != y) && ((std::fabs(y) < 1) || (std::fabs(x) < 1)))

         {

            // atanh:

            val = boost::math::atanh(val);

            inter = tanh(val);

            result = boost::math::atanh(inter);

            if(!check_complex(val, result, tanh_error))

            {

               std::cout << "Error in testing inverse complex tanh for type " << typeid(T).name() << std::endl;

               std::cout << "   val=               " << val << std::endl;

               std::cout << "   tanh(val) =        " << inter << std::endl;

               std::cout << "   atanh(tanh(val)) = " << result << std::endl;

            }

            // atan:

            if(!((x == 0) && (std::fabs(y) == 1))) // we can't test infinities here

            {

               val = std::complex<T>(x, y);

               val = boost::math::atan(val);

               inter = tan(val);

               result = boost::math::atan(inter);

               if(!check_complex(val, result, tanh_error))

               {

                  std::cout << "Error in testing inverse complex tan for type " << typeid(T).name() << std::endl;

                  std::cout << "   val=               " << val << std::endl;

                  std::cout << "   tan(val) =         " << inter << std::endl;

                  std::cout << "   atan(tan(val)) =   " << result << std::endl;

               }

            }

         }

         else

         {

            // atanh:

            inter = boost::math::atanh(val);

            result = tanh(inter);

            if(!check_complex(val, result, tanh_error))

            {

               std::cout << "Error in testing inverse complex atanh for type " << typeid(T).name() << std::endl;

               std::cout << "   val=                 " << val << std::endl;

               std::cout << "   atanh(val) =         " << inter << std::endl;

               std::cout << "   tanh(atanh(val)) =   " << result << std::endl;

            }

            // atan:

            if(!((x == 0) && (std::fabs(y) == 1))) // we can't test infinities here

            {

               inter = boost::math::atan(val);

               result = tan(inter);

               if(!check_complex(val, result, tanh_error))

               {

                  std::cout << "Error in testing inverse complex atan for type " << typeid(T).name() << std::endl;

                  std::cout << "   val=                 " << val << std::endl;

                  std::cout << "   atan(val) =          " << inter << std::endl;

                  std::cout << "   tan(atan(val)) =     " << result << std::endl;

               }

            }

         }

      }

   }

}

//

// check_spots:

// Various spot values, mostly the C99 special cases (infinites and NAN's).

// TODO: add spot checks for the Wolfram spot values.

//

template <class T>

void check_spots(const T&)

{

   typedef std::complex<T> ct;

   ct result;

   static const T two = 2.0;

   T eps = std::pow(two, 1-std::numeric_limits<T>::digits); // numeric_limits<>::epsilon way too small to be useful on Darwin.

   static const T zero = 0;

   static const T mzero = -zero;

   static const T one = 1;

   static const T pi = static_cast<T>(3.141592653589793238462643383279502884197L);

   static const T half_pi = static_cast<T>(1.57079632679489661923132169163975144L);

   static const T quarter_pi = static_cast<T>(0.78539816339744830961566084581987572L);

   static const T three_quarter_pi = static_cast<T>(2.35619449019234492884698253745962716L);

   //static const T log_two = static_cast<T>(0.69314718055994530941723212145817657L);

   T infinity = std::numeric_limits<T>::infinity();

   bool test_infinity = std::numeric_limits<T>::has_infinity;

   T nan = 0;

   bool test_nan = false;

#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x564))

   // numeric_limits reports that a quiet NaN is present

   // but an attempt to access it will terminate the program!!!!

   if(std::numeric_limits<T>::has_quiet_NaN)

      nan = std::numeric_limits<T>::quiet_NaN();

   if(boost::math::detail::test_is_nan(nan))

      test_nan = true;

#endif

#if defined(__DECCXX) && !defined(_IEEE_FP)

   // Tru64 cxx traps infinities unless the -ieee option is used:

   test_infinity = false;

#endif

   //

   // C99 spot tests for acos:

   //

   result = boost::math::acos(ct(zero));

   check_complex(ct(half_pi), result, 2);

   result = boost::math::acos(ct(mzero));

   check_complex(ct(half_pi), result, 2);

   result = boost::math::acos(ct(zero, mzero));

   check_complex(ct(half_pi), result, 2);

   result = boost::math::acos(ct(mzero, mzero));

   check_complex(ct(half_pi), result, 2);

   if(test_nan)

   {

      result = boost::math::acos(ct(zero,nan));

      BOOST_CHECK_CLOSE(result.real(), half_pi, eps*200);

      BOOST_CHECK(boost::math::detail::test_is_nan(result.imag()));

      result = boost::math::acos(ct(mzero,nan));

      BOOST_CHECK_CLOSE(result.real(), half_pi, eps*200);

      BOOST_CHECK(boost::math::detail::test_is_nan(result.imag()));

   }

   if(test_infinity)

   {

      result = boost::math::acos(ct(zero, infinity));

      BOOST_CHECK_CLOSE(result.real(), half_pi, eps*200);

      BOOST_CHECK(result.imag() == -infinity);

      result = boost::math::acos(ct(zero, -infinity));

      BOOST_CHECK_CLOSE(result.real(), half_pi, eps*200);

      BOOST_CHECK(result.imag() == infinity);

   }

   if(test_nan)

   {

      result = boost::math::acos(ct(one, nan));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.real()));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.imag()));

   }

   if(test_infinity)

   {

      result = boost::math::acos(ct(-infinity, one));

      BOOST_CHECK_CLOSE(result.real(), pi, eps*200);

      BOOST_CHECK(result.imag() == -infinity);

      result = boost::math::acos(ct(infinity, one));

      BOOST_CHECK(result.real() == 0);

      BOOST_CHECK(result.imag() == -infinity);

      result = boost::math::acos(ct(-infinity, -one));

      BOOST_CHECK_CLOSE(result.real(), pi, eps*200);

      BOOST_CHECK(result.imag() == infinity);

      result = boost::math::acos(ct(infinity, -one));

      BOOST_CHECK(result.real() == 0);

      BOOST_CHECK(result.imag() == infinity);

      result = boost::math::acos(ct(-infinity, infinity));

      BOOST_CHECK_CLOSE(result.real(), three_quarter_pi, eps*200);

      BOOST_CHECK(result.imag() == -infinity);

      result = boost::math::acos(ct(infinity, infinity));

      BOOST_CHECK_CLOSE(result.real(), quarter_pi, eps*200);

      BOOST_CHECK(result.imag() == -infinity);

      result = boost::math::acos(ct(-infinity, -infinity));

      BOOST_CHECK_CLOSE(result.real(), three_quarter_pi, eps*200);

      BOOST_CHECK(result.imag() == infinity);

      result = boost::math::acos(ct(infinity, -infinity));

      BOOST_CHECK_CLOSE(result.real(), quarter_pi, eps*200);

      BOOST_CHECK(result.imag() == infinity);

   }

   if(test_nan)

   {

      result = boost::math::acos(ct(infinity, nan));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.real()));

      BOOST_CHECK(std::fabs(result.imag()) == infinity);

      result = boost::math::acos(ct(-infinity, nan));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.real()));

      BOOST_CHECK(std::fabs(result.imag()) == infinity);

      result = boost::math::acos(ct(nan, zero));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.real()));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.imag()));

      result = boost::math::acos(ct(nan, -zero));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.real()));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.imag()));

      result = boost::math::acos(ct(nan, one));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.real()));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.imag()));

      result = boost::math::acos(ct(nan, -one));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.real()));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.imag()));

      result = boost::math::acos(ct(nan, nan));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.real()));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.imag()));

      result = boost::math::acos(ct(nan, infinity));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.real()));

      BOOST_CHECK(result.imag() == -infinity);

      result = boost::math::acos(ct(nan, -infinity));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.real()));

      BOOST_CHECK(result.imag() == infinity);

   }

   //

   // C99 spot tests for acosh:

   //

   result = boost::math::acosh(ct(zero, zero));

   BOOST_CHECK(result.real() == 0);

   BOOST_CHECK_CLOSE(result.imag(), half_pi, eps*200);

   result = boost::math::acosh(ct(zero, mzero));

   BOOST_CHECK(result.real() == 0);

   BOOST_CHECK_CLOSE(result.imag(), half_pi, eps*200);

   result = boost::math::acosh(ct(mzero, zero));

   BOOST_CHECK(result.real() == 0);

   BOOST_CHECK_CLOSE(result.imag(), half_pi, eps*200);

   result = boost::math::acosh(ct(mzero, mzero));

   BOOST_CHECK(result.real() == 0);

   BOOST_CHECK_CLOSE(result.imag(), half_pi, eps*200);

   if(test_infinity)

   {

      result = boost::math::acosh(ct(one, infinity));

      BOOST_CHECK(result.real() == infinity);

      BOOST_CHECK_CLOSE(result.imag(), half_pi, eps*200);

      result = boost::math::acosh(ct(one, -infinity));

      BOOST_CHECK(result.real() == infinity);

      BOOST_CHECK_CLOSE(result.imag(), -half_pi, eps*200);

   }

   if(test_nan)

   {

      result = boost::math::acosh(ct(one, nan));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.real()));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.imag()));

   }

   if(test_infinity)

   {

      result = boost::math::acosh(ct(-infinity, one));

      BOOST_CHECK(result.real() == infinity);

      BOOST_CHECK_CLOSE(result.imag(), pi, eps*200);

      result = boost::math::acosh(ct(infinity, one));

      BOOST_CHECK(result.real() == infinity);

      BOOST_CHECK(result.imag() == 0);

      result = boost::math::acosh(ct(-infinity, -one));

      BOOST_CHECK(result.real() == infinity);

      BOOST_CHECK_CLOSE(result.imag(), -pi, eps*200);

      result = boost::math::acosh(ct(infinity, -one));

      BOOST_CHECK(result.real() == infinity);

      BOOST_CHECK(result.imag() == 0);

      result = boost::math::acosh(ct(-infinity, infinity));

      BOOST_CHECK(result.real() == infinity);

      BOOST_CHECK_CLOSE(result.imag(), three_quarter_pi, eps*200);

      result = boost::math::acosh(ct(infinity, infinity));

      BOOST_CHECK(result.real() == infinity);

      BOOST_CHECK_CLOSE(result.imag(), quarter_pi, eps*200);

      result = boost::math::acosh(ct(-infinity, -infinity));

      BOOST_CHECK(result.real() == infinity);

      BOOST_CHECK_CLOSE(result.imag(), -three_quarter_pi, eps*200);

      result = boost::math::acosh(ct(infinity, -infinity));

      BOOST_CHECK(result.real() == infinity);

      BOOST_CHECK_CLOSE(result.imag(), -quarter_pi, eps*200);

   }

   if(test_nan)

   {

      result = boost::math::acosh(ct(infinity, nan));

      BOOST_CHECK(result.real() == infinity);

      BOOST_CHECK(boost::math::detail::test_is_nan(result.imag()));

      result = boost::math::acosh(ct(-infinity, nan));

      BOOST_CHECK(result.real() == infinity);

      BOOST_CHECK(boost::math::detail::test_is_nan(result.imag()));

      result = boost::math::acosh(ct(nan, one));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.real()));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.imag()));

      result = boost::math::acosh(ct(nan, infinity));

      BOOST_CHECK(result.real() == infinity);

      BOOST_CHECK(boost::math::detail::test_is_nan(result.imag()));

      result = boost::math::acosh(ct(nan, -one));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.real()));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.imag()));

      result = boost::math::acosh(ct(nan, -infinity));

      BOOST_CHECK(result.real() == infinity);

      BOOST_CHECK(boost::math::detail::test_is_nan(result.imag()));

      result = boost::math::acosh(ct(nan, nan));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.real()));

      BOOST_CHECK(boost::math::detail::test_is_nan(result.imag()));

   }

   //

   // C99 spot checks for asinh:

   //

   result = boost::math::asinh(ct(zero, zero));

   BOOST_CHECK(result.real() == 0);

   BOOST_CHECK(result.imag() == 0);

   result = boost::math::asinh(ct(mzero, zero));

   BOOST_CHECK(result.real() == 0);

   BOOST_CHECK(result.imag() == 0);

   result = boost::math::asinh(ct(zero, mzero));

   BOOST_CHECK(result.real() == 0);

   BOOST_CHECK(result.imag() == 0);

   result = boost::math::asinh(ct(mzero, mzero));

   BOOST_CHECK(result.real() == 0);

   BOOST_CHECK(result.imag() == 0);

   if(test_infinity)

   {

      result = boost::math::asinh(ct(one, infinity));

      BOOST_CHECK(result.real() == infinity);