[ceph.git] / ceph / src / boost / libs / math / test / test_legendre.hpp

// Copyright John Maddock 2006.
// Copyright Paul A. Bristow 2007, 2009
//  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)

#ifdef _MSC_VER
#  pragma warning (disable : 4756) // overflow in constant arithmetic
#endif

#include <boost/math/concepts/real_concept.hpp>
#define BOOST_TEST_MAIN
#include <boost/test/unit_test.hpp>
#include <boost/test/tools/floating_point_comparison.hpp>
#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/special_functions/legendre.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
#include <boost/array.hpp>
#include "functor.hpp"

#include "handle_test_result.hpp"
#include "table_type.hpp"

#ifndef SC_
#define SC_(x) static_cast<typename table_type<T>::type>(BOOST_JOIN(x, L))
#endif

template <class Real, class T>
void do_test_legendre_p(const T& data, const char* type_name, const char* test_name)
{
   typedef Real                   value_type;

   typedef value_type (*pg)(int, value_type);
   pg funcp;

#if !(defined(ERROR_REPORTING_MODE) && !defined(LEGENDRE_P_FUNCTION_TO_TEST))
#ifdef LEGENDRE_P_FUNCTION_TO_TEST
   funcp = LEGENDRE_P_FUNCTION_TO_TEST;
#elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
   funcp = boost::math::legendre_p<value_type>;
#else
   funcp = boost::math::legendre_p;
#endif

   boost::math::tools::test_result<value_type> result;

   std::cout << "Testing " << test_name << " with type " << type_name
      << "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";

   //
   // test legendre_p against data:
   //
   result = boost::math::tools::test_hetero<Real>(
      data,
      bind_func_int1<Real>(funcp, 0, 1),
      extract_result<Real>(2));
   handle_test_result(result, data[result.worst()], result.worst(), type_name, "legendre_p", test_name);
#endif

   typedef value_type (*pg2)(unsigned, value_type);
#if !(defined(ERROR_REPORTING_MODE) && !defined(LEGENDRE_Q_FUNCTION_TO_TEST))
#ifdef LEGENDRE_Q_FUNCTION_TO_TEST
   pg2 funcp2 = LEGENDRE_Q_FUNCTION_TO_TEST;
#elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
   pg2 funcp2 = boost::math::legendre_q<value_type>;
#else
   pg2 funcp2 = boost::math::legendre_q;
#endif

   //
   // test legendre_q against data:
   //
   result = boost::math::tools::test_hetero<Real>(
      data,
      bind_func_int1<Real>(funcp2, 0, 1),
      extract_result<Real>(3));
   handle_test_result(result, data[result.worst()], result.worst(), type_name, "legendre_q", test_name);

   std::cout << std::endl;
#endif
}

template <class Real, class T>
void do_test_assoc_legendre_p(const T& data, const char* type_name, const char* test_name)
{
#if !(defined(ERROR_REPORTING_MODE) && !defined(LEGENDRE_PA_FUNCTION_TO_TEST))
   typedef Real                   value_type;

   typedef value_type (*pg)(int, int, value_type);
#ifdef LEGENDRE_PA_FUNCTION_TO_TEST
   pg funcp = LEGENDRE_PA_FUNCTION_TO_TEST;
#elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
   pg funcp = boost::math::legendre_p<value_type>;
#else
   pg funcp = boost::math::legendre_p;
#endif

   boost::math::tools::test_result<value_type> result;

   std::cout << "Testing " << test_name << " with type " << type_name
      << "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";

   //
   // test legendre_p against data:
   //
   result = boost::math::tools::test_hetero<Real>(
      data,
      bind_func_int2<Real>(funcp, 0, 1, 2),
      extract_result<Real>(3));
   handle_test_result(result, data[result.worst()], result.worst(), type_name, "legendre_p (associated)", test_name);
   std::cout << std::endl;
#endif
}

template <class T>
void test_legendre_p(T, const char* name)
{
   //
   // The actual test data is rather verbose, so it's in a separate file
   //
   // The contents are as follows, each row of data contains
   // three items, input value a, input value b and erf(a, b):
   //
#  include "legendre_p.ipp"

   do_test_legendre_p<T>(legendre_p, name, "Legendre Polynomials: Small Values");

#  include "legendre_p_large.ipp"

   do_test_legendre_p<T>(legendre_p_large, name, "Legendre Polynomials: Large Values");

#  include "assoc_legendre_p.ipp"

   do_test_assoc_legendre_p<T>(assoc_legendre_p, name, "Associated Legendre Polynomials: Small Values");

}

template <class T>
void test_spots(T, const char* t)
{
   std::cout << "Testing basic sanity checks for type " << t << std::endl;
   //
   // basic sanity checks, tolerance is 100 epsilon:
   //
   T tolerance = boost::math::tools::epsilon<T>() * 100;
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(1, static_cast<T>(0.5L)), static_cast<T>(0.5L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-1, static_cast<T>(0.5L)), static_cast<T>(1L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(4, static_cast<T>(0.5L)), static_cast<T>(-0.2890625000000000000000000000000000000000L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-4, static_cast<T>(0.5L)), static_cast<T>(-0.4375000000000000000000000000000000000000L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(7, static_cast<T>(0.5L)), static_cast<T>(0.2231445312500000000000000000000000000000L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-7, static_cast<T>(0.5L)), static_cast<T>(0.3232421875000000000000000000000000000000L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(40, static_cast<T>(0.5L)), static_cast<T>(-0.09542943523261546936538467572384923220258L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-40, static_cast<T>(0.5L)), static_cast<T>(-0.1316993126940266257030910566308990611306L), tolerance);

   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(4, 2, static_cast<T>(0.5L)), static_cast<T>(4.218750000000000000000000000000000000000L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-4, 2, static_cast<T>(0.5L)), static_cast<T>(5.625000000000000000000000000000000000000L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(7, 5, static_cast<T>(0.5L)), static_cast<T>(-5696.789530152175143607977274672800795328L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-7, 4, static_cast<T>(0.5L)), static_cast<T>(465.1171875000000000000000000000000000000L), tolerance);
   if(std::numeric_limits<T>::max_exponent > std::numeric_limits<float>::max_exponent)
   {
      BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(40, 30, static_cast<T>(0.5L)), static_cast<T>(-7.855722083232252643913331343916012143461e45L), tolerance);
   }
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-40, 20, static_cast<T>(0.5L)), static_cast<T>(4.966634149702370788037088925152355134665e30L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(4, 2, static_cast<T>(-0.5L)), static_cast<T>(4.218750000000000000000000000000000000000L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-4, 2, static_cast<T>(-0.5L)), static_cast<T>(-5.625000000000000000000000000000000000000L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(7, 5, static_cast<T>(-0.5L)), static_cast<T>(-5696.789530152175143607977274672800795328L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-7, 4, static_cast<T>(-0.5L)), static_cast<T>(465.1171875000000000000000000000000000000L), tolerance);
   if(std::numeric_limits<T>::max_exponent > std::numeric_limits<float>::max_exponent)
   {
      BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(40, 30, static_cast<T>(-0.5L)), static_cast<T>(-7.855722083232252643913331343916012143461e45L), tolerance);
   }
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-40, 20, static_cast<T>(-0.5L)), static_cast<T>(-4.966634149702370788037088925152355134665e30L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(4, -2, static_cast<T>(0.5L)), static_cast<T>(0.01171875000000000000000000000000000000000L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-4, -2, static_cast<T>(0.5L)), static_cast<T>(0.04687500000000000000000000000000000000000L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(7, -5, static_cast<T>(0.5L)), static_cast<T>(0.00002378609812640364935569308025139290054701L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-7, -4, static_cast<T>(0.5L)), static_cast<T>(0.0002563476562500000000000000000000000000000L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(40, -30, static_cast<T>(0.5L)), static_cast<T>(-2.379819988646847616996471299410611801239e-48L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(-40, -20, static_cast<T>(0.5L)), static_cast<T>(4.356454600748202401657099008867502679122e-33L), tolerance);

   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_q(1, static_cast<T>(0.5L)), static_cast<T>(-0.7253469278329725771511886907693685738381L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_q(4, static_cast<T>(0.5L)), static_cast<T>(0.4401745259867706044988642951843745400835L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_q(7, static_cast<T>(0.5L)), static_cast<T>(-0.3439152932669753451878700644212067616780L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_q(40, static_cast<T>(0.5L)), static_cast<T>(0.1493671665503550095010454949479907886011L), tolerance);
   //
   // m = l-1, see https://github.com/boostorg/math/issues/453:
   //
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(0, -1, static_cast<T>(0.5)), static_cast<T>(0.5773502691896257645091487805019574556476017512701268760186023264L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(0, -1, static_cast<T>(1)), static_cast<T>(0), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(0, -1, static_cast<T>(0)), static_cast<T>(1), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(1, 0, static_cast<T>(0.5)), static_cast<T>(0.5), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(1, 0, static_cast<T>(1)), static_cast<T>(1), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(1, 0, static_cast<T>(0)), static_cast<T>(0), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(2, 1, static_cast<T>(0.5)), static_cast<T>(-1.2990381056766579701455847561294042752071039403577854710418552L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(2, 1, static_cast<T>(1)), static_cast<T>(0), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(2, 1, static_cast<T>(0)), static_cast<T>(0), tolerance);

   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(2, -2, static_cast<T>(0.5)), static_cast<T>(0.09375L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(2, -2, static_cast<T>(1)), static_cast<T>(0), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(2, -2, static_cast<T>(0)), static_cast<T>(0.125), tolerance);
}

template <class T>
void test_legendre_p_prime()
{
    T tolerance = 100*boost::math::tools::epsilon<T>();
    T x = -1;
    while (x <= 1)
    {
        // P_0'(x) = 0
        BOOST_CHECK_SMALL(::boost::math::legendre_p_prime<T>(0,  x), tolerance);
        // Reflection formula for P_{-1}(x) = P_{0}(x):
        BOOST_CHECK_SMALL(::boost::math::legendre_p_prime<T>(-1,  x), tolerance);

        // P_1(x) = x, so P_1'(x) = 1:
        BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(1,  x), static_cast<T>(1), tolerance);
        BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-2,  x), static_cast<T>(1), tolerance);

        // P_2(x) = 3x^2/2 + k => P_2'(x) = 3x
        BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(2,  x), 3*x, tolerance);
        BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-3,  x), 3*x, tolerance);

        // P_3(x) = (5x^3 - 3x)/2 => P_3'(x) = (15x^2 - 3)/2:
        T xsq = x*x;
        BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(3,  x), (15*xsq - 3)/2, tolerance);
        BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-4,  x), (15*xsq -3)/2, tolerance);

        // P_4(x) = (35x^4 - 30x^2 +3)/8 => P_4'(x) = (5x/2)*(7x^2 - 3)
        T expected = 5*x*(7*xsq - 3)/2;
        BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(4,  x), expected, tolerance);
        BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-5,  x), expected, tolerance);

        // P_5(x) = (63x^5 - 70x^3 + 15x)/8 => P_5'(x) = (315*x^4 - 210*x^2 + 15)/8
        T x4 = xsq*xsq;
        expected = (315*x4 - 210*xsq + 15)/8;
        BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(5,  x), expected, tolerance);
        BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-6,  x), expected, tolerance);

        // 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
        expected = 21*x*(33*x4 - 30*xsq + 5)/8;
        BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(6,  x), expected, tolerance);
        BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-7,  x), expected, tolerance);

        // Mathematica: D[LegendreP[7, x],x]
        T x6 = x4*xsq;
        expected = 7*(429*x6 -495*x4 + 135*xsq - 5)/16;
        BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(7,  x), expected, tolerance);
        BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-8,  x), expected, tolerance);

        // Mathematica: D[LegendreP[8, x],x]
        // The naive polynomial evaluation algorithm is going to get worse from here out, so this will be enough.
        expected = 9*x*(715*x6 - 1001*x4 + 385*xsq - 35)/16;
        BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(8,  x), expected, tolerance);
        BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-9,  x), expected, tolerance);

        x += static_cast<T>(1)/static_cast<T>(pow(T(2), T(4)));
    }

    int n = 0;
    while (n < 5000)
    {
        T expected = n*(n+1)*boost::math::constants::half<T>();
        BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(n, (T) 1), expected, tolerance);
        BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-n - 1,  (T) 1), expected, tolerance);
        if (n & 1)
        {
            BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(n, (T) -1), expected, tolerance);
            BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-n - 1,  (T) -1), expected, tolerance);
        }
        else
        {
            BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(n, (T) -1), -expected, tolerance);
            BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-n - 1,  (T) -1), -expected, tolerance);
        }
        ++n;
    }
}

template<class Real>
void test_legendre_p_zeros()
{
    std::cout << "Testing Legendre zeros on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
    using std::sqrt;
    using std::abs;
    using boost::math::legendre_p_zeros;
    using boost::math::legendre_p;
    using boost::math::constants::third;
    Real tol = std::numeric_limits<Real>::epsilon();

    // Check the trivial cases:
    std::vector<Real> zeros = legendre_p_zeros<Real>(1);
    BOOST_MATH_ASSERT(zeros.size() == 1);
    BOOST_CHECK_SMALL(zeros[0], tol);

    zeros = legendre_p_zeros<Real>(2);
    BOOST_MATH_ASSERT(zeros.size() == 1);
    BOOST_CHECK_CLOSE_FRACTION(zeros[0], (Real) 1/ sqrt(static_cast<Real>(3)), tol);

    zeros = legendre_p_zeros<Real>(3);
    BOOST_MATH_ASSERT(zeros.size() == 2);
    BOOST_CHECK_SMALL(zeros[0], tol);
    BOOST_CHECK_CLOSE_FRACTION(zeros[1], sqrt(static_cast<Real>(3)/static_cast<Real>(5)), tol);

    zeros = legendre_p_zeros<Real>(4);
    BOOST_MATH_ASSERT(zeros.size() == 2);
    BOOST_CHECK_CLOSE_FRACTION(zeros[0], sqrt( (15-2*sqrt(static_cast<Real>(30)))/static_cast<Real>(35) ), tol);
    BOOST_CHECK_CLOSE_FRACTION(zeros[1], sqrt( (15+2*sqrt(static_cast<Real>(30)))/static_cast<Real>(35) ), tol);


    zeros = legendre_p_zeros<Real>(5);
    BOOST_MATH_ASSERT(zeros.size() == 3);
    BOOST_CHECK_SMALL(zeros[0], tol);
    BOOST_CHECK_CLOSE_FRACTION(zeros[1], third<Real>()*sqrt( (35 - 2*sqrt(static_cast<Real>(70)))/static_cast<Real>(7) ), 2*tol);
    BOOST_CHECK_CLOSE_FRACTION(zeros[2], third<Real>()*sqrt( (35 + 2*sqrt(static_cast<Real>(70)))/static_cast<Real>(7) ), 2*tol);

    // Don't take the tolerances too seriously.
    // The other test shows that the zeros are estimated more accurately than the function!
    for (unsigned n = 6; n < 130; ++n)
    {
        zeros = legendre_p_zeros<Real>(n);
        if (n & 1)
        {
            BOOST_CHECK(zeros.size() == (n-1)/2 +1);
            BOOST_CHECK_SMALL(zeros[0], tol);
        }
        else
        {
            // Zero is not a zero of the odd Legendre polynomials
            BOOST_CHECK(zeros.size() == n/2);
            BOOST_CHECK(zeros[0] > 0);
            BOOST_CHECK_SMALL(legendre_p(n, zeros[0]), 550*tol);
        }
        Real previous_zero = zeros[0];
        for (unsigned k = 1; k < zeros.size(); ++k)
        {
            Real next_zero = zeros[k];
            BOOST_CHECK(next_zero > previous_zero);

            std::string err = "Tolerance failed for (n, k) = (" + std::to_string(n) + "," + std::to_string(k) + ")\n";
            if (n < 40)
            {
                BOOST_CHECK_MESSAGE( abs(legendre_p(n, next_zero)) < 100*tol,
                                     err);
            }
            else
            {
              BOOST_CHECK_MESSAGE( abs(legendre_p(n, next_zero)) < 1000*tol,
                                   err);
            }
            previous_zero = next_zero;
        }
        // The zeros of orthogonal polynomials are contained strictly in (a, b).
        BOOST_CHECK(previous_zero < 1);
    }
    return;
}

int test_legendre_p_zeros_double_ulp(int min_x, int max_n)
{
    std::cout << "Testing ULP distance for Legendre zeros.\n";
    using std::abs;
    using boost::math::legendre_p_zeros;
    using boost::math::float_distance;
    using boost::multiprecision::cpp_bin_float_quad;

    double max_float_distance = 0;
    for (int n = min_x; n < max_n; ++n)
    {
        std::vector<double> double_zeros = legendre_p_zeros<double>(n);
        std::vector<cpp_bin_float_quad> quad_zeros   = legendre_p_zeros<cpp_bin_float_quad>(n);
        BOOST_MATH_ASSERT(quad_zeros.size() == double_zeros.size());
        for (int k = 0; k < (int)double_zeros.size(); ++k)
        {
            double d = abs(float_distance(double_zeros[k], quad_zeros[k].convert_to<double>()));
            if (d > max_float_distance)
            {
                max_float_distance = d;
            }
        }
    }

    return (int) max_float_distance;
}
Commit	Line	Data
7c673cae FG	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>
92f5a8d4	14	#include <boost/test/tools/floating_point_comparison.hpp>
7c673cae	15	#include <boost/math/special_functions/math_fwd.hpp>
b32b8144	16	#include <boost/math/special_functions/legendre.hpp>
7c673cae	17	#include <boost/math/constants/constants.hpp>
b32b8144	18	#include <boost/multiprecision/cpp_bin_float.hpp>
7c673cae FG	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);
20effc67 TL	185	//
	186	// m = l-1, see https://github.com/boostorg/math/issues/453:
	187	//
	188	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(0, -1, static_cast<T>(0.5)), static_cast<T>(0.5773502691896257645091487805019574556476017512701268760186023264L), tolerance);
	189	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(0, -1, static_cast<T>(1)), static_cast<T>(0), tolerance);
	190	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(0, -1, static_cast<T>(0)), static_cast<T>(1), tolerance);
	191	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(1, 0, static_cast<T>(0.5)), static_cast<T>(0.5), tolerance);
	192	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(1, 0, static_cast<T>(1)), static_cast<T>(1), tolerance);
	193	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(1, 0, static_cast<T>(0)), static_cast<T>(0), tolerance);
	194	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(2, 1, static_cast<T>(0.5)), static_cast<T>(-1.2990381056766579701455847561294042752071039403577854710418552L), tolerance);
	195	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(2, 1, static_cast<T>(1)), static_cast<T>(0), tolerance);
	196	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(2, 1, static_cast<T>(0)), static_cast<T>(0), tolerance);
	197
	198	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(2, -2, static_cast<T>(0.5)), static_cast<T>(0.09375L), tolerance);
	199	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(2, -2, static_cast<T>(1)), static_cast<T>(0), tolerance);
	200	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p(2, -2, static_cast<T>(0)), static_cast<T>(0.125), tolerance);
7c673cae FG	201	}
7c673cae FG	202
b32b8144 FG	203	template <class T>
	204	void test_legendre_p_prime()
	205	{
	206	T tolerance = 100*boost::math::tools::epsilon<T>();
	207	T x = -1;
	208	while (x <= 1)
	209	{
	210	// P_0'(x) = 0
	211	BOOST_CHECK_SMALL(::boost::math::legendre_p_prime<T>(0, x), tolerance);
	212	// Reflection formula for P_{-1}(x) = P_{0}(x):
	213	BOOST_CHECK_SMALL(::boost::math::legendre_p_prime<T>(-1, x), tolerance);
	214
	215	// P_1(x) = x, so P_1'(x) = 1:
	216	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(1, x), static_cast<T>(1), tolerance);
	217	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-2, x), static_cast<T>(1), tolerance);
	218
	219	// P_2(x) = 3x^2/2 + k => P_2'(x) = 3x
	220	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(2, x), 3*x, tolerance);
	221	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-3, x), 3*x, tolerance);
	222
	223	// P_3(x) = (5x^3 - 3x)/2 => P_3'(x) = (15x^2 - 3)/2:
	224	T xsq = x*x;
	225	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(3, x), (15*xsq - 3)/2, tolerance);
	226	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-4, x), (15*xsq -3)/2, tolerance);
	227
	228	// P_4(x) = (35x^4 - 30x^2 +3)/8 => P_4'(x) = (5x/2)*(7x^2 - 3)
	229	T expected = 5x(7*xsq - 3)/2;
	230	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(4, x), expected, tolerance);
	231	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-5, x), expected, tolerance);
	232
	233	// P_5(x) = (63x^5 - 70x^3 + 15x)/8 => P_5'(x) = (315x^4 - 210x^2 + 15)/8
	234	T x4 = xsq*xsq;
	235	expected = (315x4 - 210xsq + 15)/8;
	236	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(5, x), expected, tolerance);
	237	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-6, x), expected, tolerance);
	238
	239	// P_6(x) = (231x^6 -315x^4 +105x^2 -5)/16 => P_6'(x) = (6231x^5 - 4315*x^3 + 105x)/16
	240	expected = 21x(33x4 - 30xsq + 5)/8;
	241	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(6, x), expected, tolerance);
	242	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-7, x), expected, tolerance);
	243
	244	// Mathematica: D[LegendreP[7, x],x]
	245	T x6 = x4*xsq;
	246	expected = 7(429x6 -495x4 + 135xsq - 5)/16;
	247	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(7, x), expected, tolerance);
	248	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-8, x), expected, tolerance);
	249
	250	// Mathematica: D[LegendreP[8, x],x]
	251	// The naive polynomial evaluation algorithm is going to get worse from here out, so this will be enough.
	252	expected = 9x(715x6 - 1001x4 + 385*xsq - 35)/16;
	253	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(8, x), expected, tolerance);
	254	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-9, x), expected, tolerance);
	255
	256	x += static_cast<T>(1)/static_cast<T>(pow(T(2), T(4)));
	257	}
	258
	259	int n = 0;
	260	while (n < 5000)
	261	{
	262	T expected = n(n+1)boost::math::constants::half<T>();
	263	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(n, (T) 1), expected, tolerance);
	264	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-n - 1, (T) 1), expected, tolerance);
	265	if (n & 1)
	266	{
267	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(n, (T) -1), expected, tolerance);
268	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-n - 1, (T) -1), expected, tolerance);
269	}
270	else
271	{
272	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(n, (T) -1), -expected, tolerance);
273	BOOST_CHECK_CLOSE_FRACTION(::boost::math::legendre_p_prime<T>(-n - 1, (T) -1), -expected, tolerance);
274	}
275	++n;
276	}
277	}
278
279	template<class Real>
280	void test_legendre_p_zeros()
281	{
282	std::cout << "Testing Legendre zeros on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
283	using std::sqrt;
284	using std::abs;
285	using boost::math::legendre_p_zeros;
286	using boost::math::legendre_p;
287	using boost::math::constants::third;
288	Real tol = std::numeric_limits<Real>::epsilon();
289
290	// Check the trivial cases:
291	std::vector<Real> zeros = legendre_p_zeros<Real>(1);
1e59de90	292	BOOST_MATH_ASSERT(zeros.size() == 1);
b32b8144 FG	293	BOOST_CHECK_SMALL(zeros[0], tol);
	294
	295	zeros = legendre_p_zeros<Real>(2);
1e59de90	296	BOOST_MATH_ASSERT(zeros.size() == 1);
b32b8144 FG	297	BOOST_CHECK_CLOSE_FRACTION(zeros[0], (Real) 1/ sqrt(static_cast<Real>(3)), tol);
	298
	299	zeros = legendre_p_zeros<Real>(3);
1e59de90	300	BOOST_MATH_ASSERT(zeros.size() == 2);
b32b8144 FG	301	BOOST_CHECK_SMALL(zeros[0], tol);
	302	BOOST_CHECK_CLOSE_FRACTION(zeros[1], sqrt(static_cast<Real>(3)/static_cast<Real>(5)), tol);
	303
	304	zeros = legendre_p_zeros<Real>(4);
1e59de90	305	BOOST_MATH_ASSERT(zeros.size() == 2);
b32b8144 FG	306	BOOST_CHECK_CLOSE_FRACTION(zeros[0], sqrt( (15-2*sqrt(static_cast<Real>(30)))/static_cast<Real>(35) ), tol);
	307	BOOST_CHECK_CLOSE_FRACTION(zeros[1], sqrt( (15+2*sqrt(static_cast<Real>(30)))/static_cast<Real>(35) ), tol);
	308
	309
	310	zeros = legendre_p_zeros<Real>(5);
1e59de90	311	BOOST_MATH_ASSERT(zeros.size() == 3);
b32b8144 FG	312	BOOST_CHECK_SMALL(zeros[0], tol);
	313	BOOST_CHECK_CLOSE_FRACTION(zeros[1], third<Real>()sqrt( (35 - 2sqrt(static_cast<Real>(70)))/static_cast<Real>(7) ), 2*tol);
	314	BOOST_CHECK_CLOSE_FRACTION(zeros[2], third<Real>()sqrt( (35 + 2sqrt(static_cast<Real>(70)))/static_cast<Real>(7) ), 2*tol);
	315
	316	// Don't take the tolerances too seriously.
	317	// The other test shows that the zeros are estimated more accurately than the function!
92f5a8d4	318	for (unsigned n = 6; n < 130; ++n)
b32b8144 FG	319	{
	320	zeros = legendre_p_zeros<Real>(n);
	321	if (n & 1)
	322	{
	323	BOOST_CHECK(zeros.size() == (n-1)/2 +1);
	324	BOOST_CHECK_SMALL(zeros[0], tol);
	325	}
	326	else
	327	{
	328	// Zero is not a zero of the odd Legendre polynomials
	329	BOOST_CHECK(zeros.size() == n/2);
	330	BOOST_CHECK(zeros[0] > 0);
	331	BOOST_CHECK_SMALL(legendre_p(n, zeros[0]), 550*tol);
	332	}
	333	Real previous_zero = zeros[0];
92f5a8d4	334	for (unsigned k = 1; k < zeros.size(); ++k)
b32b8144 FG	335	{
	336	Real next_zero = zeros[k];
	337	BOOST_CHECK(next_zero > previous_zero);
	338
1e59de90	339	std::string err = "Tolerance failed for (n, k) = (" + std::to_string(n) + "," + std::to_string(k) + ")\n";
b32b8144 FG	340	if (n < 40)
	341	{
	342	BOOST_CHECK_MESSAGE( abs(legendre_p(n, next_zero)) < 100*tol,
	343	err);
	344	}
	345	else
	346	{
	347	BOOST_CHECK_MESSAGE( abs(legendre_p(n, next_zero)) < 1000*tol,
	348	err);
	349	}
	350	previous_zero = next_zero;
	351	}
	352	// The zeros of orthogonal polynomials are contained strictly in (a, b).
	353	BOOST_CHECK(previous_zero < 1);
	354	}
	355	return;
	356	}
	357
	358	int test_legendre_p_zeros_double_ulp(int min_x, int max_n)
	359	{
	360	std::cout << "Testing ULP distance for Legendre zeros.\n";
	361	using std::abs;
	362	using boost::math::legendre_p_zeros;
	363	using boost::math::float_distance;
	364	using boost::multiprecision::cpp_bin_float_quad;
	365
	366	double max_float_distance = 0;
	367	for (int n = min_x; n < max_n; ++n)
	368	{
	369	std::vector<double> double_zeros = legendre_p_zeros<double>(n);
	370	std::vector<cpp_bin_float_quad> quad_zeros = legendre_p_zeros<cpp_bin_float_quad>(n);
1e59de90	371	BOOST_MATH_ASSERT(quad_zeros.size() == double_zeros.size());
b32b8144 FG	372	for (int k = 0; k < (int)double_zeros.size(); ++k)
	373	{
	374	double d = abs(float_distance(double_zeros[k], quad_zeros[k].convert_to<double>()));
	375	if (d > max_float_distance)
	376	{
	377	max_float_distance = d;
	378	}
	379	}
	380	}
	381
	382	return (int) max_float_distance;
	383	}