math/special_functions/ellint_3.hpp

0001 //  Copyright (c) 2006 Xiaogang Zhang
0002 //  Copyright (c) 2006 John Maddock
0003 //  Use, modification and distribution are subject to the
0004 //  Boost Software License, Version 1.0. (See accompanying file
0005 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006 //
0007 //  History:
0008 //  XZ wrote the original of this file as part of the Google
0009 //  Summer of Code 2006.  JM modified it to fit into the
0010 //  Boost.Math conceptual framework better, and to correctly
0011 //  handle the various corner cases.
0012 //
0013
0014 #ifndef BOOST_MATH_ELLINT_3_HPP
0015 #define BOOST_MATH_ELLINT_3_HPP
0016
0017 #ifdef _MSC_VER
0018 #pragma once
0019 #endif
0020
0021 #include <boost/math/tools/config.hpp>
0022 #include <boost/math/tools/type_traits.hpp>
0023 #include <boost/math/special_functions/math_fwd.hpp>
0024 #include <boost/math/special_functions/ellint_rf.hpp>
0025 #include <boost/math/special_functions/ellint_rj.hpp>
0026 #include <boost/math/special_functions/ellint_1.hpp>
0027 #include <boost/math/special_functions/ellint_2.hpp>
0028 #include <boost/math/special_functions/log1p.hpp>
0029 #include <boost/math/special_functions/atanh.hpp>
0030 #include <boost/math/constants/constants.hpp>
0031 #include <boost/math/policies/error_handling.hpp>
0032 #include <boost/math/tools/workaround.hpp>
0033 #include <boost/math/special_functions/round.hpp>
0034
0035 // Elliptic integrals (complete and incomplete) of the third kind
0036 // Carlson, Numerische Mathematik, vol 33, 1 (1979)
0037
0038 namespace boost { namespace math {
0039
0040 namespace detail{
0041
0042 template <typename T, typename Policy>
0043 BOOST_MATH_CUDA_ENABLED T ellint_pi_imp(T v, T k, T vc, const Policy& pol);
0044
0045 // Elliptic integral (Legendre form) of the third kind
0046 template <typename T, typename Policy>
0047 BOOST_MATH_CUDA_ENABLED T ellint_pi_imp(T v, T phi, T k, T vc, const Policy& pol)
0048 {
0049    // Note vc = 1-v presumably without cancellation error.
0050    BOOST_MATH_STD_USING
0051
0052    constexpr auto function = "boost::math::ellint_3<%1%>(%1%,%1%,%1%)";
0053
0054
0055    T sphi = sin(fabs(phi));
0056    T result = 0;
0057
0058    if (k * k * sphi * sphi > 1)
0059    {
0060       return policies::raise_domain_error<T>(function, "Got k = %1%, function requires |k| <= 1", k, pol);
0061    }
0062    // Special cases first:
0063    if(v == 0)
0064    {
0065       // A&S 17.7.18 & 19
0066       return (k == 0) ? phi : ellint_f_imp(phi, k, pol);
0067    }
0068    if((v > 0) && (1 / v < (sphi * sphi)))
0069    {
0070       // Complex result is a domain error:
0071       return policies::raise_domain_error<T>(function, "Got v = %1%, but result is complex for v > 1 / sin^2(phi)", v, pol);
0072    }
0073
0074    if(v == 1)
0075    {
0076       if (k == 0)
0077          return tan(phi);
0078
0079       // http://functions.wolfram.com/08.06.03.0008.01
0080       T m = k * k;
0081       result = sqrt(1 - m * sphi * sphi) * tan(phi) - ellint_e_imp(phi, k, pol);
0082       result /= 1 - m;
0083       result += ellint_f_imp(phi, k, pol);
0084       return result;
0085    }
0086    if(phi == constants::half_pi<T>())
0087    {
0088       // Have to filter this case out before the next
0089       // special case, otherwise we might get an infinity from
0090       // tan(phi).
0091       // Also note that since we can't represent PI/2 exactly
0092       // in a T, this is a bit of a guess as to the users true
0093       // intent...
0094       //
0095       return ellint_pi_imp(v, k, vc, pol);
0096    }
0097    if((phi > constants::half_pi<T>()) || (phi < 0))
0098    {
0099       // Carlson's algorithm works only for |phi| <= pi/2,
0100       // use the integrand's periodicity to normalize phi
0101       //
0102       // Xiaogang's original code used a cast to long long here
0103       // but that fails if T has more digits than a long long,
0104       // so rewritten to use fmod instead:
0105       //
0106       // See http://functions.wolfram.com/08.06.16.0002.01
0107       //
0108       if(fabs(phi) > 1 / tools::epsilon<T>())
0109       {
0110          // Invalid for v > 1, this case is caught above since v > 1 implies 1/v < sin^2(phi)
0111          BOOST_MATH_ASSERT(v <= 1);
0112          //
0113          // Phi is so large that phi%pi is necessarily zero (or garbage),
0114          // just return the second part of the duplication formula:
0115          //
0116          result = 2 * fabs(phi) * ellint_pi_imp(v, k, vc, pol) / constants::pi<T>();
0117       }
0118       else
0119       {
0120          T rphi = boost::math::tools::fmod_workaround(T(fabs(phi)), T(constants::half_pi<T>()));
0121          T m = boost::math::round((fabs(phi) - rphi) / constants::half_pi<T>());
0122          int sign = 1;
0123          if((m != 0) && (k >= 1))
0124          {
0125             return policies::raise_domain_error<T>(function, "Got k=1 and phi=%1% but the result is complex in that domain", phi, pol);
0126          }
0127          if(boost::math::tools::fmod_workaround(m, T(2)) > T(0.5))
0128          {
0129             m += 1;
0130             sign = -1;
0131             rphi = constants::half_pi<T>() - rphi;
0132          }
0133          result = sign * ellint_pi_imp(v, rphi, k, vc, pol);
0134          if((m > 0) && (vc > 0))
0135             result += m * ellint_pi_imp(v, k, vc, pol);
0136       }
0137       return phi < 0 ? T(-result) : result;
0138    }
0139    if(k == 0)
0140    {
0141       // A&S 17.7.20:
0142       if(v < 1)
0143       {
0144          T vcr = sqrt(vc);
0145          return atan(vcr * tan(phi)) / vcr;
0146       }
0147       else
0148       {
0149          // v > 1:
0150          T vcr = sqrt(-vc);
0151          T arg = vcr * tan(phi);
0152          return (boost::math::log1p(arg, pol) - boost::math::log1p(-arg, pol)) / (2 * vcr);
0153       }
0154    }
0155    if((v < 0) && fabs(k) <= 1)
0156    {
0157       //
0158       // If we don't shift to 0 <= v <= 1 we get
0159       // cancellation errors later on.  Use
0160       // A&S 17.7.15/16 to shift to v > 0.
0161       //
0162       // Mathematica simplifies the expressions
0163       // given in A&S as follows (with thanks to
0164       // Rocco Romeo for figuring these out!):
0165       //
0166       // V = (k2 - n)/(1 - n)
0167       // Assuming[(k2 >= 0 && k2 <= 1) && n < 0, FullSimplify[Sqrt[(1 - V)*(1 - k2 / V)] / Sqrt[((1 - n)*(1 - k2 / n))]]]
0168       // Result: ((-1 + k2) n) / ((-1 + n) (-k2 + n))
0169       //
0170       // Assuming[(k2 >= 0 && k2 <= 1) && n < 0, FullSimplify[k2 / (Sqrt[-n*(k2 - n) / (1 - n)] * Sqrt[(1 - n)*(1 - k2 / n)])]]
0171       // Result : k2 / (k2 - n)
0172       //
0173       // Assuming[(k2 >= 0 && k2 <= 1) && n < 0, FullSimplify[Sqrt[1 / ((1 - n)*(1 - k2 / n))]]]
0174       // Result : Sqrt[n / ((k2 - n) (-1 + n))]
0175       //
0176       T k2 = k * k;
0177       T N = (k2 - v) / (1 - v);
0178       T Nm1 = (1 - k2) / (1 - v);
0179       T p2 = -v * N;
0180       T t;
0181       if (p2 <= tools::min_value<T>())
0182       {
0183          p2 = sqrt(-v) * sqrt(N);
0184       }
0185       else
0186          p2 = sqrt(p2);
0187       T delta = sqrt(1 - k2 * sphi * sphi);
0188       if(N > k2)
0189       {
0190          result = ellint_pi_imp(N, phi, k, Nm1, pol);
0191          result *= v / (v - 1);
0192          result *= (k2 - 1) / (v - k2);
0193       }
0194
0195       if(k != 0)
0196       {
0197          t = ellint_f_imp(phi, k, pol);
0198          t *= k2 / (k2 - v);
0199          result += t;
0200       }
0201       t = v / ((k2 - v) * (v - 1));
0202       if(t > tools::min_value<T>())
0203       {
0204          result += atan((p2 / 2) * sin(2 * phi) / delta) * sqrt(t);
0205       }
0206       else
0207       {
0208          result += atan((p2 / 2) * sin(2 * phi) / delta) * sqrt(fabs(1 / (k2 - v))) * sqrt(fabs(v / (v - 1)));
0209       }
0210       return result;
0211    }
0212    if(k == 1)
0213    {
0214       // See http://functions.wolfram.com/08.06.03.0013.01
0215       result = sqrt(v) * atanh(sqrt(v) * sin(phi), pol) - log(1 / cos(phi) + tan(phi));
0216       result /= v - 1;
0217       return result;
0218    }
0219 #if 0  // disabled but retained for future reference: see below.
0220    if(v > 1)
0221    {
0222       //
0223       // If v > 1 we can use the identity in A&S 17.7.7/8
0224       // to shift to 0 <= v <= 1.  In contrast to previous
0225       // revisions of this header, this identity does now work
0226       // but appears not to produce better error rates in
0227       // practice.  Archived here for future reference...
0228       //
0229       T k2 = k * k;
0230       T N = k2 / v;
0231       T Nm1 = (v - k2) / v;
0232       T p1 = sqrt((-vc) * (1 - k2 / v));
0233       T delta = sqrt(1 - k2 * sphi * sphi);
0234       //
0235       // These next two terms have a large amount of cancellation
0236       // so it's not clear if this relation is useable even if
0237       // the issues with phi > pi/2 can be fixed:
0238       //
0239       result = -ellint_pi_imp(N, phi, k, Nm1, pol);
0240       result += ellint_f_imp(phi, k, pol);
0241       //
0242       // This log term gives the complex result when
0243       //     n > 1/sin^2(phi)
0244       // However that case is dealt with as an error above,
0245       // so we should always get a real result here:
0246       //
0247       result += log((delta + p1 * tan(phi)) / (delta - p1 * tan(phi))) / (2 * p1);
0248       return result;
0249    }
0250 #endif
0251    //
0252    // Carlson's algorithm works only for |phi| <= pi/2,
0253    // by the time we get here phi should already have been
0254    // normalised above.
0255    //
0256    BOOST_MATH_ASSERT(fabs(phi) < constants::half_pi<T>());
0257    BOOST_MATH_ASSERT(phi >= 0);
0258    T x, y, z, p, t;
0259    T cosp = cos(phi);
0260    x = cosp * cosp;
0261    t = sphi * sphi;
0262    y = 1 - k * k * t;
0263    z = 1;
0264    if(v * t < T(0.5))
0265       p = 1 - v * t;
0266    else
0267       p = x + vc * t;
0268    result = sphi * (ellint_rf_imp(x, y, z, pol) + v * t * ellint_rj_imp(x, y, z, p, pol) / 3);
0269
0270    return result;
0271 }
0272
0273 // Complete elliptic integral (Legendre form) of the third kind
0274 template <typename T, typename Policy>
0275 BOOST_MATH_CUDA_ENABLED T ellint_pi_imp(T v, T k, T vc, const Policy& pol)
0276 {
0277     // Note arg vc = 1-v, possibly without cancellation errors
0278     BOOST_MATH_STD_USING
0279     using namespace boost::math::tools;
0280
0281     constexpr auto function = "boost::math::ellint_pi<%1%>(%1%,%1%)";
0282
0283     if (abs(k) >= 1)
0284     {
0285        return policies::raise_domain_error<T>(function, "Got k = %1%, function requires |k| <= 1", k, pol);
0286     }
0287     if(vc <= 0)
0288     {
0289        // Result is complex:
0290        return policies::raise_domain_error<T>(function, "Got v = %1%, function requires v < 1", v, pol);
0291     }
0292
0293     if(v == 0)
0294     {
0295        return (k == 0) ? boost::math::constants::pi<T>() / 2 : boost::math::ellint_1(k, pol);
0296     }
0297
0298     if(v < 0)
0299     {
0300        // Apply A&S 17.7.17:
0301        T k2 = k * k;
0302        T N = (k2 - v) / (1 - v);
0303        T Nm1 = (1 - k2) / (1 - v);
0304        T result = 0;
0305        result = boost::math::detail::ellint_pi_imp(N, k, Nm1, pol);
0306        // This next part is split in two to avoid spurious over/underflow:
0307        result *= -v / (1 - v);
0308        result *= (1 - k2) / (k2 - v);
0309        result += boost::math::ellint_1(k, pol) * k2 / (k2 - v);
0310        return result;
0311     }
0312
0313     T x = 0;
0314     T y = 1 - k * k;
0315     T z = 1;
0316     T p = vc;
0317     T value = ellint_rf_imp(x, y, z, pol) + v * ellint_rj_imp(x, y, z, p, pol) / 3;
0318
0319     return value;
0320 }
0321
0322 template <class T1, class T2, class T3>
0323 BOOST_MATH_CUDA_ENABLED inline typename tools::promote_args<T1, T2, T3>::type ellint_3(T1 k, T2 v, T3 phi, const boost::math::false_type&)
0324 {
0325    return boost::math::ellint_3(k, v, phi, policies::policy<>());
0326 }
0327
0328 template <class T1, class T2, class Policy>
0329 BOOST_MATH_CUDA_ENABLED inline typename tools::promote_args<T1, T2>::type ellint_3(T1 k, T2 v, const Policy& pol, const boost::math::true_type&)
0330 {
0331    typedef typename tools::promote_args<T1, T2>::type result_type;
0332    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0333    return policies::checked_narrowing_cast<result_type, Policy>(
0334       detail::ellint_pi_imp(
0335          static_cast<value_type>(v),
0336          static_cast<value_type>(k),
0337          static_cast<value_type>(1-v),
0338          pol), "boost::math::ellint_3<%1%>(%1%,%1%)");
0339 }
0340
0341 } // namespace detail
0342
0343 template <class T1, class T2, class T3, class Policy>
0344 BOOST_MATH_CUDA_ENABLED inline typename tools::promote_args<T1, T2, T3>::type ellint_3(T1 k, T2 v, T3 phi, const Policy&)
0345 {
0346    typedef typename tools::promote_args<T1, T2, T3>::type result_type;
0347    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0348    typedef typename policies::normalise<Policy, policies::promote_float<false>, policies::promote_double<false> >::type forwarding_policy;
0349    return policies::checked_narrowing_cast<result_type, Policy>(
0350       detail::ellint_pi_imp(
0351          static_cast<value_type>(v),
0352          static_cast<value_type>(phi),
0353          static_cast<value_type>(k),
0354          static_cast<value_type>(1-v),
0355          forwarding_policy()), "boost::math::ellint_3<%1%>(%1%,%1%,%1%)");
0356 }
0357
0358 template <class T1, class T2, class T3>
0359 BOOST_MATH_CUDA_ENABLED typename detail::ellint_3_result<T1, T2, T3>::type ellint_3(T1 k, T2 v, T3 phi)
0360 {
0361    typedef typename policies::is_policy<T3>::type tag_type;
0362    return detail::ellint_3(k, v, phi, tag_type());
0363 }
0364
0365 template <class T1, class T2>
0366 BOOST_MATH_CUDA_ENABLED inline typename tools::promote_args<T1, T2>::type ellint_3(T1 k, T2 v)
0367 {
0368    return ellint_3(k, v, policies::policy<>());
0369 }
0370
0371 }} // namespaces
0372
0373 #endif // BOOST_MATH_ELLINT_3_HPP
0374