math/special_functions/ellint_2.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 ensure
0011 //  that the code continues to work no matter how many digits
0012 //  type T has.
0013
0014 #ifndef BOOST_MATH_ELLINT_2_HPP
0015 #define BOOST_MATH_ELLINT_2_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_rd.hpp>
0024 #include <boost/math/special_functions/ellint_rg.hpp>
0025 #include <boost/math/constants/constants.hpp>
0026 #include <boost/math/policies/error_handling.hpp>
0027 #include <boost/math/tools/workaround.hpp>
0028 #include <boost/math/special_functions/round.hpp>
0029
0030 // Elliptic integrals (complete and incomplete) of the second kind
0031 // Carlson, Numerische Mathematik, vol 33, 1 (1979)
0032
0033 namespace boost { namespace math {
0034
0035 template <class T1, class T2, class Policy>
0036 typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi, const Policy& pol);
0037
0038 namespace detail{
0039
0040 template <typename T, typename Policy>
0041 T ellint_e_imp(T k, const Policy& pol, const std::integral_constant<int, 0>&);
0042 template <typename T, typename Policy>
0043 T ellint_e_imp(T k, const Policy& pol, const std::integral_constant<int, 1>&);
0044 template <typename T, typename Policy>
0045 T ellint_e_imp(T k, const Policy& pol, const std::integral_constant<int, 2>&);
0046
0047 // Elliptic integral (Legendre form) of the second kind
0048 template <typename T, typename Policy>
0049 T ellint_e_imp(T phi, T k, const Policy& pol)
0050 {
0051     BOOST_MATH_STD_USING
0052     using namespace boost::math::tools;
0053     using namespace boost::math::constants;
0054
0055     bool invert = false;
0056     if (phi == 0)
0057        return 0;
0058
0059     if(phi < 0)
0060     {
0061        phi = fabs(phi);
0062        invert = true;
0063     }
0064
0065     T result;
0066
0067     if(phi >= tools::max_value<T>())
0068     {
0069        // Need to handle infinity as a special case:
0070        result = policies::raise_overflow_error<T>("boost::math::ellint_e<%1%>(%1%,%1%)", nullptr, pol);
0071     }
0072     else if(phi > 1 / tools::epsilon<T>())
0073     {
0074        typedef std::integral_constant<int,
0075           std::is_floating_point<T>::value&& std::numeric_limits<T>::digits && (std::numeric_limits<T>::digits <= 54) ? 0 :
0076           std::is_floating_point<T>::value && std::numeric_limits<T>::digits && (std::numeric_limits<T>::digits <= 64) ? 1 : 2
0077        > precision_tag_type;
0078        // Phi is so large that phi%pi is necessarily zero (or garbage),
0079        // just return the second part of the duplication formula:
0080        result = 2 * phi * ellint_e_imp(k, pol, precision_tag_type()) / constants::pi<T>();
0081     }
0082     else if(k == 0)
0083     {
0084        return invert ? T(-phi) : phi;
0085     }
0086     else if(fabs(k) == 1)
0087     {
0088        //
0089        // For k = 1 ellipse actually turns to a line and every pi/2 in phi is exactly 1 in arc length
0090        // Periodicity though is in pi, curve follows sin(pi) for 0 <= phi <= pi/2 and then
0091        // 2 - sin(pi- phi) = 2 + sin(phi - pi) for pi/2 <= phi <= pi, so general form is:
0092        //
0093        // 2n + sin(phi - n * pi) ; |phi - n * pi| <= pi / 2
0094        //
0095        T m = boost::math::round(phi / boost::math::constants::pi<T>());
0096        T remains = phi - m * boost::math::constants::pi<T>();
0097        T value = 2 * m + sin(remains);
0098
0099        // negative arc length for negative phi
0100        return invert ? -value : value;
0101     }
0102     else
0103     {
0104        // Carlson's algorithm works only for |phi| <= pi/2,
0105        // use the integrand's periodicity to normalize phi
0106        //
0107        // Xiaogang's original code used a cast to long long here
0108        // but that fails if T has more digits than a long long,
0109        // so rewritten to use fmod instead:
0110        //
0111        T rphi = boost::math::tools::fmod_workaround(phi, T(constants::half_pi<T>()));
0112        T m = boost::math::round((phi - rphi) / constants::half_pi<T>());
0113        int s = 1;
0114        if(boost::math::tools::fmod_workaround(m, T(2)) > T(0.5))
0115        {
0116           m += 1;
0117           s = -1;
0118           rphi = constants::half_pi<T>() - rphi;
0119        }
0120        T k2 = k * k;
0121        if(boost::math::pow<3>(rphi) * k2 / 6 < tools::epsilon<T>() * fabs(rphi))
0122        {
0123           // See http://functions.wolfram.com/EllipticIntegrals/EllipticE2/06/01/03/0001/
0124           result = s * rphi;
0125        }
0126        else
0127        {
0128           // http://dlmf.nist.gov/19.25#E10
0129           T sinp = sin(rphi);
0130           if (k2 * sinp * sinp >= 1)
0131           {
0132              return policies::raise_domain_error<T>("boost::math::ellint_2<%1%>(%1%, %1%)", "The parameter k is out of range, got k = %1%", k, pol);
0133           }
0134           T cosp = cos(rphi);
0135           T c = 1 / (sinp * sinp);
0136           T cm1 = cosp * cosp / (sinp * sinp);  // c - 1
0137           result = s * ((1 - k2) * ellint_rf_imp(cm1, T(c - k2), c, pol) + k2 * (1 - k2) * ellint_rd(cm1, c, T(c - k2), pol) / 3 + k2 * sqrt(cm1 / (c * (c - k2))));
0138        }
0139        if (m != 0)
0140        {
0141           typedef std::integral_constant<int,
0142              std::is_floating_point<T>::value&& std::numeric_limits<T>::digits && (std::numeric_limits<T>::digits <= 54) ? 0 :
0143              std::is_floating_point<T>::value && std::numeric_limits<T>::digits && (std::numeric_limits<T>::digits <= 64) ? 1 : 2
0144           > precision_tag_type;
0145           result += m * ellint_e_imp(k, pol, precision_tag_type());
0146        }
0147     }
0148     return invert ? T(-result) : result;
0149 }
0150
0151 // Complete elliptic integral (Legendre form) of the second kind
0152 template <typename T, typename Policy>
0153 T ellint_e_imp(T k, const Policy& pol, std::integral_constant<int, 2> const&)
0154 {
0155     BOOST_MATH_STD_USING
0156     using namespace boost::math::tools;
0157
0158     if (abs(k) > 1)
0159     {
0160        return policies::raise_domain_error<T>("boost::math::ellint_e<%1%>(%1%)",
0161             "Got k = %1%, function requires |k| <= 1", k, pol);
0162     }
0163     if (abs(k) == 1)
0164     {
0165         return static_cast<T>(1);
0166     }
0167
0168     T x = 0;
0169     T t = k * k;
0170     T y = 1 - t;
0171     T z = 1;
0172     T value = 2 * ellint_rg_imp(x, y, z, pol);
0173
0174     return value;
0175 }
0176 //
0177 // Special versions for double and 80-bit long double precision,
0178 // double precision versions use the coefficients from:
0179 // "Fast computation of complete elliptic integrals and Jacobian elliptic functions",
0180 // Celestial Mechanics and Dynamical Astronomy, April 2012.
0181 //
0182 // Higher precision coefficients for 80-bit long doubles can be calculated
0183 // using for example:
0184 // Table[N[SeriesCoefficient[ EllipticE [ m ], { m, 875/1000, i} ], 20], {i, 0, 24}]
0185 // and checking the value of the first neglected term with:
0186 // N[SeriesCoefficient[ EllipticE [ m ], { m, 875/1000, 24} ], 20] * (2.5/100)^24
0187 //
0188 // For m > 0.9 we don't use the method of the paper above, but simply call our
0189 // existing routines.
0190 //
0191 template <typename T, typename Policy>
0192 BOOST_FORCEINLINE T ellint_e_imp(T k, const Policy& pol, std::integral_constant<int, 0> const&)
0193 {
0194    using std::abs;
0195    using namespace boost::math::tools;
0196
0197    T m = k * k;
0198    switch (static_cast<int>(20 * m))
0199    {
0200    case 0:
0201    case 1:
0202    //if (m < 0.1)
0203    {
0204       constexpr T coef[] =
0205       {
0206          static_cast<T>(1.550973351780472328),
0207          -static_cast<T>(0.400301020103198524),
0208          -static_cast<T>(0.078498619442941939),
0209          -static_cast<T>(0.034318853117591992),
0210          -static_cast<T>(0.019718043317365499),
0211          -static_cast<T>(0.013059507731993309),
0212          -static_cast<T>(0.009442372874146547),
0213          -static_cast<T>(0.007246728512402157),
0214          -static_cast<T>(0.005807424012956090),
0215          -static_cast<T>(0.004809187786009338),
0216          -static_cast<T>(0.004086399233255150)
0217       };
0218       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.05));
0219    }
0220    case 2:
0221    case 3:
0222    //else if (m < 0.2)
0223    {
0224       constexpr T coef[] =
0225       {
0226          static_cast<T>(1.510121832092819728),
0227          -static_cast<T>(0.417116333905867549),
0228          -static_cast<T>(0.090123820404774569),
0229          -static_cast<T>(0.043729944019084312),
0230          -static_cast<T>(0.027965493064761785),
0231          -static_cast<T>(0.020644781177568105),
0232          -static_cast<T>(0.016650786739707238),
0233          -static_cast<T>(0.014261960828842520),
0234          -static_cast<T>(0.012759847429264803),
0235          -static_cast<T>(0.011799303775587354),
0236          -static_cast<T>(0.011197445703074968)
0237       };
0238       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.15));
0239    }
0240    case 4:
0241    case 5:
0242    //else if (m < 0.3)
0243    {
0244       constexpr T coef[] =
0245       {
0246          static_cast<T>(1.467462209339427155),
0247          -static_cast<T>(0.436576290946337775),
0248          -static_cast<T>(0.105155557666942554),
0249          -static_cast<T>(0.057371843593241730),
0250          -static_cast<T>(0.041391627727340220),
0251          -static_cast<T>(0.034527728505280841),
0252          -static_cast<T>(0.031495443512532783),
0253          -static_cast<T>(0.030527000890325277),
0254          -static_cast<T>(0.030916984019238900),
0255          -static_cast<T>(0.032371395314758122),
0256          -static_cast<T>(0.034789960386404158)
0257       };
0258       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.25));
0259    }
0260    case 6:
0261    case 7:
0262    //else if (m < 0.4)
0263    {
0264       constexpr T coef[] =
0265       {
0266          static_cast<T>(1.422691133490879171),
0267          -static_cast<T>(0.459513519621048674),
0268          -static_cast<T>(0.125250539822061878),
0269          -static_cast<T>(0.078138545094409477),
0270          -static_cast<T>(0.064714278472050002),
0271          -static_cast<T>(0.062084339131730311),
0272          -static_cast<T>(0.065197032815572477),
0273          -static_cast<T>(0.072793895362578779),
0274          -static_cast<T>(0.084959075171781003),
0275          -static_cast<T>(0.102539850131045997),
0276          -static_cast<T>(0.127053585157696036),
0277          -static_cast<T>(0.160791120691274606)
0278       };
0279       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.35));
0280    }
0281    case 8:
0282    case 9:
0283    //else if (m < 0.5)
0284    {
0285       constexpr T coef[] =
0286       {
0287          static_cast<T>(1.375401971871116291),
0288          -static_cast<T>(0.487202183273184837),
0289          -static_cast<T>(0.153311701348540228),
0290          -static_cast<T>(0.111849444917027833),
0291          -static_cast<T>(0.108840952523135768),
0292          -static_cast<T>(0.122954223120269076),
0293          -static_cast<T>(0.152217163962035047),
0294          -static_cast<T>(0.200495323642697339),
0295          -static_cast<T>(0.276174333067751758),
0296          -static_cast<T>(0.393513114304375851),
0297          -static_cast<T>(0.575754406027879147),
0298          -static_cast<T>(0.860523235727239756),
0299          -static_cast<T>(1.308833205758540162)
0300       };
0301       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.45));
0302    }
0303    case 10:
0304    case 11:
0305    //else if (m < 0.6)
0306    {
0307       constexpr T coef[] =
0308       {
0309          static_cast<T>(1.325024497958230082),
0310          -static_cast<T>(0.521727647557566767),
0311          -static_cast<T>(0.194906430482126213),
0312          -static_cast<T>(0.171623726822011264),
0313          -static_cast<T>(0.202754652926419141),
0314          -static_cast<T>(0.278798953118534762),
0315          -static_cast<T>(0.420698457281005762),
0316          -static_cast<T>(0.675948400853106021),
0317          -static_cast<T>(1.136343121839229244),
0318          -static_cast<T>(1.976721143954398261),
0319          -static_cast<T>(3.531696773095722506),
0320          -static_cast<T>(6.446753640156048150),
0321          -static_cast<T>(11.97703130208884026)
0322       };
0323       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.55));
0324    }
0325    case 12:
0326    case 13:
0327    //else if (m < 0.7)
0328    {
0329       constexpr T coef[] =
0330       {
0331          static_cast<T>(1.270707479650149744),
0332          -static_cast<T>(0.566839168287866583),
0333          -static_cast<T>(0.262160793432492598),
0334          -static_cast<T>(0.292244173533077419),
0335          -static_cast<T>(0.440397840850423189),
0336          -static_cast<T>(0.774947641381397458),
0337          -static_cast<T>(1.498870837987561088),
0338          -static_cast<T>(3.089708310445186667),
0339          -static_cast<T>(6.667595903381001064),
0340          -static_cast<T>(14.89436036517319078),
0341          -static_cast<T>(34.18120574251449024),
0342          -static_cast<T>(80.15895841905397306),
0343          -static_cast<T>(191.3489480762984920),
0344          -static_cast<T>(463.5938853480342030),
0345          -static_cast<T>(1137.380822169360061)
0346       };
0347       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.65));
0348    }
0349    case 14:
0350    case 15:
0351    //else if (m < 0.8)
0352    {
0353       constexpr T coef[] =
0354       {
0355          static_cast<T>(1.211056027568459525),
0356          -static_cast<T>(0.630306413287455807),
0357          -static_cast<T>(0.387166409520669145),
0358          -static_cast<T>(0.592278235311934603),
0359          -static_cast<T>(1.237555584513049844),
0360          -static_cast<T>(3.032056661745247199),
0361          -static_cast<T>(8.181688221573590762),
0362          -static_cast<T>(23.55507217389693250),
0363          -static_cast<T>(71.04099935893064956),
0364          -static_cast<T>(221.8796853192349888),
0365          -static_cast<T>(712.1364793277635425),
0366          -static_cast<T>(2336.125331440396407),
0367          -static_cast<T>(7801.945954775964673),
0368          -static_cast<T>(26448.19586059191933),
0369          -static_cast<T>(90799.48341621365251),
0370          -static_cast<T>(315126.0406449163424),
0371          -static_cast<T>(1104011.344311591159)
0372       };
0373       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.75));
0374    }
0375    case 16:
0376    //else if (m < 0.85)
0377    {
0378       constexpr T coef[] =
0379       {
0380          static_cast<T>(1.161307152196282836),
0381          -static_cast<T>(0.701100284555289548),
0382          -static_cast<T>(0.580551474465437362),
0383          -static_cast<T>(1.243693061077786614),
0384          -static_cast<T>(3.679383613496634879),
0385          -static_cast<T>(12.81590924337895775),
0386          -static_cast<T>(49.25672530759985272),
0387          -static_cast<T>(202.1818735434090269),
0388          -static_cast<T>(869.8602699308701437),
0389          -static_cast<T>(3877.005847313289571),
0390          -static_cast<T>(17761.70710170939814),
0391          -static_cast<T>(83182.69029154232061),
0392          -static_cast<T>(396650.4505013548170),
0393          -static_cast<T>(1920033.413682634405)
0394       };
0395       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.825));
0396    }
0397    case 17:
0398    //else if (m < 0.90)
0399    {
0400       constexpr T coef[] =
0401       {
0402          static_cast<T>(1.124617325119752213),
0403          -static_cast<T>(0.770845056360909542),
0404          -static_cast<T>(0.844794053644911362),
0405          -static_cast<T>(2.490097309450394453),
0406          -static_cast<T>(10.23971741154384360),
0407          -static_cast<T>(49.74900546551479866),
0408          -static_cast<T>(267.0986675195705196),
0409          -static_cast<T>(1532.665883825229947),
0410          -static_cast<T>(9222.313478526091951),
0411          -static_cast<T>(57502.51612140314030),
0412          -static_cast<T>(368596.1167416106063),
0413          -static_cast<T>(2415611.088701091428),
0414          -static_cast<T>(16120097.81581656797),
0415          -static_cast<T>(109209938.5203089915),
0416          -static_cast<T>(749380758.1942496220),
0417          -static_cast<T>(5198725846.725541393),
0418          -static_cast<T>(36409256888.12139973)
0419       };
0420       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.875));
0421    }
0422    default:
0423       //
0424       // All cases where m > 0.9
0425       // including all error handling:
0426       //
0427       return ellint_e_imp(k, pol, std::integral_constant<int, 2>());
0428    }
0429 }
0430 template <typename T, typename Policy>
0431 BOOST_FORCEINLINE T ellint_e_imp(T k, const Policy& pol, std::integral_constant<int, 1> const&)
0432 {
0433    using std::abs;
0434    using namespace boost::math::tools;
0435
0436    T m = k * k;
0437    switch (static_cast<int>(20 * m))
0438    {
0439    case 0:
0440    case 1:
0441       //if (m < 0.1)
0442    {
0443       constexpr T coef[] =
0444       {
0445          1.5509733517804723277L,
0446          -0.40030102010319852390L,
0447          -0.078498619442941939212L,
0448          -0.034318853117591992417L,
0449          -0.019718043317365499309L,
0450          -0.013059507731993309191L,
0451          -0.0094423728741465473894L,
0452          -0.0072467285124021568126L,
0453          -0.0058074240129560897940L,
0454          -0.0048091877860093381762L,
0455          -0.0040863992332551506768L,
0456          -0.0035450302604139562644L,
0457          -0.0031283511188028336315L
0458       };
0459       return boost::math::tools::evaluate_polynomial(coef, m - 0.05L);
0460    }
0461    case 2:
0462    case 3:
0463       //else if (m < 0.2)
0464    {
0465       constexpr T coef[] =
0466       {
0467          1.5101218320928197276L,
0468          -0.41711633390586754922L,
0469          -0.090123820404774568894L,
0470          -0.043729944019084311555L,
0471          -0.027965493064761784548L,
0472          -0.020644781177568105268L,
0473          -0.016650786739707238037L,
0474          -0.014261960828842519634L,
0475          -0.012759847429264802627L,
0476          -0.011799303775587354169L,
0477          -0.011197445703074968018L,
0478          -0.010850368064799902735L,
0479          -0.010696133481060989818L
0480       };
0481       return boost::math::tools::evaluate_polynomial(coef, m - 0.15L);
0482    }
0483    case 4:
0484    case 5:
0485       //else if (m < 0.3L)
0486    {
0487       constexpr T coef[] =
0488       {
0489          1.4674622093394271555L,
0490          -0.43657629094633777482L,
0491          -0.10515555766694255399L,
0492          -0.057371843593241729895L,
0493          -0.041391627727340220236L,
0494          -0.034527728505280841188L,
0495          -0.031495443512532782647L,
0496          -0.030527000890325277179L,
0497          -0.030916984019238900349L,
0498          -0.032371395314758122268L,
0499          -0.034789960386404158240L,
0500          -0.038182654612387881967L,
0501          -0.042636187648900252525L,
0502          -0.048302272505241634467
0503       };
0504       return boost::math::tools::evaluate_polynomial(coef, m - 0.25L);
0505    }
0506    case 6:
0507    case 7:
0508       //else if (m < 0.4L)
0509    {
0510       constexpr T coef[] =
0511       {
0512          1.4226911334908791711L,
0513          -0.45951351962104867394L,
0514          -0.12525053982206187849L,
0515          -0.078138545094409477156L,
0516          -0.064714278472050001838L,
0517          -0.062084339131730310707L,
0518          -0.065197032815572476910L,
0519          -0.072793895362578779473L,
0520          -0.084959075171781003264L,
0521          -0.10253985013104599679L,
0522          -0.12705358515769603644L,
0523          -0.16079112069127460621L,
0524          -0.20705400012405941376L,
0525          -0.27053164884730888948L
0526       };
0527       return boost::math::tools::evaluate_polynomial(coef, m - 0.35L);
0528    }
0529    case 8:
0530    case 9:
0531       //else if (m < 0.5L)
0532    {
0533       constexpr T coef[] =
0534       {
0535          1.3754019718711162908L,
0536          -0.48720218327318483652L,
0537          -0.15331170134854022753L,
0538          -0.11184944491702783273L,
0539          -0.10884095252313576755L,
0540          -0.12295422312026907610L,
0541          -0.15221716396203504746L,
0542          -0.20049532364269733857L,
0543          -0.27617433306775175837L,
0544          -0.39351311430437585139L,
0545          -0.57575440602787914711L,
0546          -0.86052323572723975634L,
0547          -1.3088332057585401616L,
0548          -2.0200280559452241745L,
0549          -3.1566019548237606451L
0550       };
0551       return boost::math::tools::evaluate_polynomial(coef, m - 0.45L);
0552    }
0553    case 10:
0554    case 11:
0555       //else if (m < 0.6L)
0556    {
0557       constexpr T coef[] =
0558       {
0559          1.3250244979582300818L,
0560          -0.52172764755756676713L,
0561          -0.19490643048212621262L,
0562          -0.17162372682201126365L,
0563          -0.20275465292641914128L,
0564          -0.27879895311853476205L,
0565          -0.42069845728100576224L,
0566          -0.67594840085310602110L,
0567          -1.1363431218392292440L,
0568          -1.9767211439543982613L,
0569          -3.5316967730957225064L,
0570          -6.4467536401560481499L,
0571          -11.977031302088840261L,
0572          -22.581360948073964469L,
0573          -43.109479829481450573L,
0574          -83.186290908288807424L
0575       };
0576       return boost::math::tools::evaluate_polynomial(coef, m - 0.55L);
0577    }
0578    case 12:
0579    case 13:
0580       //else if (m < 0.7L)
0581    {
0582       constexpr T coef[] =
0583       {
0584          1.2707074796501497440L,
0585          -0.56683916828786658286L,
0586          -0.26216079343249259779L,
0587          -0.29224417353307741931L,
0588          -0.44039784085042318909L,
0589          -0.77494764138139745824L,
0590          -1.4988708379875610880L,
0591          -3.0897083104451866665L,
0592          -6.6675959033810010645L,
0593          -14.894360365173190775L,
0594          -34.181205742514490240L,
0595          -80.158958419053973056L,
0596          -191.34894807629849204L,
0597          -463.59388534803420301L,
0598          -1137.3808221693600606L,
0599          -2820.7073786352269339L,
0600          -7061.1382244658715621L,
0601          -17821.809331816437058L,
0602          -45307.849987201897801L
0603       };
0604       return boost::math::tools::evaluate_polynomial(coef, m - 0.65L);
0605    }
0606    case 14:
0607    case 15:
0608       //else if (m < 0.8L)
0609    {
0610       constexpr T coef[] =
0611       {
0612          1.2110560275684595248L,
0613          -0.63030641328745580709L,
0614          -0.38716640952066914514L,
0615          -0.59227823531193460257L,
0616          -1.2375555845130498445L,
0617          -3.0320566617452471986L,
0618          -8.1816882215735907624L,
0619          -23.555072173896932503L,
0620          -71.040999358930649565L,
0621          -221.87968531923498875L,
0622          -712.13647932776354253L,
0623          -2336.1253314403964072L,
0624          -7801.9459547759646726L,
0625          -26448.195860591919335L,
0626          -90799.483416213652512L,
0627          -315126.04064491634241L,
0628          -1.1040113443115911589e6L,
0629          -3.8998018348056769095e6L,
0630          -1.3876249116223745041e7L,
0631          -4.9694982823537861149e7L,
0632          -1.7900668836197342979e8L,
0633          -6.4817399873722371964e8L
0634       };
0635       return boost::math::tools::evaluate_polynomial(coef, m - 0.75L);
0636    }
0637    case 16:
0638       //else if (m < 0.85L)
0639    {
0640       constexpr T coef[] =
0641       {
0642          1.1613071521962828360L,
0643          -0.70110028455528954752L,
0644          -0.58055147446543736163L,
0645          -1.2436930610777866138L,
0646          -3.6793836134966348789L,
0647          -12.815909243378957753L,
0648          -49.256725307599852720L,
0649          -202.18187354340902693L,
0650          -869.86026993087014372L,
0651          -3877.0058473132895713L,
0652          -17761.707101709398174L,
0653          -83182.690291542320614L,
0654          -396650.45050135481698L,
0655          -1.9200334136826344054e6L,
0656          -9.4131321779500838352e6L,
0657          -4.6654858837335370627e7L,
0658          -2.3343549352617609390e8L,
0659          -1.1776928223045913454e9L,
0660          -5.9850851892915740401e9L,
0661          -3.0614702984618644983e10L
0662       };
0663       return boost::math::tools::evaluate_polynomial(coef, m - 0.825L);
0664    }
0665    case 17:
0666       //else if (m < 0.90L)
0667    {
0668       constexpr T coef[] =
0669       {
0670          1.1246173251197522132L,
0671          -0.77084505636090954218L,
0672          -0.84479405364491136236L,
0673          -2.4900973094503944527L,
0674          -10.239717411543843601L,
0675          -49.749005465514798660L,
0676          -267.09866751957051961L,
0677          -1532.6658838252299468L,
0678          -9222.3134785260919507L,
0679          -57502.516121403140303L,
0680          -368596.11674161060626L,
0681          -2.4156110887010914281e6L,
0682          -1.6120097815816567971e7L,
0683          -1.0920993852030899148e8L,
0684          -7.4938075819424962198e8L,
0685          -5.1987258467255413931e9L,
0686          -3.6409256888121399726e10L,
0687          -2.5711802891217393544e11L,
0688          -1.8290904062978796996e12L,
0689          -1.3096838781743248404e13L,
0690          -9.4325465851415135118e13L,
0691          -6.8291980829471896669e14L
0692       };
0693       return boost::math::tools::evaluate_polynomial(coef, m - 0.875L);
0694    }
0695    default:
0696       //
0697       // All cases where m > 0.9
0698       // including all error handling:
0699       //
0700       return ellint_e_imp(k, pol, std::integral_constant<int, 2>());
0701    }
0702 }
0703
0704 template <typename T, typename Policy>
0705 BOOST_FORCEINLINE typename tools::promote_args<T>::type ellint_2(T k, const Policy& pol, const std::true_type&)
0706 {
0707    typedef typename tools::promote_args<T>::type result_type;
0708    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0709    typedef std::integral_constant<int,
0710       std::is_floating_point<T>::value&& std::numeric_limits<T>::digits && (std::numeric_limits<T>::digits <= 54) ? 0 :
0711       std::is_floating_point<T>::value && std::numeric_limits<T>::digits && (std::numeric_limits<T>::digits <= 64) ? 1 : 2
0712    > precision_tag_type;
0713    return policies::checked_narrowing_cast<result_type, Policy>(detail::ellint_e_imp(static_cast<value_type>(k), pol, precision_tag_type()), "boost::math::ellint_2<%1%>(%1%)");
0714 }
0715
0716 // Elliptic integral (Legendre form) of the second kind
0717 template <class T1, class T2>
0718 BOOST_FORCEINLINE typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi, const std::false_type&)
0719 {
0720    return boost::math::ellint_2(k, phi, policies::policy<>());
0721 }
0722
0723 } // detail
0724
0725 // Complete elliptic integral (Legendre form) of the second kind
0726 template <typename T>
0727 BOOST_FORCEINLINE typename tools::promote_args<T>::type ellint_2(T k)
0728 {
0729    return ellint_2(k, policies::policy<>());
0730 }
0731
0732 // Elliptic integral (Legendre form) of the second kind
0733 template <class T1, class T2>
0734 BOOST_FORCEINLINE typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi)
0735 {
0736    typedef typename policies::is_policy<T2>::type tag_type;
0737    return detail::ellint_2(k, phi, tag_type());
0738 }
0739
0740 template <class T1, class T2, class Policy>
0741 BOOST_FORCEINLINE typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi, const Policy& pol)
0742 {
0743    typedef typename tools::promote_args<T1, T2>::type result_type;
0744    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0745    return policies::checked_narrowing_cast<result_type, Policy>(detail::ellint_e_imp(static_cast<value_type>(phi), static_cast<value_type>(k), pol), "boost::math::ellint_2<%1%>(%1%,%1%)");
0746 }
0747
0748 }} // namespaces
0749
0750 #endif // BOOST_MATH_ELLINT_2_HPP
0751