boost/math/special_functions/ellint_1.hpp

   1 //  Copyright (c) 2006 Xiaogang Zhang
   2 //  Copyright (c) 2006 John Maddock
   3 //  Use, modification and distribution are subject to the
   4 //  Boost Software License, Version 1.0. (See accompanying file
   5 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6 //
   7 //  History:
   8 //  XZ wrote the original of this file as part of the Google
   9 //  Summer of Code 2006.  JM modified it to fit into the
  10 //  Boost.Math conceptual framework better, and to ensure
  11 //  that the code continues to work no matter how many digits
  12 //  type T has.
  13
  14 #ifndef BOOST_MATH_ELLINT_1_HPP
  15 #define BOOST_MATH_ELLINT_1_HPP
  16
  17 #ifdef _MSC_VER
  18 #pragma once
  19 #endif
  20
  21 #include <boost/math/special_functions/ellint_rf.hpp>
  22 #include <boost/math/constants/constants.hpp>
  23 #include <boost/math/policies/error_handling.hpp>
  24 #include <boost/math/tools/workaround.hpp>
  25 #include <boost/math/special_functions/round.hpp>
  26
  27 // Elliptic integrals (complete and incomplete) of the first kind
  28 // Carlson, Numerische Mathematik, vol 33, 1 (1979)
  29
  30 namespace boost { namespace math {
  31
  32 template <class T1, class T2, class Policy>
  33 typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi, const Policy& pol);
  34
  35 namespace detail{
  36
  37 template <typename T, typename Policy>
  38 T ellint_k_imp(T k, const Policy& pol);
  39
  40 // Elliptic integral (Legendre form) of the first kind
  41 template <typename T, typename Policy>
  42 T ellint_f_imp(T phi, T k, const Policy& pol)
  43 {
  44     BOOST_MATH_STD_USING
  45     using namespace boost::math::tools;
  46     using namespace boost::math::constants;
  47
  48     static const char* function = "boost::math::ellint_f<%1%>(%1%,%1%)";
  49     BOOST_MATH_INSTRUMENT_VARIABLE(phi);
  50     BOOST_MATH_INSTRUMENT_VARIABLE(k);
  51     BOOST_MATH_INSTRUMENT_VARIABLE(function);
  52
  53     if (abs(k) > 1)
  54     {
  55        return policies::raise_domain_error<T>(function,
  56             "Got k = %1%, function requires |k| <= 1", k, pol);
  57     }
  58
  59     bool invert = false;
  60     if(phi < 0)
  61     {
  62        BOOST_MATH_INSTRUMENT_VARIABLE(phi);
  63        phi = fabs(phi);
  64        invert = true;
  65     }
  66
  67     T result;
  68
  69     if(phi >= tools::max_value<T>())
  70     {
  71        // Need to handle infinity as a special case:
  72        result = policies::raise_overflow_error<T>(function, 0, pol);
  73        BOOST_MATH_INSTRUMENT_VARIABLE(result);
  74     }
  75     else if(phi > 1 / tools::epsilon<T>())
  76     {
  77        // Phi is so large that phi%pi is necessarily zero (or garbage),
  78        // just return the second part of the duplication formula:
  79        result = 2 * phi * ellint_k_imp(k, pol) / constants::pi<T>();
  80        BOOST_MATH_INSTRUMENT_VARIABLE(result);
  81     }
  82     else
  83     {
  84        // Carlson's algorithm works only for |phi| <= pi/2,
  85        // use the integrand's periodicity to normalize phi
  86        //
  87        // Xiaogang's original code used a cast to long long here
  88        // but that fails if T has more digits than a long long,
  89        // so rewritten to use fmod instead:
  90        //
  91        BOOST_MATH_INSTRUMENT_CODE("pi/2 = " << constants::pi<T>() / 2);
  92        T rphi = boost::math::tools::fmod_workaround(phi, T(constants::half_pi<T>()));
  93        BOOST_MATH_INSTRUMENT_VARIABLE(rphi);
  94        T m = boost::math::round((phi - rphi) / constants::half_pi<T>());
  95        BOOST_MATH_INSTRUMENT_VARIABLE(m);
  96        int s = 1;
  97        if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5)
  98        {
  99           m += 1;
 100           s = -1;
 101           rphi = constants::half_pi<T>() - rphi;
 102           BOOST_MATH_INSTRUMENT_VARIABLE(rphi);
 103        }
 104        T sinp = sin(rphi);
 105        T cosp = cos(rphi);
 106        BOOST_MATH_INSTRUMENT_VARIABLE(sinp);
 107        BOOST_MATH_INSTRUMENT_VARIABLE(cosp);
 108        result = s * sinp * ellint_rf_imp(T(cosp * cosp), T(1 - k * k * sinp * sinp), T(1), pol);
 109        BOOST_MATH_INSTRUMENT_VARIABLE(result);
 110        if(m != 0)
 111        {
 112           result += m * ellint_k_imp(k, pol);
 113           BOOST_MATH_INSTRUMENT_VARIABLE(result);
 114        }
 115     }
 116     return invert ? T(-result) : result;
 117 }
 118
 119 // Complete elliptic integral (Legendre form) of the first kind
 120 template <typename T, typename Policy>
 121 T ellint_k_imp(T k, const Policy& pol)
 122 {
 123     BOOST_MATH_STD_USING
 124     using namespace boost::math::tools;
 125
 126     static const char* function = "boost::math::ellint_k<%1%>(%1%)";
 127
 128     if (abs(k) > 1)
 129     {
 130        return policies::raise_domain_error<T>(function,
 131             "Got k = %1%, function requires |k| <= 1", k, pol);
 132     }
 133     if (abs(k) == 1)
 134     {
 135        return policies::raise_overflow_error<T>(function, 0, pol);
 136     }
 137
 138     T x = 0;
 139     T y = 1 - k * k;
 140     T z = 1;
 141     T value = ellint_rf_imp(x, y, z, pol);
 142
 143     return value;
 144 }
 145
 146 template <typename T, typename Policy>
 147 inline typename tools::promote_args<T>::type ellint_1(T k, const Policy& pol, const mpl::true_&)
 148 {
 149    typedef typename tools::promote_args<T>::type result_type;
 150    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 151    return policies::checked_narrowing_cast<result_type, Policy>(detail::ellint_k_imp(static_cast<value_type>(k), pol), "boost::math::ellint_1<%1%>(%1%)");
 152 }
 153
 154 template <class T1, class T2>
 155 inline typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi, const mpl::false_&)
 156 {
 157    return boost::math::ellint_1(k, phi, policies::policy<>());
 158 }
 159
 160 }
 161
 162 // Complete elliptic integral (Legendre form) of the first kind
 163 template <typename T>
 164 inline typename tools::promote_args<T>::type ellint_1(T k)
 165 {
 166    return ellint_1(k, policies::policy<>());
 167 }
 168
 169 // Elliptic integral (Legendre form) of the first kind
 170 template <class T1, class T2, class Policy>
 171 inline typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi, const Policy& pol)
 172 {
 173    typedef typename tools::promote_args<T1, T2>::type result_type;
 174    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 175    return policies::checked_narrowing_cast<result_type, Policy>(detail::ellint_f_imp(static_cast<value_type>(phi), static_cast<value_type>(k), pol), "boost::math::ellint_1<%1%>(%1%,%1%)");
 176 }
 177
 178 template <class T1, class T2>
 179 inline typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi)
 180 {
 181    typedef typename policies::is_policy<T2>::type tag_type;
 182    return detail::ellint_1(k, phi, tag_type());
 183 }
 184
 185 }} // namespaces
 186
 187 #endif // BOOST_MATH_ELLINT_1_HPP
 188