libs/math/test/test_recurrence.cpp

   1 //  (C) Copyright John Maddock 2018.
   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_MODULE test_recurrences
   7
   8 #include <boost/config.hpp>
   9
  10 #ifndef BOOST_NO_CXX11_HDR_TUPLE
  11 #include <boost/multiprecision/cpp_bin_float.hpp>
  12 #include <boost/math/tools/recurrence.hpp>
  13 #include <boost/math/special_functions/bessel.hpp>
  14 #include <boost/test/included/unit_test.hpp>
  15 #include <boost/test/floating_point_comparison.hpp>
  16 //#include <boost/test/tools/floating_point_comparison.hpp>
  17 #include <boost/math/concepts/real_concept.hpp>
  18
  19 #ifdef BOOST_MSVC
  20 #pragma warning(disable:4127)
  21 #endif
  22
  23 template <class T>
  24 struct bessel_jy_recurrence
  25 {
  26    bessel_jy_recurrence(T v, T z) : v(v), z(z) {}
  27    boost::math::tuple<T, T, T> operator()(int k)const
  28    {
  29       return boost::math::tuple<T, T, T>(T(1), -2 * (v + k) / z, T(1));
  30    }
  31
  32    T v, z;
  33 };
  34
  35 template <class T>
  36 struct bessel_ik_recurrence
  37 {
  38    bessel_ik_recurrence(T v, T z) : v(v), z(z) {}
  39    boost::math::tuple<T, T, T> operator()(int k)const
  40    {
  41       return boost::math::tuple<T, T, T>(T(1), -2 * (v + k) / z, T(-1));
  42    }
  43
  44    T v, z;
  45 };
  46
  47
  48 template <class T>
  49 void test_spots(T, const char* name)
  50 {
  51    std::cout << "Running tests for type " << name << std::endl;
  52    T tol = boost::math::tools::epsilon<T>() * 5;
  53    if ((std::numeric_limits<T>::digits > 53) || (std::numeric_limits<T>::digits == 0))
  54       tol *= 5;
  55    //
  56    // Test forward recurrence on Y_v(x):
  57    //
  58    {
  59       T v = 22.25;
  60       T x = 4.125;
  61       bessel_jy_recurrence<T> coef(v, x);
  62       T prev;
  63       T first = boost::math::cyl_neumann(v - 1, x);
  64       T second = boost::math::cyl_neumann(v, x);
  65       T sixth = boost::math::tools::apply_recurrence_relation_forward(coef, 6, first, second, (int*)0, &prev);
  66       T expected1 = boost::math::cyl_neumann(v + 6, x);
  67       T expected2 = boost::math::cyl_neumann(v + 5, x);
  68       BOOST_CHECK_CLOSE_FRACTION(sixth, expected1, tol);
  69       BOOST_CHECK_CLOSE_FRACTION(prev, expected2, tol);
  70
  71       boost::math::tools::forward_recurrence_iterator< bessel_jy_recurrence<T> > it(coef, first, second);
  72       for (unsigned i = 0; i < 15; ++i)
  73       {
  74          expected1 = boost::math::cyl_neumann(v + i, x);
  75          T found = *it;
  76          BOOST_CHECK_CLOSE_FRACTION(found, expected1, tol);
  77          ++it;
  78       }
  79
  80       if (std::numeric_limits<T>::max_exponent > 300)
  81       {
  82          //
  83          // This calculates the ratio Y_v(x)/Y_v+1(x) from the recurrence relations
  84          // which are only transiently stable since Y_v is not minimal as v->-INF
  85          // but only as v->0.  We have to be sure that v is sufficiently large that
  86          // convergence is complete before we reach the origin.
  87          //
  88          v = 102.75;
  89          boost::uintmax_t max_iter = 200;
  90          T ratio = boost::math::tools::function_ratio_from_forwards_recurrence(bessel_jy_recurrence<T>(v, x), boost::math::tools::epsilon<T>(), max_iter);
  91          first = boost::math::cyl_neumann(v, x);
  92          second = boost::math::cyl_neumann(v + 1, x);
  93          BOOST_CHECK_CLOSE_FRACTION(ratio, first / second, tol);
  94
  95          boost::math::tools::forward_recurrence_iterator< bessel_jy_recurrence<T> > it2(bessel_jy_recurrence<T>(v, x), boost::math::cyl_neumann(v, x));
  96          for (unsigned i = 0; i < 15; ++i)
  97          {
  98             expected1 = boost::math::cyl_neumann(v + i, x);
  99             T found = *it2;
 100             BOOST_CHECK_CLOSE_FRACTION(found, expected1, tol);
 101             ++it2;
 102          }
 103       }
 104
 105    }
 106    //
 107    // Test backward recurrence on J_v(x):
 108    //
 109    {
 110       if ((std::numeric_limits<T>::digits > 53) || !std::numeric_limits<T>::is_specialized)
 111          tol *= 5;
 112
 113       T v = 22.25;
 114       T x = 4.125;
 115       bessel_jy_recurrence<T> coef(v, x);
 116       T prev;
 117       T first = boost::math::cyl_bessel_j(v + 1, x);
 118       T second = boost::math::cyl_bessel_j(v, x);
 119       T sixth = boost::math::tools::apply_recurrence_relation_backward(coef, 6, first, second, (int*)0, &prev);
 120       T expected1 = boost::math::cyl_bessel_j(v - 6, x);
 121       T expected2 = boost::math::cyl_bessel_j(v - 5, x);
 122       BOOST_CHECK_CLOSE_FRACTION(sixth, expected1, tol);
 123       BOOST_CHECK_CLOSE_FRACTION(prev, expected2, tol);
 124
 125       boost::math::tools::backward_recurrence_iterator< bessel_jy_recurrence<T> > it(coef, first, second);
 126       for (unsigned i = 0; i < 15; ++i)
 127       {
 128          expected1 = boost::math::cyl_bessel_j(v - i, x);
 129          T found = *it;
 130          BOOST_CHECK_CLOSE_FRACTION(found, expected1, tol);
 131          ++it;
 132       }
 133
 134       boost::uintmax_t max_iter = 200;
 135       T ratio = boost::math::tools::function_ratio_from_backwards_recurrence(bessel_jy_recurrence<T>(v, x), boost::math::tools::epsilon<T>(), max_iter);
 136       first = boost::math::cyl_bessel_j(v, x);
 137       second = boost::math::cyl_bessel_j(v - 1, x);
 138       BOOST_CHECK_CLOSE_FRACTION(ratio, first / second, tol);
 139
 140       boost::math::tools::backward_recurrence_iterator< bessel_jy_recurrence<T> > it2(bessel_jy_recurrence<T>(v, x), boost::math::cyl_bessel_j(v, x));
 141       //boost::math::tools::backward_recurrence_iterator< bessel_jy_recurrence<T> > it3(bessel_jy_recurrence<T>(v, x), boost::math::cyl_neumann(v+1, x), boost::math::cyl_neumann(v, x));
 142       for (unsigned i = 0; i < 15; ++i)
 143       {
 144          expected1 = boost::math::cyl_bessel_j(v - i, x);
 145          T found = *it2;
 146          BOOST_CHECK_CLOSE_FRACTION(found, expected1, tol);
 147          ++it2;
 148       }
 149
 150    }
 151 }
 152
 153
 154 BOOST_AUTO_TEST_CASE( test_main )
 155 {
 156    BOOST_MATH_CONTROL_FP;
 157 #if !defined(TEST) || TEST == 1
 158    test_spots(0.0F, "float");
 159    test_spots(0.0, "double");
 160 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 161    test_spots(0.0L, "long double");
 162 #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
 163    test_spots(boost::math::concepts::real_concept(0.1), "real_concept");
 164 #endif
 165 #endif
 166 #endif
 167 #if !defined(TEST) || TEST == 2 || TEST == 3
 168    test_spots(boost::multiprecision::cpp_bin_float_quad(), "cpp_bin_float_quad");
 169 #endif
 170 }
 171
 172 #else
 173
 174 int main() { return 0; }
 175
 176 #endif