[ceph.git] / ceph / src / boost / libs / math / test / test_carlson.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)

#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/constants/constants.hpp>
#include <boost/array.hpp>
#include <boost/random.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, typename T>
void do_test_ellint_rf(T& data, const char* type_name, const char* test)
{
#if !(defined(ERROR_REPORTING_MODE) && !defined(ELLINT_RF_FUNCTION_TO_TEST))
   typedef Real                   value_type;

   std::cout << "Testing: " << test << std::endl;

#ifdef ELLINT_RF_FUNCTION_TO_TEST
   value_type(*fp)(value_type, value_type, value_type) = ELLINT_RF_FUNCTION_TO_TEST;
#elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
    value_type (*fp)(value_type, value_type, value_type) = boost::math::ellint_rf<value_type, value_type, value_type>;
#else
    value_type (*fp)(value_type, value_type, value_type) = boost::math::ellint_rf;
#endif
    boost::math::tools::test_result<value_type> result;
 
    result = boost::math::tools::test_hetero<Real>(
      data, 
      bind_func<Real>(fp, 0, 1, 2),
      extract_result<Real>(3));
   handle_test_result(result, data[result.worst()], result.worst(), 
      type_name, "ellint_rf", test);

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

template <class Real, typename T>
void do_test_ellint_rc(T& data, const char* type_name, const char* test)
{
#if !(defined(ERROR_REPORTING_MODE) && !defined(ELLINT_RC_FUNCTION_TO_TEST))
   typedef Real                   value_type;

   std::cout << "Testing: " << test << std::endl;

#ifdef ELLINT_RC_FUNCTION_TO_TEST
   value_type(*fp)(value_type, value_type) = ELLINT_RC_FUNCTION_TO_TEST;
#elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
    value_type (*fp)(value_type, value_type) = boost::math::ellint_rc<value_type, value_type>;
#else
    value_type (*fp)(value_type, value_type) = boost::math::ellint_rc;
#endif
    boost::math::tools::test_result<value_type> result;
 
    result = boost::math::tools::test_hetero<Real>(
      data, 
      bind_func<Real>(fp, 0, 1),
      extract_result<Real>(2));
      handle_test_result(result, data[result.worst()], result.worst(), 
      type_name, "ellint_rc", test);

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

template <class Real, typename T>
void do_test_ellint_rj(T& data, const char* type_name, const char* test)
{
#if !(defined(ERROR_REPORTING_MODE) && !defined(ELLINT_RJ_FUNCTION_TO_TEST))
   typedef Real                   value_type;

   std::cout << "Testing: " << test << std::endl;

#ifdef ELLINT_RJ_FUNCTION_TO_TEST
   value_type(*fp)(value_type, value_type, value_type, value_type) = ELLINT_RJ_FUNCTION_TO_TEST;
#elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
    value_type (*fp)(value_type, value_type, value_type, value_type) = boost::math::ellint_rj<value_type, value_type, value_type, value_type>;
#else
    value_type (*fp)(value_type, value_type, value_type, value_type) = boost::math::ellint_rj;
#endif
    boost::math::tools::test_result<value_type> result;
 
    result = boost::math::tools::test_hetero<Real>(
      data, 
      bind_func<Real>(fp, 0, 1, 2, 3),
      extract_result<Real>(4));
      handle_test_result(result, data[result.worst()], result.worst(), 
      type_name, "ellint_rj", test);

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

template <class Real, typename T>
void do_test_ellint_rd(T& data, const char* type_name, const char* test)
{
#if !(defined(ERROR_REPORTING_MODE) && !defined(ELLINT_RD_FUNCTION_TO_TEST))
   typedef Real                   value_type;

   std::cout << "Testing: " << test << std::endl;

#ifdef ELLINT_RD_FUNCTION_TO_TEST
   value_type(*fp)(value_type, value_type, value_type) = ELLINT_RD_FUNCTION_TO_TEST;
#elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
    value_type (*fp)(value_type, value_type, value_type) = boost::math::ellint_rd<value_type, value_type, value_type>;
#else
    value_type (*fp)(value_type, value_type, value_type) = boost::math::ellint_rd;
#endif
    boost::math::tools::test_result<value_type> result;
 
    result = boost::math::tools::test_hetero<Real>(
      data, 
      bind_func<Real>(fp, 0, 1, 2),
      extract_result<Real>(3));
    handle_test_result(result, data[result.worst()], result.worst(), 
      type_name, "ellint_rd", test);

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

template <class Real, typename T>
void do_test_ellint_rg(T& data, const char* type_name, const char* test)
{
#if !(defined(ERROR_REPORTING_MODE) && !defined(ELLINT_RD_FUNCTION_TO_TEST))
   typedef Real                   value_type;

   std::cout << "Testing: " << test << std::endl;

#ifdef ELLINT_RG_FUNCTION_TO_TEST
   value_type(*fp)(value_type, value_type, value_type) = ELLINT_RG_FUNCTION_TO_TEST;
#elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
   value_type(*fp)(value_type, value_type, value_type) = boost::math::ellint_rg<value_type, value_type, value_type>;
#else
   value_type(*fp)(value_type, value_type, value_type) = boost::math::ellint_rg;
#endif
   boost::math::tools::test_result<value_type> result;

   result = boost::math::tools::test_hetero<Real>(
      data,
      bind_func<Real>(fp, 0, 1, 2),
      extract_result<Real>(3));
   handle_test_result(result, data[result.worst()], result.worst(),
      type_name, "ellint_rg", test);

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

#if !defined(TEST1) && !defined(TEST2) && !defined(TEST3) && !defined(TEST4)
#define TEST1
#define TEST2
#define TEST3
#define TEST4
#endif

#ifdef TEST1

template <typename T>
void t1(T, const char* type_name)
{
#include "ellint_rf_data.ipp"

   do_test_ellint_rf<T>(ellint_rf_data, type_name, "RF: Random data");
}

template <typename T>
void t2(T, const char* type_name)
{
#include "ellint_rf_xxx.ipp"

   do_test_ellint_rf<T>(ellint_rf_xxx, type_name, "RF: x = y = z");
}

template <typename T>
void t3(T, const char* type_name)
{
#include "ellint_rf_xyy.ipp"

   do_test_ellint_rf<T>(ellint_rf_xyy, type_name, "RF: x = y or y = z or x = z");
}

template <typename T>
void t4(T, const char* type_name)
{
#include "ellint_rf_0yy.ipp"

   do_test_ellint_rf<T>(ellint_rf_0yy, type_name, "RF: x = 0, y = z");
}

template <typename T>
void t5(T, const char* type_name)
{
#include "ellint_rf_xy0.ipp"

   do_test_ellint_rf<T>(ellint_rf_xy0, type_name, "RF: z = 0");
}

#endif
#ifdef TEST2

template <typename T>
void t6(T, const char* type_name)
{
#include "ellint_rc_data.ipp"

   do_test_ellint_rc<T>(ellint_rc_data, type_name, "RC: Random data");
}

template <typename T>
void t7(T, const char* type_name)
{
#include "ellint_rj_data.ipp"

   do_test_ellint_rj<T>(ellint_rj_data, type_name, "RJ: Random data");
}

template <typename T>
void t8(T, const char* type_name)
{
#include "ellint_rj_e4.ipp"

   do_test_ellint_rj<T>(ellint_rj_e4, type_name, "RJ: 4 Equal Values");
}

template <typename T>
void t9(T, const char* type_name)
{
#include "ellint_rj_e3.ipp"

   do_test_ellint_rj<T>(ellint_rj_e3, type_name, "RJ: 3 Equal Values");
}

template <typename T>
void t10(T, const char* type_name)
{
#include "ellint_rj_e2.ipp"

   do_test_ellint_rj<T>(ellint_rj_e2, type_name, "RJ: 2 Equal Values");
}

template <typename T>
void t11(T, const char* type_name)
{
#include "ellint_rj_zp.ipp"

   do_test_ellint_rj<T>(ellint_rj_zp, type_name, "RJ: Equal z and p");
}

#endif
#ifdef TEST3

template <typename T>
void t12(T, const char* type_name)
{
#include "ellint_rd_data.ipp"

   do_test_ellint_rd<T>(ellint_rd_data, type_name, "RD: Random data");
}

template <typename T>
void t13(T, const char* type_name)
{
#include "ellint_rd_xyy.ipp"

   do_test_ellint_rd<T>(ellint_rd_xyy, type_name, "RD: y = z");
}

template <typename T>
void t14(T, const char* type_name)
{
#include "ellint_rd_xxz.ipp"

   do_test_ellint_rd<T>(ellint_rd_xxz, type_name, "RD: x = y");
}

template <typename T>
void t15(T, const char* type_name)
{
#include "ellint_rd_0yy.ipp"

   do_test_ellint_rd<T>(ellint_rd_0yy, type_name, "RD: x = 0, y = z");
}

template <typename T>
void t16(T, const char* type_name)
{
#include "ellint_rd_xxx.ipp"

   do_test_ellint_rd<T>(ellint_rd_xxx, type_name, "RD: x = y = z");
}

template <typename T>
void t17(T, const char* type_name)
{
#include "ellint_rd_0xy.ipp"

   do_test_ellint_rd<T>(ellint_rd_0xy, type_name, "RD: x = 0");
}

#endif
#ifdef TEST4

template <typename T>
void t18(T, const char* type_name)
{
#include "ellint_rg.ipp"

   do_test_ellint_rg<T>(ellint_rg, type_name, "RG: Random Data");
}

template <typename T>
void t19(T, const char* type_name)
{
#include "ellint_rg_00x.ipp"

   do_test_ellint_rg<T>(ellint_rg_00x, type_name, "RG: two values 0");
}

template <typename T>
void t20(T, const char* type_name)
{
#include "ellint_rg_xxx.ipp"

   do_test_ellint_rg<T>(ellint_rg_xxx, type_name, "RG: All values the same or zero");
}

template <typename T>
void t21(T, const char* type_name)
{
#include "ellint_rg_xyy.ipp"

   do_test_ellint_rg<T>(ellint_rg_xyy, type_name, "RG: two values the same");
}

template <typename T>
void t22(T, const char* type_name)
{
#include "ellint_rg_xy0.ipp"

   do_test_ellint_rg<T>(ellint_rg_xy0, type_name, "RG: one value zero");
}

#endif

template <typename T>
void test_spots(T val, const char* type_name)
{
#ifndef TEST_UDT
   using namespace boost::math;
   using namespace std;
   // Spot values from Numerical Computation of Real or Complex 
   // Elliptic Integrals, B. C. Carlson: http://arxiv.org/abs/math.CA/9409227
   // RF:
   T tolerance = (std::max)(T(1e-13f), tools::epsilon<T>() * 5) * 100; // Note 5eps expressed as a persentage!!!
   T eps2 = 5 * tools::epsilon<T>();
   BOOST_CHECK_CLOSE(ellint_rf(T(1), T(2), T(0)), T(1.3110287771461), tolerance);
   BOOST_CHECK_CLOSE(ellint_rf(T(0.5), T(1), T(0)), T(1.8540746773014), tolerance);
   BOOST_CHECK_CLOSE(ellint_rf(T(2), T(3), T(4)), T(0.58408284167715), tolerance);
   // RC:
   BOOST_CHECK_CLOSE_FRACTION(ellint_rc(T(0), T(1)/4), boost::math::constants::pi<T>(), eps2);
   BOOST_CHECK_CLOSE_FRACTION(ellint_rc(T(9)/4, T(2)), boost::math::constants::ln_two<T>(), eps2);
   BOOST_CHECK_CLOSE_FRACTION(ellint_rc(T(1) / 4, T(-2)), boost::math::constants::ln_two<T>() / 3, eps2);
   // RJ:
   BOOST_CHECK_CLOSE(ellint_rj(T(0), T(1), T(2), T(3)), T(0.77688623778582), tolerance);
   BOOST_CHECK_CLOSE(ellint_rj(T(2), T(3), T(4), T(5)), T(0.14297579667157), tolerance);
   BOOST_CHECK_CLOSE(ellint_rj(T(2), T(3), T(4), T(-0.5)), T(0.24723819703052), tolerance);
   BOOST_CHECK_CLOSE(ellint_rj(T(2), T(3), T(4), T(-5)), T(-0.12711230042964), tolerance);
   // RD:
   BOOST_CHECK_CLOSE(ellint_rd(T(0), T(2), T(1)), T(1.7972103521034), tolerance);
   BOOST_CHECK_CLOSE(ellint_rd(T(2), T(3), T(4)), T(0.16510527294261), tolerance);

   // Sanity/consistency checks from Numerical Computation of Real or Complex 
   // Elliptic Integrals, B. C. Carlson: http://arxiv.org/abs/math.CA/9409227
   boost::mt19937 ran;
   boost::uniform_real<float> ur(0, 1000);
   T eps40 = 40 * tools::epsilon<T>();

   for(unsigned i = 0; i < 1000; ++i)
   {
      T x = ur(ran);
      T y = ur(ran);
      T z = ur(ran);
      T lambda = ur(ran);
      T mu = x * y / lambda;
      // RF, eq 49:
      T s1 = ellint_rf(x+lambda, y+lambda, lambda) + 
         ellint_rf(x + mu, y + mu, mu);
      T s2 = ellint_rf(x, y, T(0));
      BOOST_CHECK_CLOSE_FRACTION(s1, s2, eps40);
      // RC is degenerate case of RF:
      s1 = ellint_rc(x, y);
      s2 = ellint_rf(x, y, y);
      BOOST_CHECK_CLOSE_FRACTION(s1, s2, eps40);
      // RC, eq 50 (Note have to assume y = x):
      T mu2 = x * x / lambda;
      s1 = ellint_rc(lambda, x+lambda) 
         + ellint_rc(mu2, x + mu2);
      s2 = ellint_rc(T(0), x);
      BOOST_CHECK_CLOSE_FRACTION(s1, s2, eps40);
      /*
      T p = ????; // no closed form for a, b and p???
      s1 = ellint_rj(x+lambda, y+lambda, lambda, p+lambda)
         + ellint_rj(x+mu, y+mu, mu, p+mu);
      s2 = ellint_rj(x, y, T(0), p)
         - 3 * ellint_rc(a, b);
      */
      // RD, eq 53:
      s1 = ellint_rd(lambda, x+lambda, y+lambda)
         + ellint_rd(mu, x+mu, y+mu);
      s2 = ellint_rd(T(0), x, y)
         - 3 / (y * sqrt(x+y+lambda+mu));
      BOOST_CHECK_CLOSE_FRACTION(s1, s2, eps40);
      // RD is degenerate case of RJ:
      s1 = ellint_rd(x, y, z);
      s2 = ellint_rj(x, y, z, z);
      BOOST_CHECK_CLOSE_FRACTION(s1, s2, eps40);
   }
#endif
   //
   // Now random spot values:
   //
#ifdef TEST1
   t1(val, type_name);
   t2(val, type_name);
   t3(val, type_name);
   t4(val, type_name);
   t5(val, type_name);
#endif
#ifdef TEST2
   t6(val, type_name);
   t7(val, type_name);
   t8(val, type_name);
   t9(val, type_name);
   t10(val, type_name);
   t11(val, type_name);
#endif
#ifdef TEST3
   t12(val, type_name);
   t13(val, type_name);
   t14(val, type_name);
   t15(val, type_name);
   t16(val, type_name);
   t17(val, type_name);
#endif
#ifdef TEST4
   t18(val, type_name);
   t19(val, type_name);
   t20(val, type_name);
   t21(val, type_name);
   t22(val, type_name);
#endif
}
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	#include <boost/math/concepts/real_concept.hpp>
	8	#define BOOST_TEST_MAIN
	9	#include <boost/test/unit_test.hpp>
92f5a8d4	10	#include <boost/test/tools/floating_point_comparison.hpp>
7c673cae FG	11	#include <boost/math/special_functions/math_fwd.hpp>
	12	#include <boost/math/constants/constants.hpp>
	13	#include <boost/array.hpp>
	14	#include <boost/random.hpp>
	15	#include "functor.hpp"
	16
	17	#include "handle_test_result.hpp"
	18	#include "table_type.hpp"
	19
	20	#ifndef SC_
	21	#define SC_(x) static_cast<typename table_type<T>::type>(BOOST_JOIN(x, L))
	22	#endif
	23
	24	template <class Real, typename T>
	25	void do_test_ellint_rf(T& data, const char* type_name, const char* test)
	26	{
	27	#if !(defined(ERROR_REPORTING_MODE) && !defined(ELLINT_RF_FUNCTION_TO_TEST))
	28	typedef Real value_type;
	29
	30	std::cout << "Testing: " << test << std::endl;
	31
	32	#ifdef ELLINT_RF_FUNCTION_TO_TEST
	33	value_type(*fp)(value_type, value_type, value_type) = ELLINT_RF_FUNCTION_TO_TEST;
	34	#elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
	35	value_type (*fp)(value_type, value_type, value_type) = boost::math::ellint_rf<value_type, value_type, value_type>;
	36	#else
	37	value_type (*fp)(value_type, value_type, value_type) = boost::math::ellint_rf;
	38	#endif
	39	boost::math::tools::test_result<value_type> result;
	40
	41	result = boost::math::tools::test_hetero<Real>(
	42	data,
	43	bind_func<Real>(fp, 0, 1, 2),
	44	extract_result<Real>(3));
	45	handle_test_result(result, data[result.worst()], result.worst(),
	46	type_name, "ellint_rf", test);
	47
	48	std::cout << std::endl;
	49	#endif
	50	}
	51
	52	template <class Real, typename T>
	53	void do_test_ellint_rc(T& data, const char* type_name, const char* test)
	54	{
	55	#if !(defined(ERROR_REPORTING_MODE) && !defined(ELLINT_RC_FUNCTION_TO_TEST))
	56	typedef Real value_type;
	57
	58	std::cout << "Testing: " << test << std::endl;
	59
	60	#ifdef ELLINT_RC_FUNCTION_TO_TEST
	61	value_type(*fp)(value_type, value_type) = ELLINT_RC_FUNCTION_TO_TEST;
	62	#elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
	63	value_type (*fp)(value_type, value_type) = boost::math::ellint_rc<value_type, value_type>;
	64	#else
	65	value_type (*fp)(value_type, value_type) = boost::math::ellint_rc;
	66	#endif
	67	boost::math::tools::test_result<value_type> result;
	68
	69	result = boost::math::tools::test_hetero<Real>(
	70	data,
	71	bind_func<Real>(fp, 0, 1),
	72	extract_result<Real>(2));
	73	handle_test_result(result, data[result.worst()], result.worst(),
	74	type_name, "ellint_rc", test);
75
76	std::cout << std::endl;
77	#endif
78	}
79
80	template <class Real, typename T>
81	void do_test_ellint_rj(T& data, const char* type_name, const char* test)
82	{
83	#if !(defined(ERROR_REPORTING_MODE) && !defined(ELLINT_RJ_FUNCTION_TO_TEST))
84	typedef Real value_type;
85
86	std::cout << "Testing: " << test << std::endl;
87
88	#ifdef ELLINT_RJ_FUNCTION_TO_TEST
89	value_type(*fp)(value_type, value_type, value_type, value_type) = ELLINT_RJ_FUNCTION_TO_TEST;
90	#elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
91	value_type (*fp)(value_type, value_type, value_type, value_type) = boost::math::ellint_rj<value_type, value_type, value_type, value_type>;
92	#else
93	value_type (*fp)(value_type, value_type, value_type, value_type) = boost::math::ellint_rj;
94	#endif
95	boost::math::tools::test_result<value_type> result;
96
97	result = boost::math::tools::test_hetero<Real>(
98	data,
99	bind_func<Real>(fp, 0, 1, 2, 3),
100	extract_result<Real>(4));
101	handle_test_result(result, data[result.worst()], result.worst(),
102	type_name, "ellint_rj", test);
103
104	std::cout << std::endl;
105	#endif
106	}
107
108	template <class Real, typename T>
109	void do_test_ellint_rd(T& data, const char* type_name, const char* test)
110	{
111	#if !(defined(ERROR_REPORTING_MODE) && !defined(ELLINT_RD_FUNCTION_TO_TEST))
112	typedef Real value_type;
113
114	std::cout << "Testing: " << test << std::endl;
115
116	#ifdef ELLINT_RD_FUNCTION_TO_TEST
117	value_type(*fp)(value_type, value_type, value_type) = ELLINT_RD_FUNCTION_TO_TEST;
118	#elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
119	value_type (*fp)(value_type, value_type, value_type) = boost::math::ellint_rd<value_type, value_type, value_type>;
120	#else
121	value_type (*fp)(value_type, value_type, value_type) = boost::math::ellint_rd;
122	#endif
123	boost::math::tools::test_result<value_type> result;
124
125	result = boost::math::tools::test_hetero<Real>(
126	data,
127	bind_func<Real>(fp, 0, 1, 2),
128	extract_result<Real>(3));
129	handle_test_result(result, data[result.worst()], result.worst(),
130	type_name, "ellint_rd", test);
131
132	std::cout << std::endl;
133	#endif
134	}
135
136	template <class Real, typename T>
137	void do_test_ellint_rg(T& data, const char* type_name, const char* test)
138	{
139	#if !(defined(ERROR_REPORTING_MODE) && !defined(ELLINT_RD_FUNCTION_TO_TEST))
140	typedef Real value_type;
141
142	std::cout << "Testing: " << test << std::endl;
143
144	#ifdef ELLINT_RG_FUNCTION_TO_TEST
145	value_type(*fp)(value_type, value_type, value_type) = ELLINT_RG_FUNCTION_TO_TEST;
146	#elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
147	value_type(*fp)(value_type, value_type, value_type) = boost::math::ellint_rg<value_type, value_type, value_type>;
148	#else
149	value_type(*fp)(value_type, value_type, value_type) = boost::math::ellint_rg;
150	#endif
151	boost::math::tools::test_result<value_type> result;
152
153	result = boost::math::tools::test_hetero<Real>(
154	data,
155	bind_func<Real>(fp, 0, 1, 2),
156	extract_result<Real>(3));
157	handle_test_result(result, data[result.worst()], result.worst(),
158	type_name, "ellint_rg", test);
159
160	std::cout << std::endl;
161	#endif
162	}
163
164	#if !defined(TEST1) && !defined(TEST2) && !defined(TEST3) && !defined(TEST4)
165	#define TEST1
166	#define TEST2
167	#define TEST3
168	#define TEST4
169	#endif
170
171	#ifdef TEST1
172
173	template <typename T>
174	void t1(T, const char* type_name)
175	{
176	#include "ellint_rf_data.ipp"
177
178	do_test_ellint_rf<T>(ellint_rf_data, type_name, "RF: Random data");
179	}
180
181	template <typename T>
182	void t2(T, const char* type_name)
183	{
184	#include "ellint_rf_xxx.ipp"
185
186	do_test_ellint_rf<T>(ellint_rf_xxx, type_name, "RF: x = y = z");
187	}
188
189	template <typename T>
190	void t3(T, const char* type_name)
191	{
192	#include "ellint_rf_xyy.ipp"
193
194	do_test_ellint_rf<T>(ellint_rf_xyy, type_name, "RF: x = y or y = z or x = z");
195	}
196
197	template <typename T>
198	void t4(T, const char* type_name)
199	{
200	#include "ellint_rf_0yy.ipp"
201
202	do_test_ellint_rf<T>(ellint_rf_0yy, type_name, "RF: x = 0, y = z");
203	}
204
205	template <typename T>
206	void t5(T, const char* type_name)
207	{
208	#include "ellint_rf_xy0.ipp"
209
210	do_test_ellint_rf<T>(ellint_rf_xy0, type_name, "RF: z = 0");
211	}
212
213	#endif
214	#ifdef TEST2
215
216	template <typename T>
217	void t6(T, const char* type_name)
218	{
219	#include "ellint_rc_data.ipp"
220
221	do_test_ellint_rc<T>(ellint_rc_data, type_name, "RC: Random data");
222	}
223
224	template <typename T>
225	void t7(T, const char* type_name)
226	{
227	#include "ellint_rj_data.ipp"
228
229	do_test_ellint_rj<T>(ellint_rj_data, type_name, "RJ: Random data");
230	}
231
232	template <typename T>
233	void t8(T, const char* type_name)
234	{
235	#include "ellint_rj_e4.ipp"
236
237	do_test_ellint_rj<T>(ellint_rj_e4, type_name, "RJ: 4 Equal Values");
238	}
239
240	template <typename T>
241	void t9(T, const char* type_name)
242	{
243	#include "ellint_rj_e3.ipp"
244
245	do_test_ellint_rj<T>(ellint_rj_e3, type_name, "RJ: 3 Equal Values");
246	}
247
248	template <typename T>
249	void t10(T, const char* type_name)
250	{
251	#include "ellint_rj_e2.ipp"
252
253	do_test_ellint_rj<T>(ellint_rj_e2, type_name, "RJ: 2 Equal Values");
254	}
255
256	template <typename T>
257	void t11(T, const char* type_name)
258	{
259	#include "ellint_rj_zp.ipp"
260
261	do_test_ellint_rj<T>(ellint_rj_zp, type_name, "RJ: Equal z and p");
262	}
263
264	#endif
265	#ifdef TEST3
266
267	template <typename T>
268	void t12(T, const char* type_name)
269	{
270	#include "ellint_rd_data.ipp"
271
272	do_test_ellint_rd<T>(ellint_rd_data, type_name, "RD: Random data");
273	}
274
275	template <typename T>
276	void t13(T, const char* type_name)
277	{
278	#include "ellint_rd_xyy.ipp"
279
280	do_test_ellint_rd<T>(ellint_rd_xyy, type_name, "RD: y = z");
281	}
282
283	template <typename T>
284	void t14(T, const char* type_name)
285	{
286	#include "ellint_rd_xxz.ipp"
287
288	do_test_ellint_rd<T>(ellint_rd_xxz, type_name, "RD: x = y");
289	}
290
291	template <typename T>
292	void t15(T, const char* type_name)
293	{
294	#include "ellint_rd_0yy.ipp"
295
296	do_test_ellint_rd<T>(ellint_rd_0yy, type_name, "RD: x = 0, y = z");
297	}
298
299	template <typename T>
300	void t16(T, const char* type_name)
301	{
302	#include "ellint_rd_xxx.ipp"
303
304	do_test_ellint_rd<T>(ellint_rd_xxx, type_name, "RD: x = y = z");
305	}
306
307	template <typename T>
308	void t17(T, const char* type_name)
309	{
310	#include "ellint_rd_0xy.ipp"
311
312	do_test_ellint_rd<T>(ellint_rd_0xy, type_name, "RD: x = 0");
313	}
314
315	#endif
316	#ifdef TEST4
317
318	template <typename T>
319	void t18(T, const char* type_name)
320	{
321	#include "ellint_rg.ipp"
322
323	do_test_ellint_rg<T>(ellint_rg, type_name, "RG: Random Data");
324	}
325
326	template <typename T>
327	void t19(T, const char* type_name)
328	{
329	#include "ellint_rg_00x.ipp"
330
331	do_test_ellint_rg<T>(ellint_rg_00x, type_name, "RG: two values 0");
332	}
333
334	template <typename T>
335	void t20(T, const char* type_name)
336	{
337	#include "ellint_rg_xxx.ipp"
338
339	do_test_ellint_rg<T>(ellint_rg_xxx, type_name, "RG: All values the same or zero");
340	}
341
342	template <typename T>
343	void t21(T, const char* type_name)
344	{
345	#include "ellint_rg_xyy.ipp"
346
347	do_test_ellint_rg<T>(ellint_rg_xyy, type_name, "RG: two values the same");
348	}
349
350	template <typename T>
351	void t22(T, const char* type_name)
352	{
353	#include "ellint_rg_xy0.ipp"
354
355	do_test_ellint_rg<T>(ellint_rg_xy0, type_name, "RG: one value zero");
356	}
357
358	#endif
359
360	template <typename T>
361	void test_spots(T val, const char* type_name)
362	{
363	#ifndef TEST_UDT
364	using namespace boost::math;
365	using namespace std;
366	// Spot values from Numerical Computation of Real or Complex
367	// Elliptic Integrals, B. C. Carlson: http://arxiv.org/abs/math.CA/9409227
368	// RF:
369	T tolerance = (std::max)(T(1e-13f), tools::epsilon<T>() * 5) * 100; // Note 5eps expressed as a persentage!!!
370	T eps2 = 5 * tools::epsilon<T>();
371	BOOST_CHECK_CLOSE(ellint_rf(T(1), T(2), T(0)), T(1.3110287771461), tolerance);
372	BOOST_CHECK_CLOSE(ellint_rf(T(0.5), T(1), T(0)), T(1.8540746773014), tolerance);
373	BOOST_CHECK_CLOSE(ellint_rf(T(2), T(3), T(4)), T(0.58408284167715), tolerance);
374	// RC:
375	BOOST_CHECK_CLOSE_FRACTION(ellint_rc(T(0), T(1)/4), boost::math::constants::pi<T>(), eps2);
376	BOOST_CHECK_CLOSE_FRACTION(ellint_rc(T(9)/4, T(2)), boost::math::constants::ln_two<T>(), eps2);
377	BOOST_CHECK_CLOSE_FRACTION(ellint_rc(T(1) / 4, T(-2)), boost::math::constants::ln_two<T>() / 3, eps2);
378	// RJ:
379	BOOST_CHECK_CLOSE(ellint_rj(T(0), T(1), T(2), T(3)), T(0.77688623778582), tolerance);
380	BOOST_CHECK_CLOSE(ellint_rj(T(2), T(3), T(4), T(5)), T(0.14297579667157), tolerance);
381	BOOST_CHECK_CLOSE(ellint_rj(T(2), T(3), T(4), T(-0.5)), T(0.24723819703052), tolerance);
382	BOOST_CHECK_CLOSE(ellint_rj(T(2), T(3), T(4), T(-5)), T(-0.12711230042964), tolerance);
383	// RD:
384	BOOST_CHECK_CLOSE(ellint_rd(T(0), T(2), T(1)), T(1.7972103521034), tolerance);
385	BOOST_CHECK_CLOSE(ellint_rd(T(2), T(3), T(4)), T(0.16510527294261), tolerance);
386
387	// Sanity/consistency checks from Numerical Computation of Real or Complex
388	// Elliptic Integrals, B. C. Carlson: http://arxiv.org/abs/math.CA/9409227
389	boost::mt19937 ran;
390	boost::uniform_real<float> ur(0, 1000);
391	T eps40 = 40 * tools::epsilon<T>();
392
393	for(unsigned i = 0; i < 1000; ++i)
394	{
395	T x = ur(ran);
396	T y = ur(ran);
397	T z = ur(ran);
398	T lambda = ur(ran);
399	T mu = x * y / lambda;
400	// RF, eq 49:
401	T s1 = ellint_rf(x+lambda, y+lambda, lambda) +
402	ellint_rf(x + mu, y + mu, mu);
403	T s2 = ellint_rf(x, y, T(0));
404	BOOST_CHECK_CLOSE_FRACTION(s1, s2, eps40);
405	// RC is degenerate case of RF:
406	s1 = ellint_rc(x, y);
407	s2 = ellint_rf(x, y, y);
408	BOOST_CHECK_CLOSE_FRACTION(s1, s2, eps40);
409	// RC, eq 50 (Note have to assume y = x):
410	T mu2 = x * x / lambda;
411	s1 = ellint_rc(lambda, x+lambda)
412	+ ellint_rc(mu2, x + mu2);
413	s2 = ellint_rc(T(0), x);
414	BOOST_CHECK_CLOSE_FRACTION(s1, s2, eps40);
415	/*
416	T p = ????; // no closed form for a, b and p???
417	s1 = ellint_rj(x+lambda, y+lambda, lambda, p+lambda)
418	+ ellint_rj(x+mu, y+mu, mu, p+mu);
419	s2 = ellint_rj(x, y, T(0), p)
420	- 3 * ellint_rc(a, b);
421	*/
422	// RD, eq 53:
423	s1 = ellint_rd(lambda, x+lambda, y+lambda)
424	+ ellint_rd(mu, x+mu, y+mu);
425	s2 = ellint_rd(T(0), x, y)
426	- 3 / (y * sqrt(x+y+lambda+mu));
427	BOOST_CHECK_CLOSE_FRACTION(s1, s2, eps40);
428	// RD is degenerate case of RJ:
429	s1 = ellint_rd(x, y, z);
430	s2 = ellint_rj(x, y, z, z);
431	BOOST_CHECK_CLOSE_FRACTION(s1, s2, eps40);
432	}
433	#endif
434	//
435	// Now random spot values:
436	//
437	#ifdef TEST1
438	t1(val, type_name);
439	t2(val, type_name);
440	t3(val, type_name);
441	t4(val, type_name);
442	t5(val, type_name);
443	#endif
444	#ifdef TEST2
445	t6(val, type_name);
446	t7(val, type_name);
447	t8(val, type_name);
448	t9(val, type_name);
449	t10(val, type_name);
450	t11(val, type_name);
451	#endif
452	#ifdef TEST3
453	t12(val, type_name);
454	t13(val, type_name);
455	t14(val, type_name);
456	t15(val, type_name);
457	t16(val, type_name);
458	t17(val, type_name);
459	#endif
460	#ifdef TEST4
461	t18(val, type_name);
462	t19(val, type_name);
463	t20(val, type_name);
464	t21(val, type_name);
465	t22(val, type_name);
466	#endif
467	}
468