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)
7 #ifndef BOOST_MATH_JACOBI_ELLIPTIC_HPP
8 #define BOOST_MATH_JACOBI_ELLIPTIC_HPP
10 #include <boost/math/tools/precision.hpp>
11 #include <boost/math/tools/promotion.hpp>
12 #include <boost/math/policies/error_handling.hpp>
14 namespace boost{ namespace math{
18 template <class T, class Policy>
19 T jacobi_recurse(const T& x, const T& k, T anm1, T bnm1, unsigned N, T* pTn, const Policy& pol)
24 T cn = (anm1 - bnm1) / 2;
25 T an = (anm1 + bnm1) / 2;
26 if(cn < policies::get_epsilon<T, Policy>())
28 Tn = ldexp(T(1), (int)N) * x * an;
31 Tn = jacobi_recurse<T>(x, k, an, sqrt(anm1 * bnm1), N, 0, pol);
34 return (Tn + asin((cn / an) * sin(Tn))) / 2;
37 template <class T, class Policy>
38 T jacobi_imp(const T& x, const T& k, T* cn, T* dn, const Policy& pol, const char* function)
43 *cn = policies::raise_domain_error<T>(function, "Modulus k must be positive but got %1%.", k, pol);
52 snp = jacobi_imp(xp, kp, &cnp, &dnp, pol, function);
58 // Special cases first:
73 *cn = *dn = 1 / cosh(x);
77 // Asymptotic forms from A&S 16.13:
79 if(k < tools::forth_root_epsilon<T>())
84 *dn = 1 - m * su * su / 2;
85 *cn = cu + m * (x - su * cu) * su / 4;
86 return su - m * (x - su * cu) * cu / 4;
88 /* Can't get this to work to adequate precision - disabled for now...
90 // Asymptotic forms from A&S 16.15:
92 if(k > 1 - tools::root_epsilon<T>())
99 T m1 = 2 * kp - kp * kp;
100 *dn = sec + m1 * (su * cu + x) * tu * sec / 4;
101 *cn = sec - m1 * (su * cu - x) * tu * sec / 4;
103 T sn2 = m1 * (x * sec * sec - tu) / 4;
104 T sn3 = (72 * x * cu + 4 * (8 * x * x - 5) * su - 19 * sinh(3 * x) + sinh(5 * x)) * sec * sec * sec * m1 * m1 / 512;
105 return sn + sn2 - sn3;
109 T k_prime = k < 0.5 ? T(sqrt(1 - k * k)) : T(sqrt(2 * kc - kc * kc));
110 T T0 = jacobi_recurse(x, k, T(1), k_prime, 0, &T1, pol);
112 *dn = cos(T0) / cos(T1 - T0);
116 } // namespace detail
118 template <class T, class U, class V, class Policy>
119 inline typename tools::promote_args<T, U, V>::type jacobi_elliptic(T k, U theta, V* pcn, V* pdn, const Policy&)
121 BOOST_FPU_EXCEPTION_GUARD
122 typedef typename tools::promote_args<T>::type result_type;
123 typedef typename policies::evaluation<result_type, Policy>::type value_type;
124 typedef typename policies::normalise<
126 policies::promote_float<false>,
127 policies::promote_double<false>,
128 policies::discrete_quantile<>,
129 policies::assert_undefined<> >::type forwarding_policy;
131 static const char* function = "boost::math::jacobi_elliptic<%1%>(%1%)";
133 value_type sn, cn, dn;
134 sn = detail::jacobi_imp<value_type>(static_cast<value_type>(theta), static_cast<value_type>(k), &cn, &dn, forwarding_policy(), function);
136 *pcn = policies::checked_narrowing_cast<result_type, Policy>(cn, function);
138 *pdn = policies::checked_narrowing_cast<result_type, Policy>(dn, function);
139 return policies::checked_narrowing_cast<result_type, Policy>(sn, function);;
142 template <class T, class U, class V>
143 inline typename tools::promote_args<T, U, V>::type jacobi_elliptic(T k, U theta, V* pcn, V* pdn)
145 return jacobi_elliptic(k, theta, pcn, pdn, policies::policy<>());
148 template <class U, class T, class Policy>
149 inline typename tools::promote_args<T, U>::type jacobi_sn(U k, T theta, const Policy& pol)
151 typedef typename tools::promote_args<T, U>::type result_type;
152 return jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), static_cast<result_type*>(0), static_cast<result_type*>(0), pol);
155 template <class U, class T>
156 inline typename tools::promote_args<T, U>::type jacobi_sn(U k, T theta)
158 return jacobi_sn(k, theta, policies::policy<>());
161 template <class T, class U, class Policy>
162 inline typename tools::promote_args<T, U>::type jacobi_cn(T k, U theta, const Policy& pol)
164 typedef typename tools::promote_args<T, U>::type result_type;
166 jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), &cn, static_cast<result_type*>(0), pol);
170 template <class T, class U>
171 inline typename tools::promote_args<T, U>::type jacobi_cn(T k, U theta)
173 return jacobi_cn(k, theta, policies::policy<>());
176 template <class T, class U, class Policy>
177 inline typename tools::promote_args<T, U>::type jacobi_dn(T k, U theta, const Policy& pol)
179 typedef typename tools::promote_args<T, U>::type result_type;
181 jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), static_cast<result_type*>(0), &dn, pol);
185 template <class T, class U>
186 inline typename tools::promote_args<T, U>::type jacobi_dn(T k, U theta)
188 return jacobi_dn(k, theta, policies::policy<>());
191 template <class T, class U, class Policy>
192 inline typename tools::promote_args<T, U>::type jacobi_cd(T k, U theta, const Policy& pol)
194 typedef typename tools::promote_args<T, U>::type result_type;
196 jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), &cn, &dn, pol);
200 template <class T, class U>
201 inline typename tools::promote_args<T, U>::type jacobi_cd(T k, U theta)
203 return jacobi_cd(k, theta, policies::policy<>());
206 template <class T, class U, class Policy>
207 inline typename tools::promote_args<T, U>::type jacobi_dc(T k, U theta, const Policy& pol)
209 typedef typename tools::promote_args<T, U>::type result_type;
211 jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), &cn, &dn, pol);
215 template <class T, class U>
216 inline typename tools::promote_args<T, U>::type jacobi_dc(T k, U theta)
218 return jacobi_dc(k, theta, policies::policy<>());
221 template <class T, class U, class Policy>
222 inline typename tools::promote_args<T, U>::type jacobi_ns(T k, U theta, const Policy& pol)
224 typedef typename tools::promote_args<T, U>::type result_type;
225 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);
228 template <class T, class U>
229 inline typename tools::promote_args<T, U>::type jacobi_ns(T k, U theta)
231 return jacobi_ns(k, theta, policies::policy<>());
234 template <class T, class U, class Policy>
235 inline typename tools::promote_args<T, U>::type jacobi_sd(T k, U theta, const Policy& pol)
237 typedef typename tools::promote_args<T, U>::type result_type;
239 sn = jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), static_cast<result_type*>(0), &dn, pol);
243 template <class T, class U>
244 inline typename tools::promote_args<T, U>::type jacobi_sd(T k, U theta)
246 return jacobi_sd(k, theta, policies::policy<>());
249 template <class T, class U, class Policy>
250 inline typename tools::promote_args<T, U>::type jacobi_ds(T k, U theta, const Policy& pol)
252 typedef typename tools::promote_args<T, U>::type result_type;
254 sn = jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), static_cast<result_type*>(0), &dn, pol);
258 template <class T, class U>
259 inline typename tools::promote_args<T, U>::type jacobi_ds(T k, U theta)
261 return jacobi_ds(k, theta, policies::policy<>());
264 template <class T, class U, class Policy>
265 inline typename tools::promote_args<T, U>::type jacobi_nc(T k, U theta, const Policy& pol)
267 return 1 / jacobi_cn(k, theta, pol);
270 template <class T, class U>
271 inline typename tools::promote_args<T, U>::type jacobi_nc(T k, U theta)
273 return jacobi_nc(k, theta, policies::policy<>());
276 template <class T, class U, class Policy>
277 inline typename tools::promote_args<T, U>::type jacobi_nd(T k, U theta, const Policy& pol)
279 return 1 / jacobi_dn(k, theta, pol);
282 template <class T, class U>
283 inline typename tools::promote_args<T, U>::type jacobi_nd(T k, U theta)
285 return jacobi_nd(k, theta, policies::policy<>());
288 template <class T, class U, class Policy>
289 inline typename tools::promote_args<T, U>::type jacobi_sc(T k, U theta, const Policy& pol)
291 typedef typename tools::promote_args<T, U>::type result_type;
293 sn = jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), &cn, static_cast<result_type*>(0), pol);
297 template <class T, class U>
298 inline typename tools::promote_args<T, U>::type jacobi_sc(T k, U theta)
300 return jacobi_sc(k, theta, policies::policy<>());
303 template <class T, class U, class Policy>
304 inline typename tools::promote_args<T, U>::type jacobi_cs(T k, U theta, const Policy& pol)
306 typedef typename tools::promote_args<T, U>::type result_type;
308 sn = jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), &cn, static_cast<result_type*>(0), pol);
312 template <class T, class U>
313 inline typename tools::promote_args<T, U>::type jacobi_cs(T k, U theta)
315 return jacobi_cs(k, theta, policies::policy<>());
320 #endif // BOOST_MATH_JACOBI_ELLIPTIC_HPP