757 lines · plain
1// Copyright (c) 2006 Xiaogang Zhang2// Copyright (c) 2006 John Maddock3// Copyright (c) 2024 Matt Borland4// Use, modification and distribution are subject to the5// Boost Software License, Version 1.0. (See accompanying file6// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)7//8// History:9// XZ wrote the original of this file as part of the Google10// Summer of Code 2006. JM modified it to fit into the11// Boost.Math conceptual framework better, and to ensure12// that the code continues to work no matter how many digits13// type T has.14 15#ifndef BOOST_MATH_ELLINT_2_HPP16#define BOOST_MATH_ELLINT_2_HPP17 18#ifdef _MSC_VER19#pragma once20#endif21 22#include <boost/math/tools/config.hpp>23#include <boost/math/tools/numeric_limits.hpp>24#include <boost/math/tools/type_traits.hpp>25#include <boost/math/special_functions/math_fwd.hpp>26#include <boost/math/special_functions/ellint_rf.hpp>27#include <boost/math/special_functions/ellint_rd.hpp>28#include <boost/math/special_functions/ellint_rg.hpp>29#include <boost/math/constants/constants.hpp>30#include <boost/math/policies/error_handling.hpp>31#include <boost/math/tools/workaround.hpp>32#include <boost/math/special_functions/round.hpp>33 34// Elliptic integrals (complete and incomplete) of the second kind35// Carlson, Numerische Mathematik, vol 33, 1 (1979)36 37namespace boost { namespace math {38 39template <class T1, class T2, class Policy>40BOOST_MATH_GPU_ENABLED typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi, const Policy& pol);41 42namespace detail{43 44template <typename T, typename Policy>45BOOST_MATH_GPU_ENABLED BOOST_MATH_FORCEINLINE T ellint_e_imp(T k, const Policy& pol, const boost::math::integral_constant<int, 0>&);46template <typename T, typename Policy>47BOOST_MATH_GPU_ENABLED BOOST_MATH_FORCEINLINE T ellint_e_imp(T k, const Policy& pol, const boost::math::integral_constant<int, 1>&);48template <typename T, typename Policy>49BOOST_MATH_GPU_ENABLED BOOST_MATH_FORCEINLINE T ellint_e_imp(T k, const Policy& pol, const boost::math::integral_constant<int, 2>&);50 51// Elliptic integral (Legendre form) of the second kind52template <typename T, typename Policy>53BOOST_MATH_GPU_ENABLED T ellint_e_imp(T phi, T k, const Policy& pol)54{55 BOOST_MATH_STD_USING56 using namespace boost::math::tools;57 using namespace boost::math::constants;58 59 bool invert = false;60 if (phi == 0)61 return 0;62 63 if(phi < 0)64 {65 phi = fabs(phi);66 invert = true;67 }68 69 T result;70 71 if(phi >= tools::max_value<T>())72 {73 // Need to handle infinity as a special case:74 result = policies::raise_overflow_error<T>("boost::math::ellint_e<%1%>(%1%,%1%)", nullptr, pol);75 }76 else if(phi > 1 / tools::epsilon<T>())77 {78 typedef boost::math::integral_constant<int,79 boost::math::is_floating_point<T>::value&& boost::math::numeric_limits<T>::digits && (boost::math::numeric_limits<T>::digits <= 54) ? 0 :80 boost::math::is_floating_point<T>::value && boost::math::numeric_limits<T>::digits && (boost::math::numeric_limits<T>::digits <= 64) ? 1 : 281 > precision_tag_type;82 // Phi is so large that phi%pi is necessarily zero (or garbage),83 // just return the second part of the duplication formula:84 result = 2 * phi * ellint_e_imp(k, pol, precision_tag_type()) / constants::pi<T>();85 }86 else if(k == 0)87 {88 return invert ? T(-phi) : phi;89 }90 else if(fabs(k) == 1)91 {92 //93 // For k = 1 ellipse actually turns to a line and every pi/2 in phi is exactly 1 in arc length94 // Periodicity though is in pi, curve follows sin(pi) for 0 <= phi <= pi/2 and then95 // 2 - sin(pi- phi) = 2 + sin(phi - pi) for pi/2 <= phi <= pi, so general form is:96 //97 // 2n + sin(phi - n * pi) ; |phi - n * pi| <= pi / 298 //99 T m = boost::math::round(phi / boost::math::constants::pi<T>());100 T remains = phi - m * boost::math::constants::pi<T>();101 T value = 2 * m + sin(remains);102 103 // negative arc length for negative phi104 return invert ? -value : value;105 }106 else107 {108 // Carlson's algorithm works only for |phi| <= pi/2,109 // use the integrand's periodicity to normalize phi110 //111 // Xiaogang's original code used a cast to long long here112 // but that fails if T has more digits than a long long,113 // so rewritten to use fmod instead:114 //115 T rphi = boost::math::tools::fmod_workaround(phi, T(constants::half_pi<T>()));116 T m = boost::math::round((phi - rphi) / constants::half_pi<T>());117 int s = 1;118 if(boost::math::tools::fmod_workaround(m, T(2)) > T(0.5))119 {120 m += 1;121 s = -1;122 rphi = constants::half_pi<T>() - rphi;123 }124 T k2 = k * k;125 if(boost::math::pow<3>(rphi) * k2 / 6 < tools::epsilon<T>() * fabs(rphi))126 {127 // See http://functions.wolfram.com/EllipticIntegrals/EllipticE2/06/01/03/0001/128 result = s * rphi;129 }130 else131 {132 // http://dlmf.nist.gov/19.25#E10133 T sinp = sin(rphi);134 if (k2 * sinp * sinp >= 1)135 {136 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);137 }138 T cosp = cos(rphi);139 T c = 1 / (sinp * sinp);140 T cm1 = cosp * cosp / (sinp * sinp); // c - 1141 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))));142 }143 if (m != 0)144 {145 typedef boost::math::integral_constant<int,146 boost::math::is_floating_point<T>::value&& boost::math::numeric_limits<T>::digits && (boost::math::numeric_limits<T>::digits <= 54) ? 0 :147 boost::math::is_floating_point<T>::value && boost::math::numeric_limits<T>::digits && (boost::math::numeric_limits<T>::digits <= 64) ? 1 : 2148 > precision_tag_type;149 result += m * ellint_e_imp(k, pol, precision_tag_type());150 }151 }152 return invert ? T(-result) : result;153}154 155// Complete elliptic integral (Legendre form) of the second kind156template <typename T, typename Policy>157BOOST_MATH_GPU_ENABLED T ellint_e_imp(T k, const Policy& pol, boost::math::integral_constant<int, 2> const&)158{159 BOOST_MATH_STD_USING160 using namespace boost::math::tools;161 162 if (abs(k) > 1)163 {164 return policies::raise_domain_error<T>("boost::math::ellint_e<%1%>(%1%)", "Got k = %1%, function requires |k| <= 1", k, pol);165 }166 if (abs(k) == 1)167 {168 return static_cast<T>(1);169 }170 171 T x = 0;172 T t = k * k;173 T y = 1 - t;174 T z = 1;175 T value = 2 * ellint_rg_imp(x, y, z, pol);176 177 return value;178}179//180// Special versions for double and 80-bit long double precision,181// double precision versions use the coefficients from:182// "Fast computation of complete elliptic integrals and Jacobian elliptic functions",183// Celestial Mechanics and Dynamical Astronomy, April 2012.184// 185// Higher precision coefficients for 80-bit long doubles can be calculated186// using for example:187// Table[N[SeriesCoefficient[ EllipticE [ m ], { m, 875/1000, i} ], 20], {i, 0, 24}]188// and checking the value of the first neglected term with:189// N[SeriesCoefficient[ EllipticE [ m ], { m, 875/1000, 24} ], 20] * (2.5/100)^24190// 191// For m > 0.9 we don't use the method of the paper above, but simply call our192// existing routines.193//194template <typename T, typename Policy>195BOOST_MATH_GPU_ENABLED BOOST_MATH_FORCEINLINE T ellint_e_imp(T k, const Policy& pol, boost::math::integral_constant<int, 0> const&)196{197 BOOST_MATH_STD_USING198 using namespace boost::math::tools;199 200 T m = k * k;201 switch (static_cast<int>(20 * m))202 {203 case 0:204 case 1:205 //if (m < 0.1)206 {207 constexpr T coef[] =208 {209 static_cast<T>(1.550973351780472328),210 -static_cast<T>(0.400301020103198524),211 -static_cast<T>(0.078498619442941939),212 -static_cast<T>(0.034318853117591992),213 -static_cast<T>(0.019718043317365499),214 -static_cast<T>(0.013059507731993309),215 -static_cast<T>(0.009442372874146547),216 -static_cast<T>(0.007246728512402157),217 -static_cast<T>(0.005807424012956090),218 -static_cast<T>(0.004809187786009338),219 -static_cast<T>(0.004086399233255150)220 };221 return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.05));222 }223 case 2:224 case 3:225 //else if (m < 0.2)226 {227 constexpr T coef[] =228 {229 static_cast<T>(1.510121832092819728),230 -static_cast<T>(0.417116333905867549),231 -static_cast<T>(0.090123820404774569),232 -static_cast<T>(0.043729944019084312),233 -static_cast<T>(0.027965493064761785),234 -static_cast<T>(0.020644781177568105),235 -static_cast<T>(0.016650786739707238),236 -static_cast<T>(0.014261960828842520),237 -static_cast<T>(0.012759847429264803),238 -static_cast<T>(0.011799303775587354),239 -static_cast<T>(0.011197445703074968)240 };241 return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.15));242 }243 case 4:244 case 5:245 //else if (m < 0.3)246 {247 constexpr T coef[] =248 {249 static_cast<T>(1.467462209339427155),250 -static_cast<T>(0.436576290946337775),251 -static_cast<T>(0.105155557666942554),252 -static_cast<T>(0.057371843593241730),253 -static_cast<T>(0.041391627727340220),254 -static_cast<T>(0.034527728505280841),255 -static_cast<T>(0.031495443512532783),256 -static_cast<T>(0.030527000890325277),257 -static_cast<T>(0.030916984019238900),258 -static_cast<T>(0.032371395314758122),259 -static_cast<T>(0.034789960386404158)260 };261 return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.25));262 }263 case 6:264 case 7:265 //else if (m < 0.4)266 {267 constexpr T coef[] =268 {269 static_cast<T>(1.422691133490879171),270 -static_cast<T>(0.459513519621048674),271 -static_cast<T>(0.125250539822061878),272 -static_cast<T>(0.078138545094409477),273 -static_cast<T>(0.064714278472050002),274 -static_cast<T>(0.062084339131730311),275 -static_cast<T>(0.065197032815572477),276 -static_cast<T>(0.072793895362578779),277 -static_cast<T>(0.084959075171781003),278 -static_cast<T>(0.102539850131045997),279 -static_cast<T>(0.127053585157696036),280 -static_cast<T>(0.160791120691274606)281 };282 return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.35));283 }284 case 8:285 case 9:286 //else if (m < 0.5)287 {288 constexpr T coef[] =289 {290 static_cast<T>(1.375401971871116291),291 -static_cast<T>(0.487202183273184837),292 -static_cast<T>(0.153311701348540228),293 -static_cast<T>(0.111849444917027833),294 -static_cast<T>(0.108840952523135768),295 -static_cast<T>(0.122954223120269076),296 -static_cast<T>(0.152217163962035047),297 -static_cast<T>(0.200495323642697339),298 -static_cast<T>(0.276174333067751758),299 -static_cast<T>(0.393513114304375851),300 -static_cast<T>(0.575754406027879147),301 -static_cast<T>(0.860523235727239756),302 -static_cast<T>(1.308833205758540162)303 };304 return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.45));305 }306 case 10:307 case 11:308 //else if (m < 0.6)309 {310 constexpr T coef[] =311 {312 static_cast<T>(1.325024497958230082),313 -static_cast<T>(0.521727647557566767),314 -static_cast<T>(0.194906430482126213),315 -static_cast<T>(0.171623726822011264),316 -static_cast<T>(0.202754652926419141),317 -static_cast<T>(0.278798953118534762),318 -static_cast<T>(0.420698457281005762),319 -static_cast<T>(0.675948400853106021),320 -static_cast<T>(1.136343121839229244),321 -static_cast<T>(1.976721143954398261),322 -static_cast<T>(3.531696773095722506),323 -static_cast<T>(6.446753640156048150),324 -static_cast<T>(11.97703130208884026)325 };326 return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.55));327 }328 case 12:329 case 13:330 //else if (m < 0.7)331 {332 constexpr T coef[] =333 {334 static_cast<T>(1.270707479650149744),335 -static_cast<T>(0.566839168287866583),336 -static_cast<T>(0.262160793432492598),337 -static_cast<T>(0.292244173533077419),338 -static_cast<T>(0.440397840850423189),339 -static_cast<T>(0.774947641381397458),340 -static_cast<T>(1.498870837987561088),341 -static_cast<T>(3.089708310445186667),342 -static_cast<T>(6.667595903381001064),343 -static_cast<T>(14.89436036517319078),344 -static_cast<T>(34.18120574251449024),345 -static_cast<T>(80.15895841905397306),346 -static_cast<T>(191.3489480762984920),347 -static_cast<T>(463.5938853480342030),348 -static_cast<T>(1137.380822169360061)349 };350 return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.65));351 }352 case 14:353 case 15:354 //else if (m < 0.8)355 {356 constexpr T coef[] =357 {358 static_cast<T>(1.211056027568459525),359 -static_cast<T>(0.630306413287455807),360 -static_cast<T>(0.387166409520669145),361 -static_cast<T>(0.592278235311934603),362 -static_cast<T>(1.237555584513049844),363 -static_cast<T>(3.032056661745247199),364 -static_cast<T>(8.181688221573590762),365 -static_cast<T>(23.55507217389693250),366 -static_cast<T>(71.04099935893064956),367 -static_cast<T>(221.8796853192349888),368 -static_cast<T>(712.1364793277635425),369 -static_cast<T>(2336.125331440396407),370 -static_cast<T>(7801.945954775964673),371 -static_cast<T>(26448.19586059191933),372 -static_cast<T>(90799.48341621365251),373 -static_cast<T>(315126.0406449163424),374 -static_cast<T>(1104011.344311591159)375 };376 return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.75));377 }378 case 16:379 //else if (m < 0.85)380 {381 constexpr T coef[] =382 {383 static_cast<T>(1.161307152196282836),384 -static_cast<T>(0.701100284555289548),385 -static_cast<T>(0.580551474465437362),386 -static_cast<T>(1.243693061077786614),387 -static_cast<T>(3.679383613496634879),388 -static_cast<T>(12.81590924337895775),389 -static_cast<T>(49.25672530759985272),390 -static_cast<T>(202.1818735434090269),391 -static_cast<T>(869.8602699308701437),392 -static_cast<T>(3877.005847313289571),393 -static_cast<T>(17761.70710170939814),394 -static_cast<T>(83182.69029154232061),395 -static_cast<T>(396650.4505013548170),396 -static_cast<T>(1920033.413682634405)397 };398 return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.825));399 }400 case 17:401 //else if (m < 0.90)402 {403 constexpr T coef[] =404 {405 static_cast<T>(1.124617325119752213),406 -static_cast<T>(0.770845056360909542),407 -static_cast<T>(0.844794053644911362),408 -static_cast<T>(2.490097309450394453),409 -static_cast<T>(10.23971741154384360),410 -static_cast<T>(49.74900546551479866),411 -static_cast<T>(267.0986675195705196),412 -static_cast<T>(1532.665883825229947),413 -static_cast<T>(9222.313478526091951),414 -static_cast<T>(57502.51612140314030),415 -static_cast<T>(368596.1167416106063),416 -static_cast<T>(2415611.088701091428),417 -static_cast<T>(16120097.81581656797),418 -static_cast<T>(109209938.5203089915),419 -static_cast<T>(749380758.1942496220),420 -static_cast<T>(5198725846.725541393),421 -static_cast<T>(36409256888.12139973)422 };423 return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.875));424 }425 default:426 //427 // All cases where m > 0.9428 // including all error handling:429 //430 return ellint_e_imp(k, pol, boost::math::integral_constant<int, 2>());431 }432}433template <typename T, typename Policy>434BOOST_MATH_GPU_ENABLED BOOST_MATH_FORCEINLINE T ellint_e_imp(T k, const Policy& pol, boost::math::integral_constant<int, 1> const&)435{436 BOOST_MATH_STD_USING437 using namespace boost::math::tools;438 439 T m = k * k;440 switch (static_cast<int>(20 * m))441 {442 case 0:443 case 1:444 //if (m < 0.1)445 {446 constexpr T coef[] =447 {448 1.5509733517804723277L,449 -0.40030102010319852390L,450 -0.078498619442941939212L,451 -0.034318853117591992417L,452 -0.019718043317365499309L,453 -0.013059507731993309191L,454 -0.0094423728741465473894L,455 -0.0072467285124021568126L,456 -0.0058074240129560897940L,457 -0.0048091877860093381762L,458 -0.0040863992332551506768L,459 -0.0035450302604139562644L,460 -0.0031283511188028336315L461 };462 return boost::math::tools::evaluate_polynomial(coef, m - 0.05L);463 }464 case 2:465 case 3:466 //else if (m < 0.2)467 {468 constexpr T coef[] =469 {470 1.5101218320928197276L,471 -0.41711633390586754922L,472 -0.090123820404774568894L,473 -0.043729944019084311555L,474 -0.027965493064761784548L,475 -0.020644781177568105268L,476 -0.016650786739707238037L,477 -0.014261960828842519634L,478 -0.012759847429264802627L,479 -0.011799303775587354169L,480 -0.011197445703074968018L,481 -0.010850368064799902735L,482 -0.010696133481060989818L483 };484 return boost::math::tools::evaluate_polynomial(coef, m - 0.15L);485 }486 case 4:487 case 5:488 //else if (m < 0.3L)489 {490 constexpr T coef[] =491 {492 1.4674622093394271555L,493 -0.43657629094633777482L,494 -0.10515555766694255399L,495 -0.057371843593241729895L,496 -0.041391627727340220236L,497 -0.034527728505280841188L,498 -0.031495443512532782647L,499 -0.030527000890325277179L,500 -0.030916984019238900349L,501 -0.032371395314758122268L,502 -0.034789960386404158240L,503 -0.038182654612387881967L,504 -0.042636187648900252525L,505 -0.048302272505241634467506 };507 return boost::math::tools::evaluate_polynomial(coef, m - 0.25L);508 }509 case 6:510 case 7:511 //else if (m < 0.4L)512 {513 constexpr T coef[] =514 {515 1.4226911334908791711L,516 -0.45951351962104867394L,517 -0.12525053982206187849L,518 -0.078138545094409477156L,519 -0.064714278472050001838L,520 -0.062084339131730310707L,521 -0.065197032815572476910L,522 -0.072793895362578779473L,523 -0.084959075171781003264L,524 -0.10253985013104599679L,525 -0.12705358515769603644L,526 -0.16079112069127460621L,527 -0.20705400012405941376L,528 -0.27053164884730888948L529 };530 return boost::math::tools::evaluate_polynomial(coef, m - 0.35L);531 }532 case 8:533 case 9:534 //else if (m < 0.5L)535 {536 constexpr T coef[] =537 {538 1.3754019718711162908L,539 -0.48720218327318483652L,540 -0.15331170134854022753L,541 -0.11184944491702783273L,542 -0.10884095252313576755L,543 -0.12295422312026907610L,544 -0.15221716396203504746L,545 -0.20049532364269733857L,546 -0.27617433306775175837L,547 -0.39351311430437585139L,548 -0.57575440602787914711L,549 -0.86052323572723975634L,550 -1.3088332057585401616L,551 -2.0200280559452241745L,552 -3.1566019548237606451L553 };554 return boost::math::tools::evaluate_polynomial(coef, m - 0.45L);555 }556 case 10:557 case 11:558 //else if (m < 0.6L)559 {560 constexpr T coef[] =561 {562 1.3250244979582300818L,563 -0.52172764755756676713L,564 -0.19490643048212621262L,565 -0.17162372682201126365L,566 -0.20275465292641914128L,567 -0.27879895311853476205L,568 -0.42069845728100576224L,569 -0.67594840085310602110L,570 -1.1363431218392292440L,571 -1.9767211439543982613L,572 -3.5316967730957225064L,573 -6.4467536401560481499L,574 -11.977031302088840261L,575 -22.581360948073964469L,576 -43.109479829481450573L,577 -83.186290908288807424L578 };579 return boost::math::tools::evaluate_polynomial(coef, m - 0.55L);580 }581 case 12:582 case 13:583 //else if (m < 0.7L)584 {585 constexpr T coef[] =586 {587 1.2707074796501497440L,588 -0.56683916828786658286L,589 -0.26216079343249259779L,590 -0.29224417353307741931L,591 -0.44039784085042318909L,592 -0.77494764138139745824L,593 -1.4988708379875610880L,594 -3.0897083104451866665L,595 -6.6675959033810010645L,596 -14.894360365173190775L,597 -34.181205742514490240L,598 -80.158958419053973056L,599 -191.34894807629849204L,600 -463.59388534803420301L,601 -1137.3808221693600606L,602 -2820.7073786352269339L,603 -7061.1382244658715621L,604 -17821.809331816437058L,605 -45307.849987201897801L606 };607 return boost::math::tools::evaluate_polynomial(coef, m - 0.65L);608 }609 case 14:610 case 15:611 //else if (m < 0.8L)612 {613 constexpr T coef[] =614 {615 1.2110560275684595248L,616 -0.63030641328745580709L,617 -0.38716640952066914514L,618 -0.59227823531193460257L,619 -1.2375555845130498445L,620 -3.0320566617452471986L,621 -8.1816882215735907624L,622 -23.555072173896932503L,623 -71.040999358930649565L,624 -221.87968531923498875L,625 -712.13647932776354253L,626 -2336.1253314403964072L,627 -7801.9459547759646726L,628 -26448.195860591919335L,629 -90799.483416213652512L,630 -315126.04064491634241L,631 -1.1040113443115911589e6L,632 -3.8998018348056769095e6L,633 -1.3876249116223745041e7L,634 -4.9694982823537861149e7L,635 -1.7900668836197342979e8L,636 -6.4817399873722371964e8L637 };638 return boost::math::tools::evaluate_polynomial(coef, m - 0.75L);639 }640 case 16:641 //else if (m < 0.85L)642 {643 constexpr T coef[] =644 {645 1.1613071521962828360L,646 -0.70110028455528954752L,647 -0.58055147446543736163L,648 -1.2436930610777866138L,649 -3.6793836134966348789L,650 -12.815909243378957753L,651 -49.256725307599852720L,652 -202.18187354340902693L,653 -869.86026993087014372L,654 -3877.0058473132895713L,655 -17761.707101709398174L,656 -83182.690291542320614L,657 -396650.45050135481698L,658 -1.9200334136826344054e6L,659 -9.4131321779500838352e6L,660 -4.6654858837335370627e7L,661 -2.3343549352617609390e8L,662 -1.1776928223045913454e9L,663 -5.9850851892915740401e9L,664 -3.0614702984618644983e10L665 };666 return boost::math::tools::evaluate_polynomial(coef, m - 0.825L);667 }668 case 17:669 //else if (m < 0.90L)670 {671 constexpr T coef[] =672 {673 1.1246173251197522132L,674 -0.77084505636090954218L,675 -0.84479405364491136236L,676 -2.4900973094503944527L,677 -10.239717411543843601L,678 -49.749005465514798660L,679 -267.09866751957051961L,680 -1532.6658838252299468L,681 -9222.3134785260919507L,682 -57502.516121403140303L,683 -368596.11674161060626L,684 -2.4156110887010914281e6L,685 -1.6120097815816567971e7L,686 -1.0920993852030899148e8L,687 -7.4938075819424962198e8L,688 -5.1987258467255413931e9L,689 -3.6409256888121399726e10L,690 -2.5711802891217393544e11L,691 -1.8290904062978796996e12L,692 -1.3096838781743248404e13L,693 -9.4325465851415135118e13L,694 -6.8291980829471896669e14L695 };696 return boost::math::tools::evaluate_polynomial(coef, m - 0.875L);697 }698 default:699 //700 // All cases where m > 0.9701 // including all error handling:702 //703 return ellint_e_imp(k, pol, boost::math::integral_constant<int, 2>());704 }705}706 707template <typename T, typename Policy>708BOOST_MATH_GPU_ENABLED typename tools::promote_args<T>::type ellint_2(T k, const Policy& pol, const boost::math::true_type&)709{710 typedef typename tools::promote_args<T>::type result_type;711 typedef typename policies::evaluation<result_type, Policy>::type value_type;712 typedef boost::math::integral_constant<int,713 boost::math::is_floating_point<T>::value&& boost::math::numeric_limits<T>::digits && (boost::math::numeric_limits<T>::digits <= 54) ? 0 :714 boost::math::is_floating_point<T>::value && boost::math::numeric_limits<T>::digits && (boost::math::numeric_limits<T>::digits <= 64) ? 1 : 2715 > precision_tag_type;716 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%)");717}718 719// Elliptic integral (Legendre form) of the second kind720template <class T1, class T2>721BOOST_MATH_GPU_ENABLED typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi, const boost::math::false_type&)722{723 return boost::math::ellint_2(k, phi, policies::policy<>());724}725 726} // detail727 728// Elliptic integral (Legendre form) of the second kind729template <class T1, class T2>730BOOST_MATH_GPU_ENABLED typename tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi)731{732 typedef typename policies::is_policy<T2>::type tag_type;733 return detail::ellint_2(k, phi, tag_type());734}735 736template <class T1, class T2, class Policy>737BOOST_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!738{739 typedef typename tools::promote_args<T1, T2>::type result_type;740 typedef typename policies::evaluation<result_type, Policy>::type value_type;741 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%)");742}743 744 745// Complete elliptic integral (Legendre form) of the second kind746template <typename T>747BOOST_MATH_GPU_ENABLED typename tools::promote_args<T>::type ellint_2(T k)748{749 return ellint_2(k, policies::policy<>());750}751 752 753}} // namespaces754 755#endif // BOOST_MATH_ELLINT_2_HPP756 757