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