libs/math/test/gegenbauer_test.cpp

   1 /*
   2  * Copyright Nick Thompson, 2019
   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
   8 #include "math_unit_test.hpp"
   9 #include <numeric>
  10 #include <utility>
  11 #include <random>
  12 #include <cmath>
  13 #include <boost/core/demangle.hpp>
  14 #include <boost/math/special_functions/gegenbauer.hpp>
  15 #ifdef BOOST_HAS_FLOAT128
  16 #include <boost/multiprecision/float128.hpp>
  17 using boost::multiprecision::float128;
  18 #endif
  19
  20 using std::abs;
  21 using boost::math::gegenbauer;
  22 using boost::math::gegenbauer_derivative;
  23
  24 template<class Real>
  25 void test_parity()
  26 {
  27     std::mt19937 gen(323723);
  28     std::uniform_real_distribution<Real> xdis(-1, +1);
  29     std::uniform_real_distribution<Real> lambdadis(-0.5, 1);
  30
  31     for(unsigned n = 0; n < 50; ++n) {
  32         unsigned calls = 50;
  33         unsigned i = 0;
  34         while(i++ < calls) {
  35             Real x = xdis(gen);
  36             Real lambda = lambdadis(gen);
  37             if (n & 1) {
  38                 CHECK_ULP_CLOSE(gegenbauer(n, lambda, -x), -gegenbauer(n, lambda, x), 0);
  39             }
  40             else {
  41                 CHECK_ULP_CLOSE(gegenbauer(n, lambda, -x), gegenbauer(n, lambda, x), 0);
  42             }
  43         }
  44     }
  45 }
  46
  47 template<class Real>
  48 void test_quadratic()
  49 {
  50     Real lambda = 1/Real(4);
  51     auto c2 = [&](Real x) { return -lambda + 2*lambda*(1+lambda)*x*x; };
  52
  53     Real x = -1;
  54     Real h = 1/Real(256);
  55     while (x < 1) {
  56         Real expected = c2(x);
  57         Real computed = gegenbauer(2, lambda, x);
  58         CHECK_ULP_CLOSE(expected, computed, 0);
  59         x += h;
  60     }
  61 }
  62
  63 template<class Real>
  64 void test_cubic()
  65 {
  66     Real lambda = 1/Real(4);
  67     auto c3 = [&](Real x) { return lambda*(1+lambda)*x*(-2 + 4*(2+lambda)*x*x/3); };
  68
  69     Real x = -1;
  70     Real h = 1/Real(256);
  71     while (x < 1) {
  72         Real expected = c3(x);
  73         Real computed = gegenbauer(3, lambda, x);
  74         CHECK_ULP_CLOSE(expected, computed, 4);
  75         x += h;
  76     }
  77 }
  78
  79 template<class Real>
  80 void test_derivative()
  81 {
  82     Real lambda = 0.5;
  83     auto c3_prime = [&](Real x) { return 2*lambda*(lambda+1)*(-1 + 2*(lambda+2)*x*x); };
  84     auto c3_double_prime = [&](Real x) { return 8*lambda*(lambda+1)*(lambda+2)*x; };
  85     Real x = -1;
  86     Real h = 1/Real(256);
  87     while (x < 1) {
  88         Real expected = c3_prime(x);
  89         Real computed = gegenbauer_derivative(3, lambda, x, 1);
  90         CHECK_ULP_CLOSE(expected, computed, 1);
  91
  92         expected = c3_double_prime(x);
  93         computed = gegenbauer_derivative(3, lambda, x, 2);
  94         CHECK_ULP_CLOSE(expected, computed, 1);
  95
  96         x += h;
  97     }
  98
  99 }
 100
 101
 102
 103 int main()
 104 {
 105     test_parity<float>();
 106     test_parity<double>();
 107     test_parity<long double>();
 108
 109     test_quadratic<float>();
 110     test_quadratic<double>();
 111     test_quadratic<long double>();
 112
 113     test_cubic<double>();
 114     test_cubic<long double>();
 115
 116     test_derivative<float>();
 117     test_derivative<double>();
 118     test_derivative<long double>();
 119
 120 #ifdef BOOST_HAS_FLOAT128
 121     test_quadratic<boost::multiprecision::float128>();
 122     test_cubic<boost::multiprecision::float128>();
 123 #endif
 124
 125     return boost::math::test::report_errors();
 126 }