math/special_functions/ellint_1.hpp

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