[ceph.git] / ceph / src / boost / libs / math / include / boost / math / special_functions / jacobi_elliptic.hpp

// Copyright John Maddock 2012.
// 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)

#ifndef BOOST_MATH_JACOBI_ELLIPTIC_HPP
#define BOOST_MATH_JACOBI_ELLIPTIC_HPP

#include <boost/math/tools/precision.hpp>
#include <boost/math/tools/promotion.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/special_functions/math_fwd.hpp>

namespace boost{ namespace math{

namespace detail{

template <class T, class Policy>
T jacobi_recurse(const T& x, const T& k, T anm1, T bnm1, unsigned N, T* pTn, const Policy& pol)
{
   BOOST_MATH_STD_USING
   ++N;
   T Tn;
   T cn = (anm1 - bnm1) / 2;
   T an = (anm1 + bnm1) / 2;
   if(cn < policies::get_epsilon<T, Policy>())
   {
      Tn = ldexp(T(1), (int)N) * x * an;
   }
   else
      Tn = jacobi_recurse<T>(x, k, an, sqrt(anm1 * bnm1), N, 0, pol);
   if(pTn)
      *pTn = Tn;
   return (Tn + asin((cn / an) * sin(Tn))) / 2;
}

template <class T, class Policy>
T jacobi_imp(const T& x, const T& k, T* cn, T* dn, const Policy& pol, const char* function)
{
   BOOST_MATH_STD_USING
   if(k < 0)
   {
      *cn = policies::raise_domain_error<T>(function, "Modulus k must be positive but got %1%.", k, pol);
      *dn = *cn;
      return *cn;
   }
   if(k > 1)
   {
      T xp = x * k;
      T kp = 1 / k;
      T snp, cnp, dnp;
      snp = jacobi_imp(xp, kp, &cnp, &dnp, pol, function);
      *cn = dnp;
      *dn = cnp;
      return snp * kp;
   }
   //
   // Special cases first:
   //
   if(x == 0)
   {
      *cn = *dn = 1;
      return 0;
   }
   if(k == 0)
   {
      *cn = cos(x);
      *dn = 1;
      return sin(x);
   }
   if(k == 1)
   {
      *cn = *dn = 1 / cosh(x);
      return tanh(x);
   }
   //
   // Asymptotic forms from A&S 16.13:
   //
   if(k < tools::forth_root_epsilon<T>())
   {
      T su = sin(x);
      T cu = cos(x);
      T m = k * k;
      *dn = 1 - m * su * su / 2;
      *cn = cu + m * (x - su * cu) * su / 4;
      return su - m * (x - su * cu) * cu / 4;
   }
   /*  Can't get this to work to adequate precision - disabled for now...
   //
   // Asymptotic forms from A&S 16.15:
   //
   if(k > 1 - tools::root_epsilon<T>())
   {
      T tu = tanh(x);
      T su = sinh(x);
      T cu = cosh(x);
      T sec = 1 / cu;
      T kp = 1 - k;
      T m1 = 2 * kp - kp * kp;
      *dn = sec + m1 * (su * cu + x) * tu * sec / 4;
      *cn = sec - m1 * (su * cu - x) * tu * sec / 4;
      T sn = tu;
      T sn2 = m1 * (x * sec * sec - tu) / 4;
      T sn3 = (72 * x * cu + 4 * (8 * x * x - 5) * su - 19 * sinh(3 * x) + sinh(5 * x)) * sec * sec * sec * m1 * m1 / 512;
      return sn + sn2 - sn3;
   }*/
   T T1;
   T kc = 1 - k;
   T k_prime = k < 0.5 ? T(sqrt(1 - k * k)) : T(sqrt(2 * kc - kc * kc));
   T T0 = jacobi_recurse(x, k, T(1), k_prime, 0, &T1, pol);
   *cn = cos(T0);
   *dn = cos(T0) / cos(T1 - T0);
   return sin(T0);
}

} // namespace detail

template <class T, class U, class V, class Policy>
inline typename tools::promote_args<T, U, V>::type jacobi_elliptic(T k, U theta, V* pcn, V* pdn, const Policy&)
{
   BOOST_FPU_EXCEPTION_GUARD
   typedef typename tools::promote_args<T>::type result_type;
   typedef typename policies::evaluation<result_type, Policy>::type value_type;
   typedef typename policies::normalise<
      Policy, 
      policies::promote_float<false>, 
      policies::promote_double<false>, 
      policies::discrete_quantile<>,
      policies::assert_undefined<> >::type forwarding_policy;

   static const char* function = "boost::math::jacobi_elliptic<%1%>(%1%)";

   value_type sn, cn, dn;
   sn = detail::jacobi_imp<value_type>(static_cast<value_type>(theta), static_cast<value_type>(k), &cn, &dn, forwarding_policy(), function);
   if(pcn)
      *pcn = policies::checked_narrowing_cast<result_type, Policy>(cn, function);
   if(pdn)
      *pdn = policies::checked_narrowing_cast<result_type, Policy>(dn, function);
   return policies::checked_narrowing_cast<result_type, Policy>(sn, function);;
}

template <class T, class U, class V>
inline typename tools::promote_args<T, U, V>::type jacobi_elliptic(T k, U theta, V* pcn, V* pdn)
{
   return jacobi_elliptic(k, theta, pcn, pdn, policies::policy<>());
}

template <class U, class T, class Policy>
inline typename tools::promote_args<T, U>::type jacobi_sn(U k, T theta, const Policy& pol)
{
   typedef typename tools::promote_args<T, U>::type result_type;
   return jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), static_cast<result_type*>(0), static_cast<result_type*>(0), pol);
}

template <class U, class T>
inline typename tools::promote_args<T, U>::type jacobi_sn(U k, T theta)
{
   return jacobi_sn(k, theta, policies::policy<>());
}

template <class T, class U, class Policy>
inline typename tools::promote_args<T, U>::type jacobi_cn(T k, U theta, const Policy& pol)
{
   typedef typename tools::promote_args<T, U>::type result_type;
   result_type cn;
   jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), &cn, static_cast<result_type*>(0), pol);
   return cn;
}

template <class T, class U>
inline typename tools::promote_args<T, U>::type jacobi_cn(T k, U theta)
{
   return jacobi_cn(k, theta, policies::policy<>());
}

template <class T, class U, class Policy>
inline typename tools::promote_args<T, U>::type jacobi_dn(T k, U theta, const Policy& pol)
{
   typedef typename tools::promote_args<T, U>::type result_type;
   result_type dn;
   jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), static_cast<result_type*>(0), &dn, pol);
   return dn;
}

template <class T, class U>
inline typename tools::promote_args<T, U>::type jacobi_dn(T k, U theta)
{
   return jacobi_dn(k, theta, policies::policy<>());
}

template <class T, class U, class Policy>
inline typename tools::promote_args<T, U>::type jacobi_cd(T k, U theta, const Policy& pol)
{
   typedef typename tools::promote_args<T, U>::type result_type;
   result_type cn, dn;
   jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), &cn, &dn, pol);
   return cn / dn;
}

template <class T, class U>
inline typename tools::promote_args<T, U>::type jacobi_cd(T k, U theta)
{
   return jacobi_cd(k, theta, policies::policy<>());
}

template <class T, class U, class Policy>
inline typename tools::promote_args<T, U>::type jacobi_dc(T k, U theta, const Policy& pol)
{
   typedef typename tools::promote_args<T, U>::type result_type;
   result_type cn, dn;
   jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), &cn, &dn, pol);
   return dn / cn;
}

template <class T, class U>
inline typename tools::promote_args<T, U>::type jacobi_dc(T k, U theta)
{
   return jacobi_dc(k, theta, policies::policy<>());
}

template <class T, class U, class Policy>
inline typename tools::promote_args<T, U>::type jacobi_ns(T k, U theta, const Policy& pol)
{
   typedef typename tools::promote_args<T, U>::type result_type;
   return 1 / jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), static_cast<result_type*>(0), static_cast<result_type*>(0), pol);
}

template <class T, class U>
inline typename tools::promote_args<T, U>::type jacobi_ns(T k, U theta)
{
   return jacobi_ns(k, theta, policies::policy<>());
}

template <class T, class U, class Policy>
inline typename tools::promote_args<T, U>::type jacobi_sd(T k, U theta, const Policy& pol)
{
   typedef typename tools::promote_args<T, U>::type result_type;
   result_type sn, dn;
   sn = jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), static_cast<result_type*>(0), &dn, pol);
   return sn / dn;
}

template <class T, class U>
inline typename tools::promote_args<T, U>::type jacobi_sd(T k, U theta)
{
   return jacobi_sd(k, theta, policies::policy<>());
}

template <class T, class U, class Policy>
inline typename tools::promote_args<T, U>::type jacobi_ds(T k, U theta, const Policy& pol)
{
   typedef typename tools::promote_args<T, U>::type result_type;
   result_type sn, dn;
   sn = jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), static_cast<result_type*>(0), &dn, pol);
   return dn / sn;
}

template <class T, class U>
inline typename tools::promote_args<T, U>::type jacobi_ds(T k, U theta)
{
   return jacobi_ds(k, theta, policies::policy<>());
}

template <class T, class U, class Policy>
inline typename tools::promote_args<T, U>::type jacobi_nc(T k, U theta, const Policy& pol)
{
   return 1 / jacobi_cn(k, theta, pol);
}

template <class T, class U>
inline typename tools::promote_args<T, U>::type jacobi_nc(T k, U theta)
{
   return jacobi_nc(k, theta, policies::policy<>());
}

template <class T, class U, class Policy>
inline typename tools::promote_args<T, U>::type jacobi_nd(T k, U theta, const Policy& pol)
{
   return 1 / jacobi_dn(k, theta, pol);
}

template <class T, class U>
inline typename tools::promote_args<T, U>::type jacobi_nd(T k, U theta)
{
   return jacobi_nd(k, theta, policies::policy<>());
}

template <class T, class U, class Policy>
inline typename tools::promote_args<T, U>::type jacobi_sc(T k, U theta, const Policy& pol)
{
   typedef typename tools::promote_args<T, U>::type result_type;
   result_type sn, cn;
   sn = jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), &cn, static_cast<result_type*>(0), pol);
   return sn / cn;
}

template <class T, class U>
inline typename tools::promote_args<T, U>::type jacobi_sc(T k, U theta)
{
   return jacobi_sc(k, theta, policies::policy<>());
}

template <class T, class U, class Policy>
inline typename tools::promote_args<T, U>::type jacobi_cs(T k, U theta, const Policy& pol)
{
   typedef typename tools::promote_args<T, U>::type result_type;
   result_type sn, cn;
   sn = jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), &cn, static_cast<result_type*>(0), pol);
   return cn / sn;
}

template <class T, class U>
inline typename tools::promote_args<T, U>::type jacobi_cs(T k, U theta)
{
   return jacobi_cs(k, theta, policies::policy<>());
}

}} // namespaces

#endif // BOOST_MATH_JACOBI_ELLIPTIC_HPP
Commit	Line	Data
7c673cae FG	1	// Copyright John Maddock 2012.
	2	// Use, modification and distribution are subject to the
	3	// Boost Software License, Version 1.0.
	4	// (See accompanying file LICENSE_1_0.txt
	5	// or copy at http://www.boost.org/LICENSE_1_0.txt)
	6
	7	#ifndef BOOST_MATH_JACOBI_ELLIPTIC_HPP
	8	#define BOOST_MATH_JACOBI_ELLIPTIC_HPP
	9
	10	#include <boost/math/tools/precision.hpp>
	11	#include <boost/math/tools/promotion.hpp>
	12	#include <boost/math/policies/error_handling.hpp>
	13	#include <boost/math/special_functions/math_fwd.hpp>
	14
	15	namespace boost{ namespace math{
	16
	17	namespace detail{
	18
	19	template <class T, class Policy>
	20	T jacobi_recurse(const T& x, const T& k, T anm1, T bnm1, unsigned N, T* pTn, const Policy& pol)
	21	{
	22	BOOST_MATH_STD_USING
	23	++N;
	24	T Tn;
	25	T cn = (anm1 - bnm1) / 2;
	26	T an = (anm1 + bnm1) / 2;
	27	if(cn < policies::get_epsilon<T, Policy>())
	28	{
	29	Tn = ldexp(T(1), (int)N) * x * an;
	30	}
	31	else
	32	Tn = jacobi_recurse<T>(x, k, an, sqrt(anm1 * bnm1), N, 0, pol);
	33	if(pTn)
	34	*pTn = Tn;
	35	return (Tn + asin((cn / an) * sin(Tn))) / 2;
	36	}
	37
	38	template <class T, class Policy>
	39	T jacobi_imp(const T& x, const T& k, T* cn, T* dn, const Policy& pol, const char* function)
	40	{
	41	BOOST_MATH_STD_USING
	42	if(k < 0)
	43	{
	44	*cn = policies::raise_domain_error<T>(function, "Modulus k must be positive but got %1%.", k, pol);
	45	dn = cn;
	46	return *cn;
	47	}
	48	if(k > 1)
	49	{
	50	T xp = x * k;
	51	T kp = 1 / k;
	52	T snp, cnp, dnp;
	53	snp = jacobi_imp(xp, kp, &cnp, &dnp, pol, function);
	54	*cn = dnp;
	55	*dn = cnp;
	56	return snp * kp;
	57	}
	58	//
	59	// Special cases first:
	60	//
	61	if(x == 0)
	62	{
	63	cn = dn = 1;
	64	return 0;
65	}
66	if(k == 0)
67	{
68	*cn = cos(x);
69	*dn = 1;
70	return sin(x);
71	}
72	if(k == 1)
73	{
74	cn = dn = 1 / cosh(x);
75	return tanh(x);
76	}
77	//
78	// Asymptotic forms from A&S 16.13:
79	//
80	if(k < tools::forth_root_epsilon<T>())
81	{
82	T su = sin(x);
83	T cu = cos(x);
84	T m = k * k;
85	dn = 1 - m su * su / 2;
86	cn = cu + m (x - su * cu) * su / 4;
87	return su - m * (x - su * cu) * cu / 4;
88	}
89	/* Can't get this to work to adequate precision - disabled for now...
90	//
91	// Asymptotic forms from A&S 16.15:
92	//
93	if(k > 1 - tools::root_epsilon<T>())
94	{
95	T tu = tanh(x);
96	T su = sinh(x);
97	T cu = cosh(x);
98	T sec = 1 / cu;
99	T kp = 1 - k;
100	T m1 = 2 * kp - kp * kp;
101	dn = sec + m1 (su * cu + x) * tu * sec / 4;
102	cn = sec - m1 (su * cu - x) * tu * sec / 4;
103	T sn = tu;
104	T sn2 = m1 * (x * sec * sec - tu) / 4;
105	T sn3 = (72 * x * cu + 4 * (8 * x * x - 5) * su - 19 * sinh(3 * x) + sinh(5 * x)) * sec * sec * sec * m1 * m1 / 512;
106	return sn + sn2 - sn3;
107	}*/
108	T T1;
109	T kc = 1 - k;
110	T k_prime = k < 0.5 ? T(sqrt(1 - k * k)) : T(sqrt(2 * kc - kc * kc));
111	T T0 = jacobi_recurse(x, k, T(1), k_prime, 0, &T1, pol);
112	*cn = cos(T0);
113	*dn = cos(T0) / cos(T1 - T0);
114	return sin(T0);
115	}
116
117	} // namespace detail
118
119	template <class T, class U, class V, class Policy>
120	inline typename tools::promote_args<T, U, V>::type jacobi_elliptic(T k, U theta, V* pcn, V* pdn, const Policy&)
121	{
122	BOOST_FPU_EXCEPTION_GUARD
123	typedef typename tools::promote_args<T>::type result_type;
124	typedef typename policies::evaluation<result_type, Policy>::type value_type;
125	typedef typename policies::normalise<
126	Policy,
127	policies::promote_float<false>,
128	policies::promote_double<false>,
129	policies::discrete_quantile<>,
130	policies::assert_undefined<> >::type forwarding_policy;
131
132	static const char* function = "boost::math::jacobi_elliptic<%1%>(%1%)";
133
134	value_type sn, cn, dn;
135	sn = detail::jacobi_imp<value_type>(static_cast<value_type>(theta), static_cast<value_type>(k), &cn, &dn, forwarding_policy(), function);
136	if(pcn)
137	*pcn = policies::checked_narrowing_cast<result_type, Policy>(cn, function);
138	if(pdn)
139	*pdn = policies::checked_narrowing_cast<result_type, Policy>(dn, function);
140	return policies::checked_narrowing_cast<result_type, Policy>(sn, function);;
141	}
142
143	template <class T, class U, class V>
144	inline typename tools::promote_args<T, U, V>::type jacobi_elliptic(T k, U theta, V* pcn, V* pdn)
145	{
146	return jacobi_elliptic(k, theta, pcn, pdn, policies::policy<>());
147	}
148
149	template <class U, class T, class Policy>
150	inline typename tools::promote_args<T, U>::type jacobi_sn(U k, T theta, const Policy& pol)
151	{
152	typedef typename tools::promote_args<T, U>::type result_type;
153	return jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), static_cast<result_type>(0), static_cast<result_type>(0), pol);
154	}
155
156	template <class U, class T>
157	inline typename tools::promote_args<T, U>::type jacobi_sn(U k, T theta)
158	{
159	return jacobi_sn(k, theta, policies::policy<>());
160	}
161
162	template <class T, class U, class Policy>
163	inline typename tools::promote_args<T, U>::type jacobi_cn(T k, U theta, const Policy& pol)
164	{
165	typedef typename tools::promote_args<T, U>::type result_type;
166	result_type cn;
167	jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), &cn, static_cast<result_type*>(0), pol);
168	return cn;
169	}
170
171	template <class T, class U>
172	inline typename tools::promote_args<T, U>::type jacobi_cn(T k, U theta)
173	{
174	return jacobi_cn(k, theta, policies::policy<>());
175	}
176
177	template <class T, class U, class Policy>
178	inline typename tools::promote_args<T, U>::type jacobi_dn(T k, U theta, const Policy& pol)
179	{
180	typedef typename tools::promote_args<T, U>::type result_type;
181	result_type dn;
182	jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), static_cast<result_type*>(0), &dn, pol);
183	return dn;
184	}
185
186	template <class T, class U>
187	inline typename tools::promote_args<T, U>::type jacobi_dn(T k, U theta)
188	{
189	return jacobi_dn(k, theta, policies::policy<>());
190	}
191
192	template <class T, class U, class Policy>
193	inline typename tools::promote_args<T, U>::type jacobi_cd(T k, U theta, const Policy& pol)
194	{
195	typedef typename tools::promote_args<T, U>::type result_type;
196	result_type cn, dn;
197	jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), &cn, &dn, pol);
198	return cn / dn;
199	}
200
201	template <class T, class U>
202	inline typename tools::promote_args<T, U>::type jacobi_cd(T k, U theta)
203	{
204	return jacobi_cd(k, theta, policies::policy<>());
205	}
206
207	template <class T, class U, class Policy>
208	inline typename tools::promote_args<T, U>::type jacobi_dc(T k, U theta, const Policy& pol)
209	{
210	typedef typename tools::promote_args<T, U>::type result_type;
211	result_type cn, dn;
212	jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), &cn, &dn, pol);
213	return dn / cn;
214	}
215
216	template <class T, class U>
217	inline typename tools::promote_args<T, U>::type jacobi_dc(T k, U theta)
218	{
219	return jacobi_dc(k, theta, policies::policy<>());
220	}
221
222	template <class T, class U, class Policy>
223	inline typename tools::promote_args<T, U>::type jacobi_ns(T k, U theta, const Policy& pol)
224	{
225	typedef typename tools::promote_args<T, U>::type result_type;
226	return 1 / jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), static_cast<result_type>(0), static_cast<result_type>(0), pol);
227	}
228
229	template <class T, class U>
230	inline typename tools::promote_args<T, U>::type jacobi_ns(T k, U theta)
231	{
232	return jacobi_ns(k, theta, policies::policy<>());
233	}
234
235	template <class T, class U, class Policy>
236	inline typename tools::promote_args<T, U>::type jacobi_sd(T k, U theta, const Policy& pol)
237	{
238	typedef typename tools::promote_args<T, U>::type result_type;
239	result_type sn, dn;
240	sn = jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), static_cast<result_type*>(0), &dn, pol);
241	return sn / dn;
242	}
243
244	template <class T, class U>
245	inline typename tools::promote_args<T, U>::type jacobi_sd(T k, U theta)
246	{
247	return jacobi_sd(k, theta, policies::policy<>());
248	}
249
250	template <class T, class U, class Policy>
251	inline typename tools::promote_args<T, U>::type jacobi_ds(T k, U theta, const Policy& pol)
252	{
253	typedef typename tools::promote_args<T, U>::type result_type;
254	result_type sn, dn;
255	sn = jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), static_cast<result_type*>(0), &dn, pol);
256	return dn / sn;
257	}
258
259	template <class T, class U>
260	inline typename tools::promote_args<T, U>::type jacobi_ds(T k, U theta)
261	{
262	return jacobi_ds(k, theta, policies::policy<>());
263	}
264
265	template <class T, class U, class Policy>
266	inline typename tools::promote_args<T, U>::type jacobi_nc(T k, U theta, const Policy& pol)
267	{
268	return 1 / jacobi_cn(k, theta, pol);
269	}
270
271	template <class T, class U>
272	inline typename tools::promote_args<T, U>::type jacobi_nc(T k, U theta)
273	{
274	return jacobi_nc(k, theta, policies::policy<>());
275	}
276
277	template <class T, class U, class Policy>
278	inline typename tools::promote_args<T, U>::type jacobi_nd(T k, U theta, const Policy& pol)
279	{
280	return 1 / jacobi_dn(k, theta, pol);
281	}
282
283	template <class T, class U>
284	inline typename tools::promote_args<T, U>::type jacobi_nd(T k, U theta)
285	{
286	return jacobi_nd(k, theta, policies::policy<>());
287	}
288
289	template <class T, class U, class Policy>
290	inline typename tools::promote_args<T, U>::type jacobi_sc(T k, U theta, const Policy& pol)
291	{
292	typedef typename tools::promote_args<T, U>::type result_type;
293	result_type sn, cn;
294	sn = jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), &cn, static_cast<result_type*>(0), pol);
295	return sn / cn;
296	}
297
298	template <class T, class U>
299	inline typename tools::promote_args<T, U>::type jacobi_sc(T k, U theta)
300	{
301	return jacobi_sc(k, theta, policies::policy<>());
302	}
303
304	template <class T, class U, class Policy>
305	inline typename tools::promote_args<T, U>::type jacobi_cs(T k, U theta, const Policy& pol)
306	{
307	typedef typename tools::promote_args<T, U>::type result_type;
308	result_type sn, cn;
309	sn = jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), &cn, static_cast<result_type*>(0), pol);
310	return cn / sn;
311	}
312
313	template <class T, class U>
314	inline typename tools::promote_args<T, U>::type jacobi_cs(T k, U theta)
315	{
316	return jacobi_cs(k, theta, policies::policy<>());
317	}
318
319	}} // namespaces
320
321	#endif // BOOST_MATH_JACOBI_ELLIPTIC_HPP