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>
jacobi_elliptic(T k,U theta,V * pcn,V * pdn,const 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>
jacobi_elliptic(T k,U theta,V * pcn,V * pdn)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>
jacobi_sn(U k,T theta,const Policy & pol)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>
jacobi_sn(U k,T theta)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>
jacobi_cn(T k,U theta,const Policy & pol)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>
jacobi_cn(T k,U theta)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>
jacobi_dn(T k,U theta,const Policy & pol)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>
jacobi_dn(T k,U theta)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>
jacobi_cd(T k,U theta,const Policy & pol)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>
jacobi_cd(T k,U theta)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>
jacobi_dc(T k,U theta,const Policy & pol)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>
jacobi_dc(T k,U theta)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>
jacobi_ns(T k,U theta,const Policy & pol)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>
jacobi_ns(T k,U theta)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>
jacobi_sd(T k,U theta,const Policy & pol)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>
jacobi_sd(T k,U theta)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>
jacobi_ds(T k,U theta,const Policy & pol)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>
jacobi_ds(T k,U theta)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>
jacobi_nc(T k,U theta,const Policy & pol)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>
jacobi_nc(T k,U theta)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>
jacobi_nd(T k,U theta,const Policy & pol)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>
jacobi_nd(T k,U theta)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>
jacobi_sc(T k,U theta,const Policy & pol)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>
jacobi_sc(T k,U theta)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>
jacobi_cs(T k,U theta,const Policy & pol)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>
jacobi_cs(T k,U theta)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
322