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