// Copyright John Maddock 2012. // Use, modification and distribution are subject to the // Boost Software License, Version 1.0. // (See accompanying file LICENSE_1_0.txt // or copy at http://www.boost.org/LICENSE_1_0.txt) #ifndef BOOST_MATH_JACOBI_ELLIPTIC_HPP #define BOOST_MATH_JACOBI_ELLIPTIC_HPP #include #include #include #include namespace boost{ namespace math{ namespace detail{ template T jacobi_recurse(const T& x, const T& k, T anm1, T bnm1, unsigned N, T* pTn, const Policy& pol) { BOOST_MATH_STD_USING ++N; T Tn; T cn = (anm1 - bnm1) / 2; T an = (anm1 + bnm1) / 2; if(cn < policies::get_epsilon()) { Tn = ldexp(T(1), (int)N) * x * an; } else Tn = jacobi_recurse(x, k, an, sqrt(anm1 * bnm1), N, 0, pol); if(pTn) *pTn = Tn; return (Tn + asin((cn / an) * sin(Tn))) / 2; } template T jacobi_imp(const T& x, const T& k, T* cn, T* dn, const Policy& pol, const char* function) { BOOST_MATH_STD_USING if(k < 0) { *cn = policies::raise_domain_error(function, "Modulus k must be positive but got %1%.", k, pol); *dn = *cn; return *cn; } if(k > 1) { T xp = x * k; T kp = 1 / k; T snp, cnp, dnp; snp = jacobi_imp(xp, kp, &cnp, &dnp, pol, function); *cn = dnp; *dn = cnp; return snp * kp; } // // Special cases first: // if(x == 0) { *cn = *dn = 1; return 0; } if(k == 0) { *cn = cos(x); *dn = 1; return sin(x); } if(k == 1) { *cn = *dn = 1 / cosh(x); return tanh(x); } // // Asymptotic forms from A&S 16.13: // if(k < tools::forth_root_epsilon()) { T su = sin(x); T cu = cos(x); T m = k * k; *dn = 1 - m * su * su / 2; *cn = cu + m * (x - su * cu) * su / 4; return su - m * (x - su * cu) * cu / 4; } /* Can't get this to work to adequate precision - disabled for now... // // Asymptotic forms from A&S 16.15: // if(k > 1 - tools::root_epsilon()) { T tu = tanh(x); T su = sinh(x); T cu = cosh(x); T sec = 1 / cu; T kp = 1 - k; T m1 = 2 * kp - kp * kp; *dn = sec + m1 * (su * cu + x) * tu * sec / 4; *cn = sec - m1 * (su * cu - x) * tu * sec / 4; T sn = tu; T sn2 = m1 * (x * sec * sec - tu) / 4; T sn3 = (72 * x * cu + 4 * (8 * x * x - 5) * su - 19 * sinh(3 * x) + sinh(5 * x)) * sec * sec * sec * m1 * m1 / 512; return sn + sn2 - sn3; }*/ T T1; T kc = 1 - k; T k_prime = k < 0.5 ? T(sqrt(1 - k * k)) : T(sqrt(2 * kc - kc * kc)); T T0 = jacobi_recurse(x, k, T(1), k_prime, 0, &T1, pol); *cn = cos(T0); *dn = cos(T0) / cos(T1 - T0); return sin(T0); } } // namespace detail template inline typename tools::promote_args::type jacobi_elliptic(T k, U theta, V* pcn, V* pdn, const Policy&) { BOOST_FPU_EXCEPTION_GUARD typedef typename tools::promote_args::type result_type; typedef typename policies::evaluation::type value_type; typedef typename policies::normalise< Policy, policies::promote_float, policies::promote_double, policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; static const char* function = "boost::math::jacobi_elliptic<%1%>(%1%)"; value_type sn, cn, dn; sn = detail::jacobi_imp(static_cast(theta), static_cast(k), &cn, &dn, forwarding_policy(), function); if(pcn) *pcn = policies::checked_narrowing_cast(cn, function); if(pdn) *pdn = policies::checked_narrowing_cast(dn, function); return policies::checked_narrowing_cast(sn, function);; } template inline typename tools::promote_args::type jacobi_elliptic(T k, U theta, V* pcn, V* pdn) { return jacobi_elliptic(k, theta, pcn, pdn, policies::policy<>()); } template inline typename tools::promote_args::type jacobi_sn(U k, T theta, const Policy& pol) { typedef typename tools::promote_args::type result_type; return jacobi_elliptic(static_cast(k), static_cast(theta), static_cast(0), static_cast(0), pol); } template inline typename tools::promote_args::type jacobi_sn(U k, T theta) { return jacobi_sn(k, theta, policies::policy<>()); } template inline typename tools::promote_args::type jacobi_cn(T k, U theta, const Policy& pol) { typedef typename tools::promote_args::type result_type; result_type cn; jacobi_elliptic(static_cast(k), static_cast(theta), &cn, static_cast(0), pol); return cn; } template inline typename tools::promote_args::type jacobi_cn(T k, U theta) { return jacobi_cn(k, theta, policies::policy<>()); } template inline typename tools::promote_args::type jacobi_dn(T k, U theta, const Policy& pol) { typedef typename tools::promote_args::type result_type; result_type dn; jacobi_elliptic(static_cast(k), static_cast(theta), static_cast(0), &dn, pol); return dn; } template inline typename tools::promote_args::type jacobi_dn(T k, U theta) { return jacobi_dn(k, theta, policies::policy<>()); } template inline typename tools::promote_args::type jacobi_cd(T k, U theta, const Policy& pol) { typedef typename tools::promote_args::type result_type; result_type cn, dn; jacobi_elliptic(static_cast(k), static_cast(theta), &cn, &dn, pol); return cn / dn; } template inline typename tools::promote_args::type jacobi_cd(T k, U theta) { return jacobi_cd(k, theta, policies::policy<>()); } template inline typename tools::promote_args::type jacobi_dc(T k, U theta, const Policy& pol) { typedef typename tools::promote_args::type result_type; result_type cn, dn; jacobi_elliptic(static_cast(k), static_cast(theta), &cn, &dn, pol); return dn / cn; } template inline typename tools::promote_args::type jacobi_dc(T k, U theta) { return jacobi_dc(k, theta, policies::policy<>()); } template inline typename tools::promote_args::type jacobi_ns(T k, U theta, const Policy& pol) { typedef typename tools::promote_args::type result_type; return 1 / jacobi_elliptic(static_cast(k), static_cast(theta), static_cast(0), static_cast(0), pol); } template inline typename tools::promote_args::type jacobi_ns(T k, U theta) { return jacobi_ns(k, theta, policies::policy<>()); } template inline typename tools::promote_args::type jacobi_sd(T k, U theta, const Policy& pol) { typedef typename tools::promote_args::type result_type; result_type sn, dn; sn = jacobi_elliptic(static_cast(k), static_cast(theta), static_cast(0), &dn, pol); return sn / dn; } template inline typename tools::promote_args::type jacobi_sd(T k, U theta) { return jacobi_sd(k, theta, policies::policy<>()); } template inline typename tools::promote_args::type jacobi_ds(T k, U theta, const Policy& pol) { typedef typename tools::promote_args::type result_type; result_type sn, dn; sn = jacobi_elliptic(static_cast(k), static_cast(theta), static_cast(0), &dn, pol); return dn / sn; } template inline typename tools::promote_args::type jacobi_ds(T k, U theta) { return jacobi_ds(k, theta, policies::policy<>()); } template inline typename tools::promote_args::type jacobi_nc(T k, U theta, const Policy& pol) { return 1 / jacobi_cn(k, theta, pol); } template inline typename tools::promote_args::type jacobi_nc(T k, U theta) { return jacobi_nc(k, theta, policies::policy<>()); } template inline typename tools::promote_args::type jacobi_nd(T k, U theta, const Policy& pol) { return 1 / jacobi_dn(k, theta, pol); } template inline typename tools::promote_args::type jacobi_nd(T k, U theta) { return jacobi_nd(k, theta, policies::policy<>()); } template inline typename tools::promote_args::type jacobi_sc(T k, U theta, const Policy& pol) { typedef typename tools::promote_args::type result_type; result_type sn, cn; sn = jacobi_elliptic(static_cast(k), static_cast(theta), &cn, static_cast(0), pol); return sn / cn; } template inline typename tools::promote_args::type jacobi_sc(T k, U theta) { return jacobi_sc(k, theta, policies::policy<>()); } template inline typename tools::promote_args::type jacobi_cs(T k, U theta, const Policy& pol) { typedef typename tools::promote_args::type result_type; result_type sn, cn; sn = jacobi_elliptic(static_cast(k), static_cast(theta), &cn, static_cast(0), pol); return cn / sn; } template inline typename tools::promote_args::type jacobi_cs(T k, U theta) { return jacobi_cs(k, theta, policies::policy<>()); } }} // namespaces #endif // BOOST_MATH_JACOBI_ELLIPTIC_HPP