math/special_functions/jacobi_elliptic.hpp

0001 // Copyright John Maddock 2012.
0002 // Use, modification and distribution are subject to the
0003 // Boost Software License, Version 1.0.
0004 // (See accompanying file LICENSE_1_0.txt
0005 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_MATH_JACOBI_ELLIPTIC_HPP
0008 #define BOOST_MATH_JACOBI_ELLIPTIC_HPP
0009
0010 #include <boost/math/tools/precision.hpp>
0011 #include <boost/math/tools/promotion.hpp>
0012 #include <boost/math/policies/error_handling.hpp>
0013 #include <boost/math/special_functions/math_fwd.hpp>
0014
0015 namespace boost{ namespace math{
0016
0017 namespace detail{
0018
0019 template <class T, class Policy>
0020 T jacobi_recurse(const T& x, const T& k, T anm1, T bnm1, unsigned N, T* pTn, const Policy& pol)
0021 {
0022    BOOST_MATH_STD_USING
0023    ++N;
0024    T Tn;
0025    T cn = (anm1 - bnm1) / 2;
0026    T an = (anm1 + bnm1) / 2;
0027    if(cn < policies::get_epsilon<T, Policy>())
0028    {
0029       Tn = ldexp(T(1), (int)N) * x * an;
0030    }
0031    else
0032       Tn = jacobi_recurse<T>(x, k, an, sqrt(anm1 * bnm1), N, 0, pol);
0033    if(pTn)
0034       *pTn = Tn;
0035    return (Tn + asin((cn / an) * sin(Tn))) / 2;
0036 }
0037
0038 template <class T, class Policy>
0039 T jacobi_imp(const T& x, const T& k, T* cn, T* dn, const Policy& pol, const char* function)
0040 {
0041    BOOST_MATH_STD_USING
0042    if(k < 0)
0043    {
0044       return *dn = *cn = policies::raise_domain_error<T>(function, "Modulus k must be positive but got %1%.", k, pol);
0045    }
0046    if(k > 1)
0047    {
0048       T xp = x * k;
0049       T kp = 1 / k;
0050       T snp, cnp, dnp;
0051       snp = jacobi_imp(xp, kp, &cnp, &dnp, pol, function);
0052       *cn = dnp;
0053       *dn = cnp;
0054       return snp * kp;
0055    }
0056    //
0057    // Special cases first:
0058    //
0059    if(x == 0)
0060    {
0061       *cn = *dn = 1;
0062       return 0;
0063    }
0064    if(k == 0)
0065    {
0066       *cn = cos(x);
0067       *dn = 1;
0068       return sin(x);
0069    }
0070    if(k == 1)
0071    {
0072       *cn = *dn = 1 / cosh(x);
0073       return tanh(x);
0074    }
0075    //
0076    // Asymptotic forms from A&S 16.13:
0077    //
0078    if(k < tools::forth_root_epsilon<T>())
0079    {
0080       T su = sin(x);
0081       T cu = cos(x);
0082       T m = k * k;
0083       *dn = 1 - m * su * su / 2;
0084       *cn = cu + m * (x - su * cu) * su / 4;
0085       return su - m * (x - su * cu) * cu / 4;
0086    }
0087    /*  Can't get this to work to adequate precision - disabled for now...
0088    //
0089    // Asymptotic forms from A&S 16.15:
0090    //
0091    if(k > 1 - tools::root_epsilon<T>())
0092    {
0093       T tu = tanh(x);
0094       T su = sinh(x);
0095       T cu = cosh(x);
0096       T sec = 1 / cu;
0097       T kp = 1 - k;
0098       T m1 = 2 * kp - kp * kp;
0099       *dn = sec + m1 * (su * cu + x) * tu * sec / 4;
0100       *cn = sec - m1 * (su * cu - x) * tu * sec / 4;
0101       T sn = tu;
0102       T sn2 = m1 * (x * sec * sec - tu) / 4;
0103       T sn3 = (72 * x * cu + 4 * (8 * x * x - 5) * su - 19 * sinh(3 * x) + sinh(5 * x)) * sec * sec * sec * m1 * m1 / 512;
0104       return sn + sn2 - sn3;
0105    }*/
0106    T T1;
0107    T kc = 1 - k;
0108    T k_prime = k < T(0.5) ? T(sqrt(1 - k * k)) : T(sqrt(2 * kc - kc * kc));
0109    T T0 = jacobi_recurse(x, k, T(1), k_prime, 0, &T1, pol);
0110    *cn = cos(T0);
0111    *dn = cos(T0) / cos(T1 - T0);
0112    return sin(T0);
0113 }
0114
0115 } // namespace detail
0116
0117 template <class T, class U, class V, class Policy>
0118 inline typename tools::promote_args<T, U, V>::type jacobi_elliptic(T k, U theta, V* pcn, V* pdn, const Policy&)
0119 {
0120    BOOST_FPU_EXCEPTION_GUARD
0121    typedef typename tools::promote_args<T>::type result_type;
0122    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0123    typedef typename policies::normalise<
0124       Policy,
0125       policies::promote_float<false>,
0126       policies::promote_double<false>,
0127       policies::discrete_quantile<>,
0128       policies::assert_undefined<> >::type forwarding_policy;
0129
0130    static const char* function = "boost::math::jacobi_elliptic<%1%>(%1%)";
0131
0132    value_type sn, cn, dn;
0133    sn = detail::jacobi_imp<value_type>(static_cast<value_type>(theta), static_cast<value_type>(k), &cn, &dn, forwarding_policy(), function);
0134    if(pcn)
0135       *pcn = policies::checked_narrowing_cast<result_type, Policy>(cn, function);
0136    if(pdn)
0137       *pdn = policies::checked_narrowing_cast<result_type, Policy>(dn, function);
0138    return policies::checked_narrowing_cast<result_type, Policy>(sn, function);
0139 }
0140
0141 template <class T, class U, class V>
0142 inline typename tools::promote_args<T, U, V>::type jacobi_elliptic(T k, U theta, V* pcn, V* pdn)
0143 {
0144    return jacobi_elliptic(k, theta, pcn, pdn, policies::policy<>());
0145 }
0146
0147 template <class U, class T, class Policy>
0148 inline typename tools::promote_args<T, U>::type jacobi_sn(U k, T theta, const Policy& pol)
0149 {
0150    typedef typename tools::promote_args<T, U>::type result_type;
0151    return jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), static_cast<result_type*>(nullptr), static_cast<result_type*>(nullptr), pol);
0152 }
0153
0154 template <class U, class T>
0155 inline typename tools::promote_args<T, U>::type jacobi_sn(U k, T theta)
0156 {
0157    return jacobi_sn(k, theta, policies::policy<>());
0158 }
0159
0160 template <class T, class U, class Policy>
0161 inline typename tools::promote_args<T, U>::type jacobi_cn(T k, U theta, const Policy& pol)
0162 {
0163    typedef typename tools::promote_args<T, U>::type result_type;
0164    result_type cn;
0165    jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), &cn, static_cast<result_type*>(nullptr), pol);
0166    return cn;
0167 }
0168
0169 template <class T, class U>
0170 inline typename tools::promote_args<T, U>::type jacobi_cn(T k, U theta)
0171 {
0172    return jacobi_cn(k, theta, policies::policy<>());
0173 }
0174
0175 template <class T, class U, class Policy>
0176 inline typename tools::promote_args<T, U>::type jacobi_dn(T k, U theta, const Policy& pol)
0177 {
0178    typedef typename tools::promote_args<T, U>::type result_type;
0179    result_type dn;
0180    jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), static_cast<result_type*>(nullptr), &dn, pol);
0181    return dn;
0182 }
0183
0184 template <class T, class U>
0185 inline typename tools::promote_args<T, U>::type jacobi_dn(T k, U theta)
0186 {
0187    return jacobi_dn(k, theta, policies::policy<>());
0188 }
0189
0190 template <class T, class U, class Policy>
0191 inline typename tools::promote_args<T, U>::type jacobi_cd(T k, U theta, const Policy& pol)
0192 {
0193    typedef typename tools::promote_args<T, U>::type result_type;
0194    result_type cn, dn;
0195    jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), &cn, &dn, pol);
0196    return cn / dn;
0197 }
0198
0199 template <class T, class U>
0200 inline typename tools::promote_args<T, U>::type jacobi_cd(T k, U theta)
0201 {
0202    return jacobi_cd(k, theta, policies::policy<>());
0203 }
0204
0205 template <class T, class U, class Policy>
0206 inline typename tools::promote_args<T, U>::type jacobi_dc(T k, U theta, const Policy& pol)
0207 {
0208    typedef typename tools::promote_args<T, U>::type result_type;
0209    result_type cn, dn;
0210    jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), &cn, &dn, pol);
0211    return dn / cn;
0212 }
0213
0214 template <class T, class U>
0215 inline typename tools::promote_args<T, U>::type jacobi_dc(T k, U theta)
0216 {
0217    return jacobi_dc(k, theta, policies::policy<>());
0218 }
0219
0220 template <class T, class U, class Policy>
0221 inline typename tools::promote_args<T, U>::type jacobi_ns(T k, U theta, const Policy& pol)
0222 {
0223    typedef typename tools::promote_args<T, U>::type result_type;
0224    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);
0225 }
0226
0227 template <class T, class U>
0228 inline typename tools::promote_args<T, U>::type jacobi_ns(T k, U theta)
0229 {
0230    return jacobi_ns(k, theta, policies::policy<>());
0231 }
0232
0233 template <class T, class U, class Policy>
0234 inline typename tools::promote_args<T, U>::type jacobi_sd(T k, U theta, const Policy& pol)
0235 {
0236    typedef typename tools::promote_args<T, U>::type result_type;
0237    result_type sn, dn;
0238    sn = jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), static_cast<result_type*>(nullptr), &dn, pol);
0239    return sn / dn;
0240 }
0241
0242 template <class T, class U>
0243 inline typename tools::promote_args<T, U>::type jacobi_sd(T k, U theta)
0244 {
0245    return jacobi_sd(k, theta, policies::policy<>());
0246 }
0247
0248 template <class T, class U, class Policy>
0249 inline typename tools::promote_args<T, U>::type jacobi_ds(T k, U theta, const Policy& pol)
0250 {
0251    typedef typename tools::promote_args<T, U>::type result_type;
0252    result_type sn, dn;
0253    sn = jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), static_cast<result_type*>(nullptr), &dn, pol);
0254    return dn / sn;
0255 }
0256
0257 template <class T, class U>
0258 inline typename tools::promote_args<T, U>::type jacobi_ds(T k, U theta)
0259 {
0260    return jacobi_ds(k, theta, policies::policy<>());
0261 }
0262
0263 template <class T, class U, class Policy>
0264 inline typename tools::promote_args<T, U>::type jacobi_nc(T k, U theta, const Policy& pol)
0265 {
0266    return 1 / jacobi_cn(k, theta, pol);
0267 }
0268
0269 template <class T, class U>
0270 inline typename tools::promote_args<T, U>::type jacobi_nc(T k, U theta)
0271 {
0272    return jacobi_nc(k, theta, policies::policy<>());
0273 }
0274
0275 template <class T, class U, class Policy>
0276 inline typename tools::promote_args<T, U>::type jacobi_nd(T k, U theta, const Policy& pol)
0277 {
0278    return 1 / jacobi_dn(k, theta, pol);
0279 }
0280
0281 template <class T, class U>
0282 inline typename tools::promote_args<T, U>::type jacobi_nd(T k, U theta)
0283 {
0284    return jacobi_nd(k, theta, policies::policy<>());
0285 }
0286
0287 template <class T, class U, class Policy>
0288 inline typename tools::promote_args<T, U>::type jacobi_sc(T k, U theta, const Policy& pol)
0289 {
0290    typedef typename tools::promote_args<T, U>::type result_type;
0291    result_type sn, cn;
0292    sn = jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), &cn, static_cast<result_type*>(nullptr), pol);
0293    return sn / cn;
0294 }
0295
0296 template <class T, class U>
0297 inline typename tools::promote_args<T, U>::type jacobi_sc(T k, U theta)
0298 {
0299    return jacobi_sc(k, theta, policies::policy<>());
0300 }
0301
0302 template <class T, class U, class Policy>
0303 inline typename tools::promote_args<T, U>::type jacobi_cs(T k, U theta, const Policy& pol)
0304 {
0305    typedef typename tools::promote_args<T, U>::type result_type;
0306    result_type sn, cn;
0307    sn = jacobi_elliptic(static_cast<result_type>(k), static_cast<result_type>(theta), &cn, static_cast<result_type*>(nullptr), pol);
0308    return cn / sn;
0309 }
0310
0311 template <class T, class U>
0312 inline typename tools::promote_args<T, U>::type jacobi_cs(T k, U theta)
0313 {
0314    return jacobi_cs(k, theta, policies::policy<>());
0315 }
0316
0317 }} // namespaces
0318
0319 #endif // BOOST_MATH_JACOBI_ELLIPTIC_HPP