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/special_functions/math_fwd.hpp>
0022 #include <boost/math/special_functions/ellint_rf.hpp>
0023 #include <boost/math/special_functions/ellint_rj.hpp>
0024 #include <boost/math/special_functions/ellint_1.hpp>
0025 #include <boost/math/special_functions/ellint_2.hpp>
0026 #include <boost/math/special_functions/log1p.hpp>
0027 #include <boost/math/special_functions/atanh.hpp>
0028 #include <boost/math/constants/constants.hpp>
0029 #include <boost/math/policies/error_handling.hpp>
0030 #include <boost/math/tools/workaround.hpp>
0031 #include <boost/math/special_functions/round.hpp>
0032
0033 // Elliptic integrals (complete and incomplete) of the third kind
0034 // Carlson, Numerische Mathematik, vol 33, 1 (1979)
0035
0036 namespace boost { namespace math {
0037
0038 namespace detail{
0039
0040 template <typename T, typename Policy>
0041 T ellint_pi_imp(T v, T k, T vc, const Policy& pol);
0042
0043 // Elliptic integral (Legendre form) of the third kind
0044 template <typename T, typename Policy>
0045 T ellint_pi_imp(T v, T phi, T k, T vc, const Policy& pol)
0046 {
0047    // Note vc = 1-v presumably without cancellation error.
0048    BOOST_MATH_STD_USING
0049
0050    static const char* function = "boost::math::ellint_3<%1%>(%1%,%1%,%1%)";
0051
0052
0053    T sphi = sin(fabs(phi));
0054    T result = 0;
0055
0056    if (k * k * sphi * sphi > 1)
0057    {
0058       return policies::raise_domain_error<T>(function,
0059          "Got k = %1%, function requires |k| <= 1", k, pol);
0060    }
0061    // Special cases first:
0062    if(v == 0)
0063    {
0064       // A&S 17.7.18 & 19
0065       return (k == 0) ? phi : ellint_f_imp(phi, k, pol);
0066    }
0067    if((v > 0) && (1 / v < (sphi * sphi)))
0068    {
0069       // Complex result is a domain error:
0070       return policies::raise_domain_error<T>(function,
0071          "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          if(v > 1)
0111             return policies::raise_domain_error<T>(
0112             function,
0113             "Got v = %1%, but this is only supported for 0 <= phi <= pi/2", v, pol);
0114          //
0115          // Phi is so large that phi%pi is necessarily zero (or garbage),
0116          // just return the second part of the duplication formula:
0117          //
0118          result = 2 * fabs(phi) * ellint_pi_imp(v, k, vc, pol) / constants::pi<T>();
0119       }
0120       else
0121       {
0122          T rphi = boost::math::tools::fmod_workaround(T(fabs(phi)), T(constants::half_pi<T>()));
0123          T m = boost::math::round((fabs(phi) - rphi) / constants::half_pi<T>());
0124          int sign = 1;
0125          if((m != 0) && (k >= 1))
0126          {
0127             return policies::raise_domain_error<T>(function, "Got k=1 and phi=%1% but the result is complex in that domain", phi, pol);
0128          }
0129          if(boost::math::tools::fmod_workaround(m, T(2)) > T(0.5))
0130          {
0131             m += 1;
0132             sign = -1;
0133             rphi = constants::half_pi<T>() - rphi;
0134          }
0135          result = sign * ellint_pi_imp(v, rphi, k, vc, pol);
0136          if((m > 0) && (vc > 0))
0137             result += m * ellint_pi_imp(v, k, vc, pol);
0138       }
0139       return phi < 0 ? T(-result) : result;
0140    }
0141    if(k == 0)
0142    {
0143       // A&S 17.7.20:
0144       if(v < 1)
0145       {
0146          T vcr = sqrt(vc);
0147          return atan(vcr * tan(phi)) / vcr;
0148       }
0149       else
0150       {
0151          // v > 1:
0152          T vcr = sqrt(-vc);
0153          T arg = vcr * tan(phi);
0154          return (boost::math::log1p(arg, pol) - boost::math::log1p(-arg, pol)) / (2 * vcr);
0155       }
0156    }
0157    if((v < 0) && fabs(k) <= 1)
0158    {
0159       //
0160       // If we don't shift to 0 <= v <= 1 we get
0161       // cancellation errors later on.  Use
0162       // A&S 17.7.15/16 to shift to v > 0.
0163       //
0164       // Mathematica simplifies the expressions
0165       // given in A&S as follows (with thanks to
0166       // Rocco Romeo for figuring these out!):
0167       //
0168       // V = (k2 - n)/(1 - n)
0169       // Assuming[(k2 >= 0 && k2 <= 1) && n < 0, FullSimplify[Sqrt[(1 - V)*(1 - k2 / V)] / Sqrt[((1 - n)*(1 - k2 / n))]]]
0170       // Result: ((-1 + k2) n) / ((-1 + n) (-k2 + n))
0171       //
0172       // Assuming[(k2 >= 0 && k2 <= 1) && n < 0, FullSimplify[k2 / (Sqrt[-n*(k2 - n) / (1 - n)] * Sqrt[(1 - n)*(1 - k2 / n)])]]
0173       // Result : k2 / (k2 - n)
0174       //
0175       // Assuming[(k2 >= 0 && k2 <= 1) && n < 0, FullSimplify[Sqrt[1 / ((1 - n)*(1 - k2 / n))]]]
0176       // Result : Sqrt[n / ((k2 - n) (-1 + n))]
0177       //
0178       T k2 = k * k;
0179       T N = (k2 - v) / (1 - v);
0180       T Nm1 = (1 - k2) / (1 - v);
0181       T p2 = -v * N;
0182       T t;
0183       if(p2 <= tools::min_value<T>())
0184          p2 = sqrt(-v) * sqrt(N);
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 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     static const char* function = "boost::math::ellint_pi<%1%>(%1%,%1%)";
0282
0283     if (abs(k) >= 1)
0284     {
0285        return policies::raise_domain_error<T>(function,
0286             "Got k = %1%, function requires |k| <= 1", k, pol);
0287     }
0288     if(vc <= 0)
0289     {
0290        // Result is complex:
0291        return policies::raise_domain_error<T>(function,
0292             "Got v = %1%, function requires v < 1", v, pol);
0293     }
0294
0295     if(v == 0)
0296     {
0297        return (k == 0) ? boost::math::constants::pi<T>() / 2 : boost::math::ellint_1(k, pol);
0298     }
0299
0300     if(v < 0)
0301     {
0302        // Apply A&S 17.7.17:
0303        T k2 = k * k;
0304        T N = (k2 - v) / (1 - v);
0305        T Nm1 = (1 - k2) / (1 - v);
0306        T result = 0;
0307        result = boost::math::detail::ellint_pi_imp(N, k, Nm1, pol);
0308        // This next part is split in two to avoid spurious over/underflow:
0309        result *= -v / (1 - v);
0310        result *= (1 - k2) / (k2 - v);
0311        result += boost::math::ellint_1(k, pol) * k2 / (k2 - v);
0312        return result;
0313     }
0314
0315     T x = 0;
0316     T y = 1 - k * k;
0317     T z = 1;
0318     T p = vc;
0319     T value = ellint_rf_imp(x, y, z, pol) + v * ellint_rj_imp(x, y, z, p, pol) / 3;
0320
0321     return value;
0322 }
0323
0324 template <class T1, class T2, class T3>
0325 inline typename tools::promote_args<T1, T2, T3>::type ellint_3(T1 k, T2 v, T3 phi, const std::false_type&)
0326 {
0327    return boost::math::ellint_3(k, v, phi, policies::policy<>());
0328 }
0329
0330 template <class T1, class T2, class Policy>
0331 inline typename tools::promote_args<T1, T2>::type ellint_3(T1 k, T2 v, const Policy& pol, const std::true_type&)
0332 {
0333    typedef typename tools::promote_args<T1, T2>::type result_type;
0334    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0335    return policies::checked_narrowing_cast<result_type, Policy>(
0336       detail::ellint_pi_imp(
0337          static_cast<value_type>(v),
0338          static_cast<value_type>(k),
0339          static_cast<value_type>(1-v),
0340          pol), "boost::math::ellint_3<%1%>(%1%,%1%)");
0341 }
0342
0343 } // namespace detail
0344
0345 template <class T1, class T2, class T3, class Policy>
0346 inline typename tools::promote_args<T1, T2, T3>::type ellint_3(T1 k, T2 v, T3 phi, const Policy&)
0347 {
0348    typedef typename tools::promote_args<T1, T2, T3>::type result_type;
0349    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0350    typedef typename policies::normalise<Policy, policies::promote_float<false>, policies::promote_double<false> >::type forwarding_policy;
0351    return policies::checked_narrowing_cast<result_type, Policy>(
0352       detail::ellint_pi_imp(
0353          static_cast<value_type>(v),
0354          static_cast<value_type>(phi),
0355          static_cast<value_type>(k),
0356          static_cast<value_type>(1-v),
0357          forwarding_policy()), "boost::math::ellint_3<%1%>(%1%,%1%,%1%)");
0358 }
0359
0360 template <class T1, class T2, class T3>
0361 typename detail::ellint_3_result<T1, T2, T3>::type ellint_3(T1 k, T2 v, T3 phi)
0362 {
0363    typedef typename policies::is_policy<T3>::type tag_type;
0364    return detail::ellint_3(k, v, phi, tag_type());
0365 }
0366
0367 template <class T1, class T2>
0368 inline typename tools::promote_args<T1, T2>::type ellint_3(T1 k, T2 v)
0369 {
0370    return ellint_3(k, v, policies::policy<>());
0371 }
0372
0373 }} // namespaces
0374
0375 #endif // BOOST_MATH_ELLINT_3_HPP
0376