math/special_functions/ellint_1.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_1_HPP
0015 #define BOOST_MATH_ELLINT_1_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/constants/constants.hpp>
0024 #include <boost/math/policies/error_handling.hpp>
0025 #include <boost/math/tools/workaround.hpp>
0026 #include <boost/math/special_functions/round.hpp>
0027
0028 // Elliptic integrals (complete and incomplete) of the first kind
0029 // Carlson, Numerische Mathematik, vol 33, 1 (1979)
0030
0031 namespace boost { namespace math {
0032
0033 template <class T1, class T2, class Policy>
0034 typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi, const Policy& pol);
0035
0036 namespace detail{
0037
0038 template <typename T, typename Policy>
0039 T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 0> const&);
0040 template <typename T, typename Policy>
0041 T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 1> const&);
0042 template <typename T, typename Policy>
0043 T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 2> const&);
0044 template <typename T, typename Policy>
0045 T ellint_k_imp(T k, const Policy& pol, T one_minus_k2);
0046
0047 // Elliptic integral (Legendre form) of the first kind
0048 template <typename T, typename Policy>
0049 T ellint_f_imp(T phi, T k, const Policy& pol, T one_minus_k2)
0050 {
0051     BOOST_MATH_STD_USING
0052     using namespace boost::math::tools;
0053     using namespace boost::math::constants;
0054
0055     static const char* function = "boost::math::ellint_f<%1%>(%1%,%1%)";
0056     BOOST_MATH_INSTRUMENT_VARIABLE(phi);
0057     BOOST_MATH_INSTRUMENT_VARIABLE(k);
0058     BOOST_MATH_INSTRUMENT_VARIABLE(function);
0059
0060     bool invert = false;
0061     if(phi < 0)
0062     {
0063        BOOST_MATH_INSTRUMENT_VARIABLE(phi);
0064        phi = fabs(phi);
0065        invert = true;
0066     }
0067
0068     T result;
0069
0070     if(phi >= tools::max_value<T>())
0071     {
0072        // Need to handle infinity as a special case:
0073        result = policies::raise_overflow_error<T>(function, nullptr, pol);
0074        BOOST_MATH_INSTRUMENT_VARIABLE(result);
0075     }
0076     else if(phi > 1 / tools::epsilon<T>())
0077     {
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_k_imp(k, pol, one_minus_k2) / constants::pi<T>();
0081        BOOST_MATH_INSTRUMENT_VARIABLE(result);
0082     }
0083     else
0084     {
0085        // Carlson's algorithm works only for |phi| <= pi/2,
0086        // use the integrand's periodicity to normalize phi
0087        //
0088        // Xiaogang's original code used a cast to long long here
0089        // but that fails if T has more digits than a long long,
0090        // so rewritten to use fmod instead:
0091        //
0092        BOOST_MATH_INSTRUMENT_CODE("pi/2 = " << constants::pi<T>() / 2);
0093        T rphi = boost::math::tools::fmod_workaround(phi, T(constants::half_pi<T>()));
0094        BOOST_MATH_INSTRUMENT_VARIABLE(rphi);
0095        T m = boost::math::round((phi - rphi) / constants::half_pi<T>());
0096        BOOST_MATH_INSTRUMENT_VARIABLE(m);
0097        int s = 1;
0098        if(boost::math::tools::fmod_workaround(m, T(2)) > T(0.5))
0099        {
0100           m += 1;
0101           s = -1;
0102           rphi = constants::half_pi<T>() - rphi;
0103           BOOST_MATH_INSTRUMENT_VARIABLE(rphi);
0104        }
0105        T sinp = sin(rphi);
0106        sinp *= sinp;
0107        if (sinp * k * k >= 1)
0108        {
0109           return policies::raise_domain_error<T>(function,
0110              "Got k^2 * sin^2(phi) = %1%, but the function requires this < 1", sinp * k * k, pol);
0111        }
0112        T cosp = cos(rphi);
0113        cosp *= cosp;
0114        BOOST_MATH_INSTRUMENT_VARIABLE(sinp);
0115        BOOST_MATH_INSTRUMENT_VARIABLE(cosp);
0116        if(sinp > tools::min_value<T>())
0117        {
0118           BOOST_MATH_ASSERT(rphi != 0); // precondition, can't be true if sin(rphi) != 0.
0119           //
0120           // Use http://dlmf.nist.gov/19.25#E5, note that
0121           // c-1 simplifies to cot^2(rphi) which avoids cancellation.
0122           // Likewise c - k^2 is the same as (c - 1) + (1 - k^2).
0123           //
0124           T c = 1 / sinp;
0125           T c_minus_one = cosp / sinp;
0126           T arg2;
0127           if (k != 0)
0128           {
0129              T cross = fabs(c / (k * k));
0130              if ((cross > 0.9f) && (cross < 1.1f))
0131                 arg2 = c_minus_one + one_minus_k2;
0132              else
0133                 arg2 = c - k * k;
0134           }
0135           else
0136              arg2 = c;
0137           result = static_cast<T>(s * ellint_rf_imp(c_minus_one, arg2, c, pol));
0138        }
0139        else
0140           result = s * sin(rphi);
0141        BOOST_MATH_INSTRUMENT_VARIABLE(result);
0142        if(m != 0)
0143        {
0144           result += m * ellint_k_imp(k, pol, one_minus_k2);
0145           BOOST_MATH_INSTRUMENT_VARIABLE(result);
0146        }
0147     }
0148     return invert ? T(-result) : result;
0149 }
0150
0151 template <typename T, typename Policy>
0152 inline T ellint_f_imp(T phi, T k, const Policy& pol)
0153 {
0154    return ellint_f_imp(phi, k, pol, T(1 - k * k));
0155 }
0156
0157 // Complete elliptic integral (Legendre form) of the first kind
0158 template <typename T, typename Policy>
0159 T ellint_k_imp(T k, const Policy& pol, T one_minus_k2)
0160 {
0161     BOOST_MATH_STD_USING
0162     using namespace boost::math::tools;
0163
0164     static const char* function = "boost::math::ellint_k<%1%>(%1%)";
0165
0166     if (abs(k) > 1)
0167     {
0168        return policies::raise_domain_error<T>(function, "Got k = %1%, function requires |k| <= 1", k, pol);
0169     }
0170     if (abs(k) == 1)
0171     {
0172        return policies::raise_overflow_error<T>(function, nullptr, pol);
0173     }
0174
0175     T x = 0;
0176     T z = 1;
0177     T value = ellint_rf_imp(x, one_minus_k2, z, pol);
0178
0179     return value;
0180 }
0181 template <typename T, typename Policy>
0182 inline T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 2> const&)
0183 {
0184    return ellint_k_imp(k, pol, T(1 - k * k));
0185 }
0186
0187 //
0188 // Special versions for double and 80-bit long double precision,
0189 // double precision versions use the coefficients from:
0190 // "Fast computation of complete elliptic integrals and Jacobian elliptic functions",
0191 // Celestial Mechanics and Dynamical Astronomy, April 2012.
0192 //
0193 // Higher precision coefficients for 80-bit long doubles can be calculated
0194 // using for example:
0195 // Table[N[SeriesCoefficient[ EllipticK [ m ], { m, 875/1000, i} ], 20], {i, 0, 24}]
0196 // and checking the value of the first neglected term with:
0197 // N[SeriesCoefficient[ EllipticK [ m ], { m, 875/1000, 24} ], 20] * (2.5/100)^24
0198 //
0199 // For m > 0.9 we don't use the method of the paper above, but simply call our
0200 // existing routines.  The routine used in the above paper was tried (and is
0201 // archived in the code below), but was found to have slightly higher error rates.
0202 //
0203 template <typename T, typename Policy>
0204 BOOST_MATH_FORCEINLINE T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 0> const&)
0205 {
0206    using std::abs;
0207    using namespace boost::math::tools;
0208
0209    T m = k * k;
0210
0211    switch (static_cast<int>(m * 20))
0212    {
0213    case 0:
0214    case 1:
0215       //if (m < 0.1)
0216    {
0217       constexpr T coef[] =
0218       {
0219          static_cast<T>(1.591003453790792180),
0220          static_cast<T>(0.416000743991786912),
0221          static_cast<T>(0.245791514264103415),
0222          static_cast<T>(0.179481482914906162),
0223          static_cast<T>(0.144556057087555150),
0224          static_cast<T>(0.123200993312427711),
0225          static_cast<T>(0.108938811574293531),
0226          static_cast<T>(0.098853409871592910),
0227          static_cast<T>(0.091439629201749751),
0228          static_cast<T>(0.085842591595413900),
0229          static_cast<T>(0.081541118718303215),
0230          static_cast<T>(0.078199656811256481910)
0231       };
0232       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.05));
0233    }
0234    case 2:
0235    case 3:
0236       //else if (m < 0.2)
0237    {
0238       constexpr T coef[] =
0239       {
0240          static_cast<T>(1.635256732264579992),
0241          static_cast<T>(0.471190626148732291),
0242          static_cast<T>(0.309728410831499587),
0243          static_cast<T>(0.252208311773135699),
0244          static_cast<T>(0.226725623219684650),
0245          static_cast<T>(0.215774446729585976),
0246          static_cast<T>(0.213108771877348910),
0247          static_cast<T>(0.216029124605188282),
0248          static_cast<T>(0.223255831633057896),
0249          static_cast<T>(0.234180501294209925),
0250          static_cast<T>(0.248557682972264071),
0251          static_cast<T>(0.266363809892617521)
0252       };
0253       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.15));
0254    }
0255    case 4:
0256    case 5:
0257       //else if (m < 0.3)
0258    {
0259       constexpr T coef[] =
0260       {
0261          static_cast<T>(1.685750354812596043),
0262          static_cast<T>(0.541731848613280329),
0263          static_cast<T>(0.401524438390690257),
0264          static_cast<T>(0.369642473420889090),
0265          static_cast<T>(0.376060715354583645),
0266          static_cast<T>(0.405235887085125919),
0267          static_cast<T>(0.453294381753999079),
0268          static_cast<T>(0.520518947651184205),
0269          static_cast<T>(0.609426039204995055),
0270          static_cast<T>(0.724263522282908870),
0271          static_cast<T>(0.871013847709812357),
0272          static_cast<T>(1.057652872753547036)
0273       };
0274       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.25));
0275    }
0276    case 6:
0277    case 7:
0278       //else if (m < 0.4)
0279    {
0280       constexpr T coef[] =
0281       {
0282          static_cast<T>(1.744350597225613243),
0283          static_cast<T>(0.634864275371935304),
0284          static_cast<T>(0.539842564164445538),
0285          static_cast<T>(0.571892705193787391),
0286          static_cast<T>(0.670295136265406100),
0287          static_cast<T>(0.832586590010977199),
0288          static_cast<T>(1.073857448247933265),
0289          static_cast<T>(1.422091460675497751),
0290          static_cast<T>(1.920387183402304829),
0291          static_cast<T>(2.632552548331654201),
0292          static_cast<T>(3.652109747319039160),
0293          static_cast<T>(5.115867135558865806),
0294          static_cast<T>(7.224080007363877411)
0295       };
0296       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.35));
0297    }
0298    case 8:
0299    case 9:
0300       //else if (m < 0.5)
0301    {
0302       constexpr T coef[] =
0303       {
0304          static_cast<T>(1.813883936816982644),
0305          static_cast<T>(0.763163245700557246),
0306          static_cast<T>(0.761928605321595831),
0307          static_cast<T>(0.951074653668427927),
0308          static_cast<T>(1.315180671703161215),
0309          static_cast<T>(1.928560693477410941),
0310          static_cast<T>(2.937509342531378755),
0311          static_cast<T>(4.594894405442878062),
0312          static_cast<T>(7.330071221881720772),
0313          static_cast<T>(11.87151259742530180),
0314          static_cast<T>(19.45851374822937738),
0315          static_cast<T>(32.20638657246426863),
0316          static_cast<T>(53.73749198700554656),
0317          static_cast<T>(90.27388602940998849)
0318       };
0319       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.45));
0320    }
0321    case 10:
0322    case 11:
0323       //else if (m < 0.6)
0324    {
0325       constexpr T coef[] =
0326       {
0327          static_cast<T>(1.898924910271553526),
0328          static_cast<T>(0.950521794618244435),
0329          static_cast<T>(1.151077589959015808),
0330          static_cast<T>(1.750239106986300540),
0331          static_cast<T>(2.952676812636875180),
0332          static_cast<T>(5.285800396121450889),
0333          static_cast<T>(9.832485716659979747),
0334          static_cast<T>(18.78714868327559562),
0335          static_cast<T>(36.61468615273698145),
0336          static_cast<T>(72.45292395127771801),
0337          static_cast<T>(145.1079577347069102),
0338          static_cast<T>(293.4786396308497026),
0339          static_cast<T>(598.3851815055010179),
0340          static_cast<T>(1228.420013075863451),
0341          static_cast<T>(2536.529755382764488)
0342       };
0343       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.55));
0344    }
0345    case 12:
0346    case 13:
0347       //else if (m < 0.7)
0348    {
0349       constexpr T coef[] =
0350       {
0351          static_cast<T>(2.007598398424376302),
0352          static_cast<T>(1.248457231212347337),
0353          static_cast<T>(1.926234657076479729),
0354          static_cast<T>(3.751289640087587680),
0355          static_cast<T>(8.119944554932045802),
0356          static_cast<T>(18.66572130873555361),
0357          static_cast<T>(44.60392484291437063),
0358          static_cast<T>(109.5092054309498377),
0359          static_cast<T>(274.2779548232413480),
0360          static_cast<T>(697.5598008606326163),
0361          static_cast<T>(1795.716014500247129),
0362          static_cast<T>(4668.381716790389910),
0363          static_cast<T>(12235.76246813664335),
0364          static_cast<T>(32290.17809718320818),
0365          static_cast<T>(85713.07608195964685),
0366          static_cast<T>(228672.1890493117096),
0367          static_cast<T>(612757.2711915852774)
0368       };
0369       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.65));
0370    }
0371    case 14:
0372    case 15:
0373       //else if (m < static_cast<T>(0.8))
0374    {
0375       constexpr T coef[] =
0376       {
0377          static_cast<T>(2.156515647499643235),
0378          static_cast<T>(1.791805641849463243),
0379          static_cast<T>(3.826751287465713147),
0380          static_cast<T>(10.38672468363797208),
0381          static_cast<T>(31.40331405468070290),
0382          static_cast<T>(100.9237039498695416),
0383          static_cast<T>(337.3268282632272897),
0384          static_cast<T>(1158.707930567827917),
0385          static_cast<T>(4060.990742193632092),
0386          static_cast<T>(14454.00184034344795),
0387          static_cast<T>(52076.66107599404803),
0388          static_cast<T>(189493.6591462156887),
0389          static_cast<T>(695184.5762413896145),
0390          static_cast<T>(2567994.048255284686),
0391          static_cast<T>(9541921.966748386322),
0392          static_cast<T>(35634927.44218076174),
0393          static_cast<T>(133669298.4612040871),
0394          static_cast<T>(503352186.6866284541),
0395          static_cast<T>(1901975729.538660119),
0396          static_cast<T>(7208915015.330103756)
0397       };
0398       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.75));
0399    }
0400    case 16:
0401       //else if (m < static_cast<T>(0.85))
0402    {
0403       constexpr T coef[] =
0404       {
0405          static_cast<T>(2.318122621712510589),
0406          static_cast<T>(2.616920150291232841),
0407          static_cast<T>(7.897935075731355823),
0408          static_cast<T>(30.50239715446672327),
0409          static_cast<T>(131.4869365523528456),
0410          static_cast<T>(602.9847637356491617),
0411          static_cast<T>(2877.024617809972641),
0412          static_cast<T>(14110.51991915180325),
0413          static_cast<T>(70621.44088156540229),
0414          static_cast<T>(358977.2665825309926),
0415          static_cast<T>(1847238.263723971684),
0416          static_cast<T>(9600515.416049214109),
0417          static_cast<T>(50307677.08502366879),
0418          static_cast<T>(265444188.6527127967),
0419          static_cast<T>(1408862325.028702687),
0420          static_cast<T>(7515687935.373774627)
0421       };
0422       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.825));
0423    }
0424    case 17:
0425       //else if (m < static_cast<T>(0.90))
0426    {
0427       constexpr T coef[] =
0428       {
0429          static_cast<T>(2.473596173751343912),
0430          static_cast<T>(3.727624244118099310),
0431          static_cast<T>(15.60739303554930496),
0432          static_cast<T>(84.12850842805887747),
0433          static_cast<T>(506.9818197040613935),
0434          static_cast<T>(3252.277058145123644),
0435          static_cast<T>(21713.24241957434256),
0436          static_cast<T>(149037.0451890932766),
0437          static_cast<T>(1043999.331089990839),
0438          static_cast<T>(7427974.817042038995),
0439          static_cast<T>(53503839.67558661151),
0440          static_cast<T>(389249886.9948708474),
0441          static_cast<T>(2855288351.100810619),
0442          static_cast<T>(21090077038.76684053),
0443          static_cast<T>(156699833947.7902014),
0444          static_cast<T>(1170222242422.439893),
0445          static_cast<T>(8777948323668.937971),
0446          static_cast<T>(66101242752484.95041),
0447          static_cast<T>(499488053713388.7989),
0448          static_cast<T>(37859743397240299.20)
0449       };
0450       return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.875));
0451    }
0452    default:
0453       //
0454       // This handles all cases where m > 0.9,
0455       // including all error handling:
0456       //
0457       return ellint_k_imp(k, pol, std::integral_constant<int, 2>());
0458 #if 0
0459    else
0460    {
0461       T lambda_prime = (1 - sqrt(k)) / (2 * (1 + sqrt(k)));
0462       T k_prime = ellint_k(sqrt((1 - k) * (1 + k))); // K(m')
0463       T lambda_prime_4th = boost::math::pow<4>(lambda_prime);
0464       T q_prime = ((((((20910 * lambda_prime_4th) + 1707) * lambda_prime_4th + 150) * lambda_prime_4th + 15) * lambda_prime_4th + 2) * lambda_prime_4th + 1) * lambda_prime;
0465       /*T q_prime_2 = lambda_prime
0466          + 2 * boost::math::pow<5>(lambda_prime)
0467          + 15 * boost::math::pow<9>(lambda_prime)
0468          + 150 * boost::math::pow<13>(lambda_prime)
0469          + 1707 * boost::math::pow<17>(lambda_prime)
0470          + 20910 * boost::math::pow<21>(lambda_prime);*/
0471       return -log(q_prime) * k_prime / boost::math::constants::pi<T>();
0472    }
0473 #endif
0474    }
0475 }
0476 template <typename T, typename Policy>
0477 BOOST_MATH_FORCEINLINE T ellint_k_imp(T k, const Policy& pol, std::integral_constant<int, 1> const&)
0478 {
0479    using std::abs;
0480    using namespace boost::math::tools;
0481
0482    T m = k * k;
0483    switch (static_cast<int>(m * 20))
0484    {
0485    case 0:
0486    case 1:
0487    {
0488       constexpr T coef[] =
0489       {
0490          1.5910034537907921801L,
0491          0.41600074399178691174L,
0492          0.24579151426410341536L,
0493          0.17948148291490616181L,
0494          0.14455605708755514976L,
0495          0.12320099331242771115L,
0496          0.10893881157429353105L,
0497          0.098853409871592910399L,
0498          0.091439629201749751268L,
0499          0.085842591595413899672L,
0500          0.081541118718303214749L,
0501          0.078199656811256481910L,
0502          0.075592617535422415648L,
0503          0.073562939365441925050L
0504       };
0505       return boost::math::tools::evaluate_polynomial(coef, m - 0.05L);
0506    }
0507    case 2:
0508    case 3:
0509    {
0510       constexpr T coef[] =
0511       {
0512          1.6352567322645799924L,
0513          0.47119062614873229055L,
0514          0.30972841083149958708L,
0515          0.25220831177313569923L,
0516          0.22672562321968464974L,
0517          0.21577444672958597588L,
0518          0.21310877187734890963L,
0519          0.21602912460518828154L,
0520          0.22325583163305789567L,
0521          0.23418050129420992492L,
0522          0.24855768297226407136L,
0523          0.26636380989261752077L,
0524          0.28772845215611466775L,
0525          0.31290024539780334906L,
0526          0.34223105446381299902L
0527       };
0528       return boost::math::tools::evaluate_polynomial(coef, m - 0.15L);
0529    }
0530    case 4:
0531    case 5:
0532    {
0533       constexpr T coef[] =
0534       {
0535          1.6857503548125960429L,
0536          0.54173184861328032882L,
0537          0.40152443839069025682L,
0538          0.36964247342088908995L,
0539          0.37606071535458364462L,
0540          0.40523588708512591863L,
0541          0.45329438175399907924L,
0542          0.52051894765118420473L,
0543          0.60942603920499505544L,
0544          0.72426352228290886975L,
0545          0.87101384770981235737L,
0546          1.0576528727535470365L,
0547          1.2945970872087764321L,
0548          1.5953368253888783747L,
0549          1.9772844873556364793L,
0550          2.4628890581910021287L
0551       };
0552       return boost::math::tools::evaluate_polynomial(coef, m - 0.25L);
0553    }
0554    case 6:
0555    case 7:
0556    {
0557       constexpr T coef[] =
0558       {
0559          1.7443505972256132429L,
0560          0.63486427537193530383L,
0561          0.53984256416444553751L,
0562          0.57189270519378739093L,
0563          0.67029513626540610034L,
0564          0.83258659001097719939L,
0565          1.0738574482479332654L,
0566          1.4220914606754977514L,
0567          1.9203871834023048288L,
0568          2.6325525483316542006L,
0569          3.6521097473190391602L,
0570          5.1158671355588658061L,
0571          7.2240800073638774108L,
0572          10.270306349944787227L,
0573          14.685616935355757348L,
0574          21.104114212004582734L,
0575          30.460132808575799413L,
0576       };
0577       return boost::math::tools::evaluate_polynomial(coef, m - 0.35L);
0578    }
0579    case 8:
0580    case 9:
0581    {
0582       constexpr T coef[] =
0583       {
0584          1.8138839368169826437L,
0585          0.76316324570055724607L,
0586          0.76192860532159583095L,
0587          0.95107465366842792679L,
0588          1.3151806717031612153L,
0589          1.9285606934774109412L,
0590          2.9375093425313787550L,
0591          4.5948944054428780618L,
0592          7.3300712218817207718L,
0593          11.871512597425301798L,
0594          19.458513748229377383L,
0595          32.206386572464268628L,
0596          53.737491987005546559L,
0597          90.273886029409988491L,
0598          152.53312130253275268L,
0599          259.02388747148299086L,
0600          441.78537518096201946L,
0601          756.39903981567380952L
0602       };
0603       return boost::math::tools::evaluate_polynomial(coef, m - 0.45L);
0604    }
0605    case 10:
0606    case 11:
0607    {
0608       constexpr T coef[] =
0609       {
0610          1.8989249102715535257L,
0611          0.95052179461824443490L,
0612          1.1510775899590158079L,
0613          1.7502391069863005399L,
0614          2.9526768126368751802L,
0615          5.2858003961214508892L,
0616          9.8324857166599797471L,
0617          18.787148683275595622L,
0618          36.614686152736981447L,
0619          72.452923951277718013L,
0620          145.10795773470691023L,
0621          293.47863963084970259L,
0622          598.38518150550101790L,
0623          1228.4200130758634505L,
0624          2536.5297553827644880L,
0625          5263.9832725075189576L,
0626          10972.138126273491753L,
0627          22958.388550988306870L,
0628          48203.103373625406989L
0629       };
0630       return boost::math::tools::evaluate_polynomial(coef, m - 0.55L);
0631    }
0632    case 12:
0633    case 13:
0634    {
0635       constexpr T coef[] =
0636       {
0637          2.0075983984243763017L,
0638          1.2484572312123473371L,
0639          1.9262346570764797287L,
0640          3.7512896400875876798L,
0641          8.1199445549320458022L,
0642          18.665721308735553611L,
0643          44.603924842914370633L,
0644          109.50920543094983774L,
0645          274.27795482324134804L,
0646          697.55980086063261629L,
0647          1795.7160145002471293L,
0648          4668.3817167903899100L,
0649          12235.762468136643348L,
0650          32290.178097183208178L,
0651          85713.076081959646847L,
0652          228672.18904931170958L,
0653          612757.27119158527740L,
0654          1.6483233976504668314e6L,
0655          4.4492251046211960936e6L,
0656          1.2046317340783185238e7L,
0657          3.2705187507963254185e7L
0658       };
0659       return boost::math::tools::evaluate_polynomial(coef, m - 0.65L);
0660    }
0661    case 14:
0662    case 15:
0663    {
0664       constexpr T coef[] =
0665       {
0666          2.1565156474996432354L,
0667          1.7918056418494632425L,
0668          3.8267512874657131470L,
0669          10.386724683637972080L,
0670          31.403314054680702901L,
0671          100.92370394986954165L,
0672          337.32682826322728966L,
0673          1158.7079305678279173L,
0674          4060.9907421936320917L,
0675          14454.001840343447947L,
0676          52076.661075994048028L,
0677          189493.65914621568866L,
0678          695184.57624138961450L,
0679          2.5679940482552846861e6L,
0680          9.5419219667483863221e6L,
0681          3.5634927442180761743e7L,
0682          1.3366929846120408712e8L,
0683          5.0335218668662845411e8L,
0684          1.9019757295386601192e9L,
0685          7.2089150153301037563e9L,
0686          2.7398741806339510931e10L,
0687          1.0439286724885300495e11L,
0688          3.9864875581513728207e11L,
0689          1.5254661585564745591e12L,
0690          5.8483259088850315936e12
0691       };
0692       return boost::math::tools::evaluate_polynomial(coef, m - 0.75L);
0693    }
0694    case 16:
0695    {
0696       constexpr T coef[] =
0697       {
0698          2.3181226217125105894L,
0699          2.6169201502912328409L,
0700          7.8979350757313558232L,
0701          30.502397154466723270L,
0702          131.48693655235284561L,
0703          602.98476373564916170L,
0704          2877.0246178099726410L,
0705          14110.519919151803247L,
0706          70621.440881565402289L,
0707          358977.26658253099258L,
0708          1.8472382637239716844e6L,
0709          9.6005154160492141090e6L,
0710          5.0307677085023668786e7L,
0711          2.6544418865271279673e8L,
0712          1.4088623250287026866e9L,
0713          7.5156879353737746270e9L,
0714          4.0270783964955246149e10L,
0715          2.1662089325801126339e11L,
0716          1.1692489201929996116e12L,
0717          6.3306543358985679881e12
0718       };
0719       return boost::math::tools::evaluate_polynomial(coef, m - 0.825L);
0720    }
0721    case 17:
0722    {
0723       constexpr T coef[] =
0724       {
0725          2.4735961737513439120L,
0726          3.7276242441180993105L,
0727          15.607393035549304964L,
0728          84.128508428058877470L,
0729          506.98181970406139349L,
0730          3252.2770581451236438L,
0731          21713.242419574342564L,
0732          149037.04518909327662L,
0733          1.0439993310899908390e6L,
0734          7.4279748170420389947e6L,
0735          5.3503839675586611510e7L,
0736          3.8924988699487084738e8L,
0737          2.8552883511008106195e9L,
0738          2.1090077038766840525e10L,
0739          1.5669983394779020136e11L,
0740          1.1702222424224398927e12L,
0741          8.7779483236689379709e12L,
0742          6.6101242752484950408e13L,
0743          4.9948805371338879891e14L,
0744          3.7859743397240299201e15L,
0745          2.8775996123036112296e16L,
0746          2.1926346839925760143e17L,
0747          1.6744985438468349361e18L,
0748          1.2814410112866546052e19L,
0749          9.8249807041031260167e19
0750       };
0751       return boost::math::tools::evaluate_polynomial(coef, m - 0.875L);
0752    }
0753    default:
0754       //
0755       // All cases where m > 0.9
0756       // including all error handling:
0757       //
0758       return ellint_k_imp(k, pol, std::integral_constant<int, 2>());
0759    }
0760 }
0761
0762 template <typename T, typename Policy>
0763 BOOST_MATH_FORCEINLINE typename tools::promote_args<T>::type ellint_1(T k, const Policy& pol, const std::true_type&)
0764 {
0765    typedef typename tools::promote_args<T>::type result_type;
0766    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0767    typedef std::integral_constant<int,
0768 #if defined(__clang_major__) && (__clang_major__ == 7)
0769       2
0770 #else
0771       std::is_floating_point<T>::value && std::numeric_limits<T>::digits && (std::numeric_limits<T>::digits <= 54) ? 0 :
0772       std::is_floating_point<T>::value && std::numeric_limits<T>::digits && (std::numeric_limits<T>::digits <= 64) ? 1 : 2
0773 #endif
0774    > precision_tag_type;
0775    return policies::checked_narrowing_cast<result_type, Policy>(detail::ellint_k_imp(static_cast<value_type>(k), pol, precision_tag_type()), "boost::math::ellint_1<%1%>(%1%)");
0776 }
0777
0778 template <class T1, class T2>
0779 BOOST_MATH_FORCEINLINE typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi, const std::false_type&)
0780 {
0781    return boost::math::ellint_1(k, phi, policies::policy<>());
0782 }
0783
0784 }
0785
0786 // Complete elliptic integral (Legendre form) of the first kind
0787 template <typename T>
0788 BOOST_MATH_FORCEINLINE typename tools::promote_args<T>::type ellint_1(T k)
0789 {
0790    return ellint_1(k, policies::policy<>());
0791 }
0792
0793 // Elliptic integral (Legendre form) of the first kind
0794 template <class T1, class T2, class Policy>
0795 BOOST_MATH_FORCEINLINE typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi, const Policy& pol)  // LCOV_EXCL_LINE gcc misses this but sees the function body, strange!
0796 {
0797    typedef typename tools::promote_args<T1, T2>::type result_type;
0798    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0799    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%)");
0800 }
0801
0802 template <class T1, class T2>
0803 BOOST_MATH_FORCEINLINE typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi)
0804 {
0805    typedef typename policies::is_policy<T2>::type tag_type;
0806    return detail::ellint_1(k, phi, tag_type());
0807 }
0808
0809 }} // namespaces
0810
0811 #endif // BOOST_MATH_ELLINT_1_HPP
0812