boost/math/special_functions/jacobi_elliptic.hpp

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