Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-09-15 08:40:10 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!