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