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%)", "Got k = %1%, function requires |k| <= 1", k, pol);
0161     }
0162     if (abs(k) == 1)
0163     {
0164         return static_cast<T>(1);
0165     }
0166
0167     T x = 0;
0168     T t = k * k;
0169     T y = 1 - t;
0170     T z = 1;
0171     T value = 2 * ellint_rg_imp(x, y, z, pol);
0172
0173     return value;
0174 }
0175 //
0176 // Special versions for double and 80-bit long double precision,
0177 // double precision versions use the coefficients from:
0178 // "Fast computation of complete elliptic integrals and Jacobian elliptic functions",
0179 // Celestial Mechanics and Dynamical Astronomy, April 2012.
0180 //
0181 // Higher precision coefficients for 80-bit long doubles can be calculated
0182 // using for example:
0183 // Table[N[SeriesCoefficient[ EllipticE [ m ], { m, 875/1000, i} ], 20], {i, 0, 24}]
0184 // and checking the value of the first neglected term with:
0185 // N[SeriesCoefficient[ EllipticE [ m ], { m, 875/1000, 24} ], 20] * (2.5/100)^24
0186 //
0187 // For m > 0.9 we don't use the method of the paper above, but simply call our
0188 // existing routines.
0189 //
0190 template <typename T, typename Policy>
0191 BOOST_MATH_FORCEINLINE T ellint_e_imp(T k, const Policy& pol, std::integral_constant<int, 0> const&)
0192 {
0193    using std::abs;
0194    using namespace boost::math::tools;
0195
0196    T m = k * k;
0197    switch (static_cast<int>(20 * m))
0198    {
0199    case 0:
0200    case 1:
0201    //if (m < 0.1)
0202    {
0203       constexpr T coef[] =
0204       {
0205          static_cast<T>(1.550973351780472328),
0206          -static_cast<T>(0.400301020103198524),
0207          -static_cast<T>(0.078498619442941939),
0208          -static_cast<T>(0.034318853117591992),
0209          -static_cast<T>(0.019718043317365499),
0210          -static_cast<T>(0.013059507731993309),
0211          -static_cast<T>(0.009442372874146547),
0212          -static_cast<T>(0.007246728512402157),
0213          -static_cast<T>(0.005807424012956090),
0214          -static_cast<T>(0.004809187786009338),
0215          -static_cast<T>(0.004086399233255150)
0216       };
0217       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.05));
0218    }
0219    case 2:
0220    case 3:
0221    //else if (m < 0.2)
0222    {
0223       constexpr T coef[] =
0224       {
0225          static_cast<T>(1.510121832092819728),
0226          -static_cast<T>(0.417116333905867549),
0227          -static_cast<T>(0.090123820404774569),
0228          -static_cast<T>(0.043729944019084312),
0229          -static_cast<T>(0.027965493064761785),
0230          -static_cast<T>(0.020644781177568105),
0231          -static_cast<T>(0.016650786739707238),
0232          -static_cast<T>(0.014261960828842520),
0233          -static_cast<T>(0.012759847429264803),
0234          -static_cast<T>(0.011799303775587354),
0235          -static_cast<T>(0.011197445703074968)
0236       };
0237       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.15));
0238    }
0239    case 4:
0240    case 5:
0241    //else if (m < 0.3)
0242    {
0243       constexpr T coef[] =
0244       {
0245          static_cast<T>(1.467462209339427155),
0246          -static_cast<T>(0.436576290946337775),
0247          -static_cast<T>(0.105155557666942554),
0248          -static_cast<T>(0.057371843593241730),
0249          -static_cast<T>(0.041391627727340220),
0250          -static_cast<T>(0.034527728505280841),
0251          -static_cast<T>(0.031495443512532783),
0252          -static_cast<T>(0.030527000890325277),
0253          -static_cast<T>(0.030916984019238900),
0254          -static_cast<T>(0.032371395314758122),
0255          -static_cast<T>(0.034789960386404158)
0256       };
0257       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.25));
0258    }
0259    case 6:
0260    case 7:
0261    //else if (m < 0.4)
0262    {
0263       constexpr T coef[] =
0264       {
0265          static_cast<T>(1.422691133490879171),
0266          -static_cast<T>(0.459513519621048674),
0267          -static_cast<T>(0.125250539822061878),
0268          -static_cast<T>(0.078138545094409477),
0269          -static_cast<T>(0.064714278472050002),
0270          -static_cast<T>(0.062084339131730311),
0271          -static_cast<T>(0.065197032815572477),
0272          -static_cast<T>(0.072793895362578779),
0273          -static_cast<T>(0.084959075171781003),
0274          -static_cast<T>(0.102539850131045997),
0275          -static_cast<T>(0.127053585157696036),
0276          -static_cast<T>(0.160791120691274606)
0277       };
0278       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.35));
0279    }
0280    case 8:
0281    case 9:
0282    //else if (m < 0.5)
0283    {
0284       constexpr T coef[] =
0285       {
0286          static_cast<T>(1.375401971871116291),
0287          -static_cast<T>(0.487202183273184837),
0288          -static_cast<T>(0.153311701348540228),
0289          -static_cast<T>(0.111849444917027833),
0290          -static_cast<T>(0.108840952523135768),
0291          -static_cast<T>(0.122954223120269076),
0292          -static_cast<T>(0.152217163962035047),
0293          -static_cast<T>(0.200495323642697339),
0294          -static_cast<T>(0.276174333067751758),
0295          -static_cast<T>(0.393513114304375851),
0296          -static_cast<T>(0.575754406027879147),
0297          -static_cast<T>(0.860523235727239756),
0298          -static_cast<T>(1.308833205758540162)
0299       };
0300       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.45));
0301    }
0302    case 10:
0303    case 11:
0304    //else if (m < 0.6)
0305    {
0306       constexpr T coef[] =
0307       {
0308          static_cast<T>(1.325024497958230082),
0309          -static_cast<T>(0.521727647557566767),
0310          -static_cast<T>(0.194906430482126213),
0311          -static_cast<T>(0.171623726822011264),
0312          -static_cast<T>(0.202754652926419141),
0313          -static_cast<T>(0.278798953118534762),
0314          -static_cast<T>(0.420698457281005762),
0315          -static_cast<T>(0.675948400853106021),
0316          -static_cast<T>(1.136343121839229244),
0317          -static_cast<T>(1.976721143954398261),
0318          -static_cast<T>(3.531696773095722506),
0319          -static_cast<T>(6.446753640156048150),
0320          -static_cast<T>(11.97703130208884026)
0321       };
0322       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.55));
0323    }
0324    case 12:
0325    case 13:
0326    //else if (m < 0.7)
0327    {
0328       constexpr T coef[] =
0329       {
0330          static_cast<T>(1.270707479650149744),
0331          -static_cast<T>(0.566839168287866583),
0332          -static_cast<T>(0.262160793432492598),
0333          -static_cast<T>(0.292244173533077419),
0334          -static_cast<T>(0.440397840850423189),
0335          -static_cast<T>(0.774947641381397458),
0336          -static_cast<T>(1.498870837987561088),
0337          -static_cast<T>(3.089708310445186667),
0338          -static_cast<T>(6.667595903381001064),
0339          -static_cast<T>(14.89436036517319078),
0340          -static_cast<T>(34.18120574251449024),
0341          -static_cast<T>(80.15895841905397306),
0342          -static_cast<T>(191.3489480762984920),
0343          -static_cast<T>(463.5938853480342030),
0344          -static_cast<T>(1137.380822169360061)
0345       };
0346       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.65));
0347    }
0348    case 14:
0349    case 15:
0350    //else if (m < 0.8)
0351    {
0352       constexpr T coef[] =
0353       {
0354          static_cast<T>(1.211056027568459525),
0355          -static_cast<T>(0.630306413287455807),
0356          -static_cast<T>(0.387166409520669145),
0357          -static_cast<T>(0.592278235311934603),
0358          -static_cast<T>(1.237555584513049844),
0359          -static_cast<T>(3.032056661745247199),
0360          -static_cast<T>(8.181688221573590762),
0361          -static_cast<T>(23.55507217389693250),
0362          -static_cast<T>(71.04099935893064956),
0363          -static_cast<T>(221.8796853192349888),
0364          -static_cast<T>(712.1364793277635425),
0365          -static_cast<T>(2336.125331440396407),
0366          -static_cast<T>(7801.945954775964673),
0367          -static_cast<T>(26448.19586059191933),
0368          -static_cast<T>(90799.48341621365251),
0369          -static_cast<T>(315126.0406449163424),
0370          -static_cast<T>(1104011.344311591159)
0371       };
0372       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.75));
0373    }
0374    case 16:
0375    //else if (m < 0.85)
0376    {
0377       constexpr T coef[] =
0378       {
0379          static_cast<T>(1.161307152196282836),
0380          -static_cast<T>(0.701100284555289548),
0381          -static_cast<T>(0.580551474465437362),
0382          -static_cast<T>(1.243693061077786614),
0383          -static_cast<T>(3.679383613496634879),
0384          -static_cast<T>(12.81590924337895775),
0385          -static_cast<T>(49.25672530759985272),
0386          -static_cast<T>(202.1818735434090269),
0387          -static_cast<T>(869.8602699308701437),
0388          -static_cast<T>(3877.005847313289571),
0389          -static_cast<T>(17761.70710170939814),
0390          -static_cast<T>(83182.69029154232061),
0391          -static_cast<T>(396650.4505013548170),
0392          -static_cast<T>(1920033.413682634405)
0393       };
0394       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.825));
0395    }
0396    case 17:
0397    //else if (m < 0.90)
0398    {
0399       constexpr T coef[] =
0400       {
0401          static_cast<T>(1.124617325119752213),
0402          -static_cast<T>(0.770845056360909542),
0403          -static_cast<T>(0.844794053644911362),
0404          -static_cast<T>(2.490097309450394453),
0405          -static_cast<T>(10.23971741154384360),
0406          -static_cast<T>(49.74900546551479866),
0407          -static_cast<T>(267.0986675195705196),
0408          -static_cast<T>(1532.665883825229947),
0409          -static_cast<T>(9222.313478526091951),
0410          -static_cast<T>(57502.51612140314030),
0411          -static_cast<T>(368596.1167416106063),
0412          -static_cast<T>(2415611.088701091428),
0413          -static_cast<T>(16120097.81581656797),
0414          -static_cast<T>(109209938.5203089915),
0415          -static_cast<T>(749380758.1942496220),
0416          -static_cast<T>(5198725846.725541393),
0417          -static_cast<T>(36409256888.12139973)
0418       };
0419       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.875));
0420    }
0421    default:
0422       //
0423       // All cases where m > 0.9
0424       // including all error handling:
0425       //
0426       return ellint_e_imp(k, pol, std::integral_constant<int, 2>());
0427    }
0428 }
0429 template <typename T, typename Policy>
0430 BOOST_MATH_FORCEINLINE T ellint_e_imp(T k, const Policy& pol, std::integral_constant<int, 1> const&)
0431 {
0432    using std::abs;
0433    using namespace boost::math::tools;
0434
0435    T m = k * k;
0436    switch (static_cast<int>(20 * m))
0437    {
0438    case 0:
0439    case 1:
0440       //if (m < 0.1)
0441    {
0442       constexpr T coef[] =
0443       {
0444          1.5509733517804723277L,
0445          -0.40030102010319852390L,
0446          -0.078498619442941939212L,
0447          -0.034318853117591992417L,
0448          -0.019718043317365499309L,
0449          -0.013059507731993309191L,
0450          -0.0094423728741465473894L,
0451          -0.0072467285124021568126L,
0452          -0.0058074240129560897940L,
0453          -0.0048091877860093381762L,
0454          -0.0040863992332551506768L,
0455          -0.0035450302604139562644L,
0456          -0.0031283511188028336315L
0457       };
0458       return boost::math::tools::evaluate_polynomial(coef, m - 0.05L);
0459    }
0460    case 2:
0461    case 3:
0462       //else if (m < 0.2)
0463    {
0464       constexpr T coef[] =
0465       {
0466          1.5101218320928197276L,
0467          -0.41711633390586754922L,
0468          -0.090123820404774568894L,
0469          -0.043729944019084311555L,
0470          -0.027965493064761784548L,
0471          -0.020644781177568105268L,
0472          -0.016650786739707238037L,
0473          -0.014261960828842519634L,
0474          -0.012759847429264802627L,
0475          -0.011799303775587354169L,
0476          -0.011197445703074968018L,
0477          -0.010850368064799902735L,
0478          -0.010696133481060989818L
0479       };
0480       return boost::math::tools::evaluate_polynomial(coef, m - 0.15L);
0481    }
0482    case 4:
0483    case 5:
0484       //else if (m < 0.3L)
0485    {
0486       constexpr T coef[] =
0487       {
0488          1.4674622093394271555L,
0489          -0.43657629094633777482L,
0490          -0.10515555766694255399L,
0491          -0.057371843593241729895L,
0492          -0.041391627727340220236L,
0493          -0.034527728505280841188L,
0494          -0.031495443512532782647L,
0495          -0.030527000890325277179L,
0496          -0.030916984019238900349L,
0497          -0.032371395314758122268L,
0498          -0.034789960386404158240L,
0499          -0.038182654612387881967L,
0500          -0.042636187648900252525L,
0501          -0.048302272505241634467
0502       };
0503       return boost::math::tools::evaluate_polynomial(coef, m - 0.25L);
0504    }
0505    case 6:
0506    case 7:
0507       //else if (m < 0.4L)
0508    {
0509       constexpr T coef[] =
0510       {
0511          1.4226911334908791711L,
0512          -0.45951351962104867394L,
0513          -0.12525053982206187849L,
0514          -0.078138545094409477156L,
0515          -0.064714278472050001838L,
0516          -0.062084339131730310707L,
0517          -0.065197032815572476910L,
0518          -0.072793895362578779473L,
0519          -0.084959075171781003264L,
0520          -0.10253985013104599679L,
0521          -0.12705358515769603644L,
0522          -0.16079112069127460621L,
0523          -0.20705400012405941376L,
0524          -0.27053164884730888948L
0525       };
0526       return boost::math::tools::evaluate_polynomial(coef, m - 0.35L);
0527    }
0528    case 8:
0529    case 9:
0530       //else if (m < 0.5L)
0531    {
0532       constexpr T coef[] =
0533       {
0534          1.3754019718711162908L,
0535          -0.48720218327318483652L,
0536          -0.15331170134854022753L,
0537          -0.11184944491702783273L,
0538          -0.10884095252313576755L,
0539          -0.12295422312026907610L,
0540          -0.15221716396203504746L,
0541          -0.20049532364269733857L,
0542          -0.27617433306775175837L,
0543          -0.39351311430437585139L,
0544          -0.57575440602787914711L,
0545          -0.86052323572723975634L,
0546          -1.3088332057585401616L,
0547          -2.0200280559452241745L,
0548          -3.1566019548237606451L
0549       };
0550       return boost::math::tools::evaluate_polynomial(coef, m - 0.45L);
0551    }
0552    case 10:
0553    case 11:
0554       //else if (m < 0.6L)
0555    {
0556       constexpr T coef[] =
0557       {
0558          1.3250244979582300818L,
0559          -0.52172764755756676713L,
0560          -0.19490643048212621262L,
0561          -0.17162372682201126365L,
0562          -0.20275465292641914128L,
0563          -0.27879895311853476205L,
0564          -0.42069845728100576224L,
0565          -0.67594840085310602110L,
0566          -1.1363431218392292440L,
0567          -1.9767211439543982613L,
0568          -3.5316967730957225064L,
0569          -6.4467536401560481499L,
0570          -11.977031302088840261L,
0571          -22.581360948073964469L,
0572          -43.109479829481450573L,
0573          -83.186290908288807424L
0574       };
0575       return boost::math::tools::evaluate_polynomial(coef, m - 0.55L);
0576    }
0577    case 12:
0578    case 13:
0579       //else if (m < 0.7L)
0580    {
0581       constexpr T coef[] =
0582       {
0583          1.2707074796501497440L,
0584          -0.56683916828786658286L,
0585          -0.26216079343249259779L,
0586          -0.29224417353307741931L,
0587          -0.44039784085042318909L,
0588          -0.77494764138139745824L,
0589          -1.4988708379875610880L,
0590          -3.0897083104451866665L,
0591          -6.6675959033810010645L,
0592          -14.894360365173190775L,
0593          -34.181205742514490240L,
0594          -80.158958419053973056L,
0595          -191.34894807629849204L,
0596          -463.59388534803420301L,
0597          -1137.3808221693600606L,
0598          -2820.7073786352269339L,
0599          -7061.1382244658715621L,
0600          -17821.809331816437058L,
0601          -45307.849987201897801L
0602       };
0603       return boost::math::tools::evaluate_polynomial(coef, m - 0.65L);
0604    }
0605    case 14:
0606    case 15:
0607       //else if (m < 0.8L)
0608    {
0609       constexpr T coef[] =
0610       {
0611          1.2110560275684595248L,
0612          -0.63030641328745580709L,
0613          -0.38716640952066914514L,
0614          -0.59227823531193460257L,
0615          -1.2375555845130498445L,
0616          -3.0320566617452471986L,
0617          -8.1816882215735907624L,
0618          -23.555072173896932503L,
0619          -71.040999358930649565L,
0620          -221.87968531923498875L,
0621          -712.13647932776354253L,
0622          -2336.1253314403964072L,
0623          -7801.9459547759646726L,
0624          -26448.195860591919335L,
0625          -90799.483416213652512L,
0626          -315126.04064491634241L,
0627          -1.1040113443115911589e6L,
0628          -3.8998018348056769095e6L,
0629          -1.3876249116223745041e7L,
0630          -4.9694982823537861149e7L,
0631          -1.7900668836197342979e8L,
0632          -6.4817399873722371964e8L
0633       };
0634       return boost::math::tools::evaluate_polynomial(coef, m - 0.75L);
0635    }
0636    case 16:
0637       //else if (m < 0.85L)
0638    {
0639       constexpr T coef[] =
0640       {
0641          1.1613071521962828360L,
0642          -0.70110028455528954752L,
0643          -0.58055147446543736163L,
0644          -1.2436930610777866138L,
0645          -3.6793836134966348789L,
0646          -12.815909243378957753L,
0647          -49.256725307599852720L,
0648          -202.18187354340902693L,
0649          -869.86026993087014372L,
0650          -3877.0058473132895713L,
0651          -17761.707101709398174L,
0652          -83182.690291542320614L,
0653          -396650.45050135481698L,
0654          -1.9200334136826344054e6L,
0655          -9.4131321779500838352e6L,
0656          -4.6654858837335370627e7L,
0657          -2.3343549352617609390e8L,
0658          -1.1776928223045913454e9L,
0659          -5.9850851892915740401e9L,
0660          -3.0614702984618644983e10L
0661       };
0662       return boost::math::tools::evaluate_polynomial(coef, m - 0.825L);
0663    }
0664    case 17:
0665       //else if (m < 0.90L)
0666    {
0667       constexpr T coef[] =
0668       {
0669          1.1246173251197522132L,
0670          -0.77084505636090954218L,
0671          -0.84479405364491136236L,
0672          -2.4900973094503944527L,
0673          -10.239717411543843601L,
0674          -49.749005465514798660L,
0675          -267.09866751957051961L,
0676          -1532.6658838252299468L,
0677          -9222.3134785260919507L,
0678          -57502.516121403140303L,
0679          -368596.11674161060626L,
0680          -2.4156110887010914281e6L,
0681          -1.6120097815816567971e7L,
0682          -1.0920993852030899148e8L,
0683          -7.4938075819424962198e8L,
0684          -5.1987258467255413931e9L,
0685          -3.6409256888121399726e10L,
0686          -2.5711802891217393544e11L,
0687          -1.8290904062978796996e12L,
0688          -1.3096838781743248404e13L,
0689          -9.4325465851415135118e13L,
0690          -6.8291980829471896669e14L
0691       };
0692       return boost::math::tools::evaluate_polynomial(coef, m - 0.875L);
0693    }
0694    default:
0695       //
0696       // All cases where m > 0.9
0697       // including all error handling:
0698       //
0699       return ellint_e_imp(k, pol, std::integral_constant<int, 2>());
0700    }
0701 }
0702
0703 template <typename T, typename Policy>
0704 BOOST_MATH_FORCEINLINE typename tools::promote_args<T>::type ellint_2(T k, const Policy& pol, const std::true_type&)
0705 {
0706    typedef typename tools::promote_args<T>::type result_type;
0707    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0708    typedef std::integral_constant<int,
0709       std::is_floating_point<T>::value&& std::numeric_limits<T>::digits && (std::numeric_limits<T>::digits <= 54) ? 0 :
0710       std::is_floating_point<T>::value && std::numeric_limits<T>::digits && (std::numeric_limits<T>::digits <= 64) ? 1 : 2
0711    > precision_tag_type;
0712    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%)");
0713 }
0714
0715 // Elliptic integral (Legendre form) of the second kind
0716 template <class T1, class T2>
0717 BOOST_MATH_FORCEINLINE typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi, const std::false_type&)
0718 {
0719    return boost::math::ellint_2(k, phi, policies::policy<>());
0720 }
0721
0722 } // detail
0723
0724 // Complete elliptic integral (Legendre form) of the second kind
0725 template <typename T>
0726 BOOST_MATH_FORCEINLINE typename tools::promote_args<T>::type ellint_2(T k)
0727 {
0728    return ellint_2(k, policies::policy<>());
0729 }
0730
0731 // Elliptic integral (Legendre form) of the second kind
0732 template <class T1, class T2>
0733 BOOST_MATH_FORCEINLINE typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi)
0734 {
0735    typedef typename policies::is_policy<T2>::type tag_type;
0736    return detail::ellint_2(k, phi, tag_type());
0737 }
0738
0739 template <class T1, class T2, class Policy>
0740 BOOST_MATH_FORCEINLINE typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi, const Policy& pol)  // LCOV_EXCL_LINE gcc misses this but sees the function body, strange!
0741 {
0742    typedef typename tools::promote_args<T1, T2>::type result_type;
0743    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0744    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%)");
0745 }
0746
0747 }} // namespaces
0748
0749 #endif // BOOST_MATH_ELLINT_2_HPP
0750