libs/math/test/test_legendre.hpp

   1 // Copyright John Maddock 2006.
   2 // Copyright Paul A. Bristow 2007, 2009
   3 //  Use, modification and distribution are subject to the
   4 //  Boost Software License, Version 1.0. (See accompanying file
   5 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6
   7 #ifdef _MSC_VER
   8 #  pragma warning (disable : 4756) // overflow in constant arithmetic
   9 #endif
  10
  11 #include <boost/math/concepts/real_concept.hpp>
  12 #define BOOST_TEST_MAIN
  13 #include <boost/test/unit_test.hpp>
  14 #include <boost/test/tools/floating_point_comparison.hpp>
  15 #include <boost/math/special_functions/math_fwd.hpp>
  16 #include <boost/math/special_functions/legendre.hpp>
  17 #include <boost/math/constants/constants.hpp>
  18 #include <boost/multiprecision/cpp_bin_float.hpp>
  19 #include <boost/array.hpp>
  20 #include "functor.hpp"
  21
  22 #include "handle_test_result.hpp"
  23 #include "table_type.hpp"
  24
  25 #ifndef SC_
  26 #define SC_(x) static_cast<typename table_type<T>::type>(BOOST_JOIN(x, L))
  27 #endif
  28
  29 template <class Real, class T>
  30 void do_test_legendre_p(const T& data, const char* type_name, const char* test_name)
  31 {
  32    typedef Real                   value_type;
  33
  34    typedef value_type (*pg)(int, value_type);
  35    pg funcp;
  36
  37 #if !(defined(ERROR_REPORTING_MODE) && !defined(LEGENDRE_P_FUNCTION_TO_TEST))
  38 #ifdef LEGENDRE_P_FUNCTION_TO_TEST
  39    funcp = LEGENDRE_P_FUNCTION_TO_TEST;
  40 #elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
  41    funcp = boost::math::legendre_p<value_type>;
  42 #else
  43    funcp = boost::math::legendre_p;
  44 #endif
  45
  46    boost::math::tools::test_result<value_type> result;
  47
  48    std::cout << "Testing " << test_name << " with type " << type_name
  49       << "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
  50
  51    //
  52    // test legendre_p against data:
  53    //
  54    result = boost::math::tools::test_hetero<Real>(
  55       data,
  56       bind_func_int1<Real>(funcp, 0, 1),
  57       extract_result<Real>(2));
  58    handle_test_result(result, data[result.worst()], result.worst(), type_name, "legendre_p", test_name);
  59 #endif
  60
  61    typedef value_type (*pg2)(unsigned, value_type);
  62 #if !(defined(ERROR_REPORTING_MODE) && !defined(LEGENDRE_Q_FUNCTION_TO_TEST))
  63 #ifdef LEGENDRE_Q_FUNCTION_TO_TEST
  64    pg2 funcp2 = LEGENDRE_Q_FUNCTION_TO_TEST;
  65 #elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
  66    pg2 funcp2 = boost::math::legendre_q<value_type>;
  67 #else
  68    pg2 funcp2 = boost::math::legendre_q;
  69 #endif
  70
  71    //
  72    // test legendre_q against data:
  73    //
  74    result = boost::math::tools::test_hetero<Real>(
  75       data,
  76       bind_func_int1<Real>(funcp2, 0, 1),
  77       extract_result<Real>(3));
  78    handle_test_result(result, data[result.worst()], result.worst(), type_name, "legendre_q", test_name);
  79
  80    std::cout << std::endl;
  81 #endif
  82 }
  83
  84 template <class Real, class T>
  85 void do_test_assoc_legendre_p(const T& data, const char* type_name, const char* test_name)
  86 {
  87 #if !(defined(ERROR_REPORTING_MODE) && !defined(LEGENDRE_PA_FUNCTION_TO_TEST))
  88    typedef Real                   value_type;
  89
  90    typedef value_type (*pg)(int, int, value_type);
  91 #ifdef LEGENDRE_PA_FUNCTION_TO_TEST
  92    pg funcp = LEGENDRE_PA_FUNCTION_TO_TEST;
  93 #elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
  94    pg funcp = boost::math::legendre_p<value_type>;
  95 #else
  96    pg funcp = boost::math::legendre_p;
  97 #endif
  98
  99    boost::math::tools::test_result<value_type> result;
 100
 101    std::cout << "Testing " << test_name << " with type " << type_name
 102       << "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
 103
 104    //
 105    // test legendre_p against data:
 106    //
 107    result = boost::math::tools::test_hetero<Real>(
 108       data,
 109       bind_func_int2<Real>(funcp, 0, 1, 2),
 110       extract_result<Real>(3));
 111    handle_test_result(result, data[result.worst()], result.worst(), type_name, "legendre_p (associated)", test_name);
 112    std::cout << std::endl;
 113 #endif
 114 }
 115
 116 template <class T>
 117 void test_legendre_p(T, const char* name)
 118 {
 119    //
 120    // The actual test data is rather verbose, so it's in a separate file
 121    //
 122    // The contents are as follows, each row of data contains
 123    // three items, input value a, input value b and erf(a, b):
 124    //
 125 #  include "legendre_p.ipp"
 126
 127    do_test_legendre_p<T>(legendre_p, name, "Legendre Polynomials: Small Values");
 128
 129 #  include "legendre_p_large.ipp"
 130
 131    do_test_legendre_p<T>(legendre_p_large, name, "Legendre Polynomials: Large Values");
 132
 133 #  include "assoc_legendre_p.ipp"
 134
 135    do_test_assoc_legendre_p<T>(assoc_legendre_p, name, "Associated Legendre Polynomials: Small Values");
 136
 137 }
 138
 139 template <class T>
 140 void test_spots(T, const char* t)
 141 {
 142    std::cout << "Testing basic sanity checks for type " << t << std::endl;
 143    //
 144    // basic sanity checks, tolerance is 100 epsilon:
 145    //
 146    T tolerance = boost::math::tools::epsilon<T>() * 100;
 147    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(1, static_cast<T>(0.5L)), static_cast<T>(0.5L), tolerance);
 148    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-1, static_cast<T>(0.5L)), static_cast<T>(1L), tolerance);
 149    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(4, static_cast<T>(0.5L)), static_cast<T>(-0.2890625000000000000000000000000000000000L), tolerance);
 150    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-4, static_cast<T>(0.5L)), static_cast<T>(-0.4375000000000000000000000000000000000000L), tolerance);
 151    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(7, static_cast<T>(0.5L)), static_cast<T>(0.2231445312500000000000000000000000000000L), tolerance);
 152    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-7, static_cast<T>(0.5L)), static_cast<T>(0.3232421875000000000000000000000000000000L), tolerance);
 153    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(40, static_cast<T>(0.5L)), static_cast<T>(-0.09542943523261546936538467572384923220258L), tolerance);
 154    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-40, static_cast<T>(0.5L)), static_cast<T>(-0.1316993126940266257030910566308990611306L), tolerance);
 155
 156    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(4, 2, static_cast<T>(0.5L)), static_cast<T>(4.218750000000000000000000000000000000000L), tolerance);
 157    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-4, 2, static_cast<T>(0.5L)), static_cast<T>(5.625000000000000000000000000000000000000L), tolerance);
 158    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(7, 5, static_cast<T>(0.5L)), static_cast<T>(-5696.789530152175143607977274672800795328L), tolerance);
 159    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-7, 4, static_cast<T>(0.5L)), static_cast<T>(465.1171875000000000000000000000000000000L), tolerance);
 160    if(std::numeric_limits<T>::max_exponent > std::numeric_limits<float>::max_exponent)
 161    {
 162       BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(40, 30, static_cast<T>(0.5L)), static_cast<T>(-7.855722083232252643913331343916012143461e45L), tolerance);
 163    }
 164    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-40, 20, static_cast<T>(0.5L)), static_cast<T>(4.966634149702370788037088925152355134665e30L), tolerance);
 165    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(4, 2, static_cast<T>(-0.5L)), static_cast<T>(4.218750000000000000000000000000000000000L), tolerance);
 166    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-4, 2, static_cast<T>(-0.5L)), static_cast<T>(-5.625000000000000000000000000000000000000L), tolerance);
 167    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(7, 5, static_cast<T>(-0.5L)), static_cast<T>(-5696.789530152175143607977274672800795328L), tolerance);
 168    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-7, 4, static_cast<T>(-0.5L)), static_cast<T>(465.1171875000000000000000000000000000000L), tolerance);
 169    if(std::numeric_limits<T>::max_exponent > std::numeric_limits<float>::max_exponent)
 170    {
 171       BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(40, 30, static_cast<T>(-0.5L)), static_cast<T>(-7.855722083232252643913331343916012143461e45L), tolerance);
 172    }
 173    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-40, 20, static_cast<T>(-0.5L)), static_cast<T>(-4.966634149702370788037088925152355134665e30L), tolerance);
 174    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(4, -2, static_cast<T>(0.5L)), static_cast<T>(0.01171875000000000000000000000000000000000L), tolerance);
 175    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-4, -2, static_cast<T>(0.5L)), static_cast<T>(0.04687500000000000000000000000000000000000L), tolerance);
 176    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(7, -5, static_cast<T>(0.5L)), static_cast<T>(0.00002378609812640364935569308025139290054701L), tolerance);
 177    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-7, -4, static_cast<T>(0.5L)), static_cast<T>(0.0002563476562500000000000000000000000000000L), tolerance);
 178    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(40, -30, static_cast<T>(0.5L)), static_cast<T>(-2.379819988646847616996471299410611801239e-48L), tolerance);
 179    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-40, -20, static_cast<T>(0.5L)), static_cast<T>(4.356454600748202401657099008867502679122e-33L), tolerance);
 180
 181    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_q(1, static_cast<T>(0.5L)), static_cast<T>(-0.7253469278329725771511886907693685738381L), tolerance);
 182    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_q(4, static_cast<T>(0.5L)), static_cast<T>(0.4401745259867706044988642951843745400835L), tolerance);
 183    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_q(7, static_cast<T>(0.5L)), static_cast<T>(-0.3439152932669753451878700644212067616780L), tolerance);
 184    BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_q(40, static_cast<T>(0.5L)), static_cast<T>(0.1493671665503550095010454949479907886011L), tolerance);
 185 }
 186
 187 template <class T>
 188 void test_legendre_p_prime()
 189 {
 190     T tolerance = 100*boost::math::tools::epsilon<T>();
 191     T x = -1;
 192     while (x <= 1)
 193     {
 194         // P_0'(x) = 0
 195         BOOST_CHECK_SMALL(::boost::math::legendre_p_prime<T>(0,  x), tolerance);
 196         // Reflection formula for P_{-1}(x) = P_{0}(x):
 197         BOOST_CHECK_SMALL(::boost::math::legendre_p_prime<T>(-1,  x), tolerance);
 198
 199         // P_1(x) = x, so P_1'(x) = 1:
 200         BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(1,  x), static_cast<T>(1), tolerance);
 201         BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-2,  x), static_cast<T>(1), tolerance);
 202
 203         // P_2(x) = 3x^2/2 + k => P_2'(x) = 3x
 204         BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(2,  x), 3*x, tolerance);
 205         BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-3,  x), 3*x, tolerance);
 206
 207         // P_3(x) = (5x^3 - 3x)/2 => P_3'(x) = (15x^2 - 3)/2:
 208         T xsq = x*x;
 209         BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(3,  x), (15*xsq - 3)/2, tolerance);
 210         BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-4,  x), (15*xsq -3)/2, tolerance);
 211
 212         // P_4(x) = (35x^4 - 30x^2 +3)/8 => P_4'(x) = (5x/2)*(7x^2 - 3)
 213         T expected = 5*x*(7*xsq - 3)/2;
 214         BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(4,  x), expected, tolerance);
 215         BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-5,  x), expected, tolerance);
 216
 217         // P_5(x) = (63x^5 - 70x^3 + 15x)/8 => P_5'(x) = (315*x^4 - 210*x^2 + 15)/8
 218         T x4 = xsq*xsq;
 219         expected = (315*x4 - 210*xsq + 15)/8;
 220         BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(5,  x), expected, tolerance);
 221         BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-6,  x), expected, tolerance);
 222
 223         // P_6(x) = (231x^6 -315*x^4 +105x^2 -5)/16 => P_6'(x) = (6*231*x^5 - 4*315*x^3 + 105x)/16
 224         expected = 21*x*(33*x4 - 30*xsq + 5)/8;
 225         BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(6,  x), expected, tolerance);
 226         BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-7,  x), expected, tolerance);
 227
 228         // Mathematica: D[LegendreP[7, x],x]
 229         T x6 = x4*xsq;
 230         expected = 7*(429*x6 -495*x4 + 135*xsq - 5)/16;
 231         BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(7,  x), expected, tolerance);
 232         BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-8,  x), expected, tolerance);
 233
 234         // Mathematica: D[LegendreP[8, x],x]
 235         // The naive polynomial evaluation algorithm is going to get worse from here out, so this will be enough.
 236         expected = 9*x*(715*x6 - 1001*x4 + 385*xsq - 35)/16;
 237         BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(8,  x), expected, tolerance);
 238         BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-9,  x), expected, tolerance);
 239
 240         x += static_cast<T>(1)/static_cast<T>(pow(T(2), T(4)));
 241     }
 242
 243     int n = 0;
 244     while (n < 5000)
 245     {
 246         T expected = n*(n+1)*boost::math::constants::half<T>();
 247         BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(n, (T) 1), expected, tolerance);
 248         BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-n - 1,  (T) 1), expected, tolerance);
 249         if (n & 1)
 250         {
 251             BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(n, (T) -1), expected, tolerance);
 252             BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-n - 1,  (T) -1), expected, tolerance);
 253         }
 254         else
 255         {
 256             BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(n, (T) -1), -expected, tolerance);
 257             BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-n - 1,  (T) -1), -expected, tolerance);
 258         }
 259         ++n;
 260     }
 261 }
 262
 263 template<class Real>
 264 void test_legendre_p_zeros()
 265 {
 266     std::cout << "Testing Legendre zeros on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 267     using std::sqrt;
 268     using std::abs;
 269     using boost::math::legendre_p_zeros;
 270     using boost::math::legendre_p;
 271     using boost::math::constants::third;
 272     Real tol = std::numeric_limits<Real>::epsilon();
 273
 274     // Check the trivial cases:
 275     std::vector<Real> zeros = legendre_p_zeros<Real>(1);
 276     BOOST_ASSERT(zeros.size() == 1);
 277     BOOST_CHECK_SMALL(zeros[0], tol);
 278
 279     zeros = legendre_p_zeros<Real>(2);
 280     BOOST_ASSERT(zeros.size() == 1);
 281     BOOST_CHECK_CLOSE_FRACTION(zeros[0], (Real) 1/ sqrt(static_cast<Real>(3)), tol);
 282
 283     zeros = legendre_p_zeros<Real>(3);
 284     BOOST_ASSERT(zeros.size() == 2);
 285     BOOST_CHECK_SMALL(zeros[0], tol);
 286     BOOST_CHECK_CLOSE_FRACTION(zeros[1], sqrt(static_cast<Real>(3)/static_cast<Real>(5)), tol);
 287
 288     zeros = legendre_p_zeros<Real>(4);
 289     BOOST_ASSERT(zeros.size() == 2);
 290     BOOST_CHECK_CLOSE_FRACTION(zeros[0], sqrt( (15-2*sqrt(static_cast<Real>(30)))/static_cast<Real>(35) ), tol);
 291     BOOST_CHECK_CLOSE_FRACTION(zeros[1], sqrt( (15+2*sqrt(static_cast<Real>(30)))/static_cast<Real>(35) ), tol);
 292
 293
 294     zeros = legendre_p_zeros<Real>(5);
 295     BOOST_ASSERT(zeros.size() == 3);
 296     BOOST_CHECK_SMALL(zeros[0], tol);
 297     BOOST_CHECK_CLOSE_FRACTION(zeros[1], third<Real>()*sqrt( (35 - 2*sqrt(static_cast<Real>(70)))/static_cast<Real>(7) ), 2*tol);
 298     BOOST_CHECK_CLOSE_FRACTION(zeros[2], third<Real>()*sqrt( (35 + 2*sqrt(static_cast<Real>(70)))/static_cast<Real>(7) ), 2*tol);
 299
 300     // Don't take the tolerances too seriously.
 301     // The other test shows that the zeros are estimated more accurately than the function!
 302     for (unsigned n = 6; n < 130; ++n)
 303     {
 304         zeros = legendre_p_zeros<Real>(n);
 305         if (n & 1)
 306         {
 307             BOOST_CHECK(zeros.size() == (n-1)/2 +1);
 308             BOOST_CHECK_SMALL(zeros[0], tol);
 309         }
 310         else
 311         {
 312             // Zero is not a zero of the odd Legendre polynomials
 313             BOOST_CHECK(zeros.size() == n/2);
 314             BOOST_CHECK(zeros[0] > 0);
 315             BOOST_CHECK_SMALL(legendre_p(n, zeros[0]), 550*tol);
 316         }
 317         Real previous_zero = zeros[0];
 318         for (unsigned k = 1; k < zeros.size(); ++k)
 319         {
 320             Real next_zero = zeros[k];
 321             BOOST_CHECK(next_zero > previous_zero);
 322
 323             std::string err = "Tolerance failed for (n, k) = (" + boost::lexical_cast<std::string>(n) + "," + boost::lexical_cast<std::string>(k) + ")\n";
 324             if (n < 40)
 325             {
 326                 BOOST_CHECK_MESSAGE( abs(legendre_p(n, next_zero)) < 100*tol,
 327                                      err);
 328             }
 329             else
 330             {
 331               BOOST_CHECK_MESSAGE( abs(legendre_p(n, next_zero)) < 1000*tol,
 332                                    err);
 333             }
 334             previous_zero = next_zero;
 335         }
 336         // The zeros of orthogonal polynomials are contained strictly in (a, b).
 337         BOOST_CHECK(previous_zero < 1);
 338     }
 339     return;
 340 }
 341
 342 int test_legendre_p_zeros_double_ulp(int min_x, int max_n)
 343 {
 344     std::cout << "Testing ULP distance for Legendre zeros.\n";
 345     using std::abs;
 346     using boost::math::legendre_p_zeros;
 347     using boost::math::float_distance;
 348     using boost::multiprecision::cpp_bin_float_quad;
 349
 350     double max_float_distance = 0;
 351     for (int n = min_x; n < max_n; ++n)
 352     {
 353         std::vector<double> double_zeros = legendre_p_zeros<double>(n);
 354         std::vector<cpp_bin_float_quad> quad_zeros   = legendre_p_zeros<cpp_bin_float_quad>(n);
 355         BOOST_ASSERT(quad_zeros.size() == double_zeros.size());
 356         for (int k = 0; k < (int)double_zeros.size(); ++k)
 357         {
 358             double d = abs(float_distance(double_zeros[k], quad_zeros[k].convert_to<double>()));
 359             if (d > max_float_distance)
 360             {
 361                 max_float_distance = d;
 362             }
 363         }
 364     }
 365
 366     return (int) max_float_distance;
 367 }