816 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_1_HPP16#define BOOST_MATH_ELLINT_1_HPP17 18#ifdef _MSC_VER19#pragma once20#endif21 22#include <boost/math/tools/config.hpp>23#include <boost/math/tools/type_traits.hpp>24#include <boost/math/special_functions/math_fwd.hpp>25#include <boost/math/special_functions/ellint_rf.hpp>26#include <boost/math/constants/constants.hpp>27#include <boost/math/policies/error_handling.hpp>28#include <boost/math/tools/workaround.hpp>29#include <boost/math/special_functions/round.hpp>30 31// Elliptic integrals (complete and incomplete) of the first kind32// Carlson, Numerische Mathematik, vol 33, 1 (1979)33 34namespace boost { namespace math {35 36template <class T1, class T2, class Policy>37BOOST_MATH_GPU_ENABLED typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi, const Policy& pol);38 39namespace detail{40 41template <typename T, typename Policy>42BOOST_MATH_GPU_ENABLED BOOST_MATH_FORCEINLINE T ellint_k_imp(T k, const Policy& pol, boost::math::integral_constant<int, 0> const&);43template <typename T, typename Policy>44BOOST_MATH_GPU_ENABLED BOOST_MATH_FORCEINLINE T ellint_k_imp(T k, const Policy& pol, boost::math::integral_constant<int, 1> const&);45template <typename T, typename Policy>46BOOST_MATH_GPU_ENABLED BOOST_MATH_FORCEINLINE T ellint_k_imp(T k, const Policy& pol, boost::math::integral_constant<int, 2> const&);47template <typename T, typename Policy>48BOOST_MATH_GPU_ENABLED BOOST_MATH_FORCEINLINE T ellint_k_imp(T k, const Policy& pol, T one_minus_k2);49 50// Elliptic integral (Legendre form) of the first kind51template <typename T, typename Policy>52BOOST_MATH_GPU_ENABLED T ellint_f_imp(T phi, T k, const Policy& pol, T one_minus_k2)53{54 BOOST_MATH_STD_USING55 using namespace boost::math::tools;56 using namespace boost::math::constants;57 58 constexpr auto function = "boost::math::ellint_f<%1%>(%1%,%1%)";59 BOOST_MATH_INSTRUMENT_VARIABLE(phi);60 BOOST_MATH_INSTRUMENT_VARIABLE(k);61 BOOST_MATH_INSTRUMENT_VARIABLE(function);62 63 bool invert = false;64 if(phi < 0)65 {66 BOOST_MATH_INSTRUMENT_VARIABLE(phi);67 phi = fabs(phi);68 invert = true;69 }70 71 T result;72 73 if(phi >= tools::max_value<T>())74 {75 // Need to handle infinity as a special case:76 result = policies::raise_overflow_error<T>(function, nullptr, pol);77 BOOST_MATH_INSTRUMENT_VARIABLE(result);78 }79 else if(phi > 1 / tools::epsilon<T>())80 {81 // Phi is so large that phi%pi is necessarily zero (or garbage),82 // just return the second part of the duplication formula:83 result = 2 * phi * ellint_k_imp(k, pol, one_minus_k2) / constants::pi<T>();84 BOOST_MATH_INSTRUMENT_VARIABLE(result);85 }86 else87 {88 // Carlson's algorithm works only for |phi| <= pi/2,89 // use the integrand's periodicity to normalize phi90 //91 // Xiaogang's original code used a cast to long long here92 // but that fails if T has more digits than a long long,93 // so rewritten to use fmod instead:94 //95 BOOST_MATH_INSTRUMENT_CODE("pi/2 = " << constants::pi<T>() / 2);96 T rphi = boost::math::tools::fmod_workaround(phi, T(constants::half_pi<T>()));97 BOOST_MATH_INSTRUMENT_VARIABLE(rphi);98 T m = boost::math::round((phi - rphi) / constants::half_pi<T>());99 BOOST_MATH_INSTRUMENT_VARIABLE(m);100 int s = 1;101 if(boost::math::tools::fmod_workaround(m, T(2)) > T(0.5))102 {103 m += 1;104 s = -1;105 rphi = constants::half_pi<T>() - rphi;106 BOOST_MATH_INSTRUMENT_VARIABLE(rphi);107 }108 T sinp = sin(rphi);109 sinp *= sinp;110 if (sinp * k * k >= 1)111 {112 return policies::raise_domain_error<T>(function,113 "Got k^2 * sin^2(phi) = %1%, but the function requires this < 1", sinp * k * k, pol);114 }115 T cosp = cos(rphi);116 cosp *= cosp;117 BOOST_MATH_INSTRUMENT_VARIABLE(sinp);118 BOOST_MATH_INSTRUMENT_VARIABLE(cosp);119 if(sinp > tools::min_value<T>())120 {121 BOOST_MATH_ASSERT(rphi != 0); // precondition, can't be true if sin(rphi) != 0.122 //123 // Use http://dlmf.nist.gov/19.25#E5, note that124 // c-1 simplifies to cot^2(rphi) which avoids cancellation.125 // Likewise c - k^2 is the same as (c - 1) + (1 - k^2).126 //127 T c = 1 / sinp;128 T c_minus_one = cosp / sinp;129 T arg2;130 if (k != 0)131 {132 T cross = fabs(c / (k * k));133 if ((cross > 0.9f) && (cross < 1.1f))134 arg2 = c_minus_one + one_minus_k2;135 else136 arg2 = c - k * k;137 }138 else139 arg2 = c;140 result = static_cast<T>(s * ellint_rf_imp(c_minus_one, arg2, c, pol));141 }142 else143 result = s * sin(rphi);144 BOOST_MATH_INSTRUMENT_VARIABLE(result);145 if(m != 0)146 {147 result += m * ellint_k_imp(k, pol, one_minus_k2);148 BOOST_MATH_INSTRUMENT_VARIABLE(result);149 }150 }151 return invert ? T(-result) : result;152}153 154template <typename T, typename Policy>155BOOST_MATH_GPU_ENABLED inline T ellint_f_imp(T phi, T k, const Policy& pol)156{157 return ellint_f_imp(phi, k, pol, T(1 - k * k));158}159 160// Complete elliptic integral (Legendre form) of the first kind161template <typename T, typename Policy>162BOOST_MATH_GPU_ENABLED T ellint_k_imp(T k, const Policy& pol, T one_minus_k2)163{164 BOOST_MATH_STD_USING165 using namespace boost::math::tools;166 167 constexpr auto function = "boost::math::ellint_k<%1%>(%1%)";168 169 if (abs(k) > 1)170 {171 return policies::raise_domain_error<T>(function, "Got k = %1%, function requires |k| <= 1", k, pol);172 }173 if (abs(k) == 1)174 {175 return policies::raise_overflow_error<T>(function, nullptr, pol);176 }177 178 T x = 0;179 T z = 1;180 T value = ellint_rf_imp(x, one_minus_k2, z, pol);181 182 return value;183}184template <typename T, typename Policy>185BOOST_MATH_GPU_ENABLED inline T ellint_k_imp(T k, const Policy& pol, boost::math::integral_constant<int, 2> const&)186{187 return ellint_k_imp(k, pol, T(1 - k * k));188}189 190//191// Special versions for double and 80-bit long double precision,192// double precision versions use the coefficients from:193// "Fast computation of complete elliptic integrals and Jacobian elliptic functions",194// Celestial Mechanics and Dynamical Astronomy, April 2012.195// 196// Higher precision coefficients for 80-bit long doubles can be calculated197// using for example:198// Table[N[SeriesCoefficient[ EllipticK [ m ], { m, 875/1000, i} ], 20], {i, 0, 24}]199// and checking the value of the first neglected term with:200// N[SeriesCoefficient[ EllipticK [ m ], { m, 875/1000, 24} ], 20] * (2.5/100)^24201// 202// For m > 0.9 we don't use the method of the paper above, but simply call our203// existing routines. The routine used in the above paper was tried (and is204// archived in the code below), but was found to have slightly higher error rates.205//206template <typename T, typename Policy>207BOOST_MATH_GPU_ENABLED BOOST_MATH_FORCEINLINE T ellint_k_imp(T k, const Policy& pol, boost::math::integral_constant<int, 0> const&)208{209 BOOST_MATH_STD_USING210 using namespace boost::math::tools;211 212 T m = k * k;213 214 switch (static_cast<int>(m * 20))215 {216 case 0:217 case 1:218 //if (m < 0.1)219 {220 constexpr T coef[] =221 {222 static_cast<T>(1.591003453790792180),223 static_cast<T>(0.416000743991786912),224 static_cast<T>(0.245791514264103415),225 static_cast<T>(0.179481482914906162),226 static_cast<T>(0.144556057087555150),227 static_cast<T>(0.123200993312427711),228 static_cast<T>(0.108938811574293531),229 static_cast<T>(0.098853409871592910),230 static_cast<T>(0.091439629201749751),231 static_cast<T>(0.085842591595413900),232 static_cast<T>(0.081541118718303215),233 static_cast<T>(0.078199656811256481910)234 };235 return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.05));236 }237 case 2:238 case 3:239 //else if (m < 0.2)240 {241 constexpr T coef[] =242 {243 static_cast<T>(1.635256732264579992),244 static_cast<T>(0.471190626148732291),245 static_cast<T>(0.309728410831499587),246 static_cast<T>(0.252208311773135699),247 static_cast<T>(0.226725623219684650),248 static_cast<T>(0.215774446729585976),249 static_cast<T>(0.213108771877348910),250 static_cast<T>(0.216029124605188282),251 static_cast<T>(0.223255831633057896),252 static_cast<T>(0.234180501294209925),253 static_cast<T>(0.248557682972264071),254 static_cast<T>(0.266363809892617521)255 };256 return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.15));257 }258 case 4:259 case 5:260 //else if (m < 0.3)261 {262 constexpr T coef[] =263 {264 static_cast<T>(1.685750354812596043),265 static_cast<T>(0.541731848613280329),266 static_cast<T>(0.401524438390690257),267 static_cast<T>(0.369642473420889090),268 static_cast<T>(0.376060715354583645),269 static_cast<T>(0.405235887085125919),270 static_cast<T>(0.453294381753999079),271 static_cast<T>(0.520518947651184205),272 static_cast<T>(0.609426039204995055),273 static_cast<T>(0.724263522282908870),274 static_cast<T>(0.871013847709812357),275 static_cast<T>(1.057652872753547036)276 };277 return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.25));278 }279 case 6:280 case 7:281 //else if (m < 0.4)282 {283 constexpr T coef[] =284 {285 static_cast<T>(1.744350597225613243),286 static_cast<T>(0.634864275371935304),287 static_cast<T>(0.539842564164445538),288 static_cast<T>(0.571892705193787391),289 static_cast<T>(0.670295136265406100),290 static_cast<T>(0.832586590010977199),291 static_cast<T>(1.073857448247933265),292 static_cast<T>(1.422091460675497751),293 static_cast<T>(1.920387183402304829),294 static_cast<T>(2.632552548331654201),295 static_cast<T>(3.652109747319039160),296 static_cast<T>(5.115867135558865806),297 static_cast<T>(7.224080007363877411)298 };299 return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.35));300 }301 case 8:302 case 9:303 //else if (m < 0.5)304 {305 constexpr T coef[] =306 {307 static_cast<T>(1.813883936816982644),308 static_cast<T>(0.763163245700557246),309 static_cast<T>(0.761928605321595831),310 static_cast<T>(0.951074653668427927),311 static_cast<T>(1.315180671703161215),312 static_cast<T>(1.928560693477410941),313 static_cast<T>(2.937509342531378755),314 static_cast<T>(4.594894405442878062),315 static_cast<T>(7.330071221881720772),316 static_cast<T>(11.87151259742530180),317 static_cast<T>(19.45851374822937738),318 static_cast<T>(32.20638657246426863),319 static_cast<T>(53.73749198700554656),320 static_cast<T>(90.27388602940998849)321 };322 return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.45));323 }324 case 10:325 case 11:326 //else if (m < 0.6)327 {328 constexpr T coef[] =329 {330 static_cast<T>(1.898924910271553526),331 static_cast<T>(0.950521794618244435),332 static_cast<T>(1.151077589959015808),333 static_cast<T>(1.750239106986300540),334 static_cast<T>(2.952676812636875180),335 static_cast<T>(5.285800396121450889),336 static_cast<T>(9.832485716659979747),337 static_cast<T>(18.78714868327559562),338 static_cast<T>(36.61468615273698145),339 static_cast<T>(72.45292395127771801),340 static_cast<T>(145.1079577347069102),341 static_cast<T>(293.4786396308497026),342 static_cast<T>(598.3851815055010179),343 static_cast<T>(1228.420013075863451),344 static_cast<T>(2536.529755382764488)345 };346 return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.55));347 }348 case 12:349 case 13:350 //else if (m < 0.7)351 {352 constexpr T coef[] =353 {354 static_cast<T>(2.007598398424376302),355 static_cast<T>(1.248457231212347337),356 static_cast<T>(1.926234657076479729),357 static_cast<T>(3.751289640087587680),358 static_cast<T>(8.119944554932045802),359 static_cast<T>(18.66572130873555361),360 static_cast<T>(44.60392484291437063),361 static_cast<T>(109.5092054309498377),362 static_cast<T>(274.2779548232413480),363 static_cast<T>(697.5598008606326163),364 static_cast<T>(1795.716014500247129),365 static_cast<T>(4668.381716790389910),366 static_cast<T>(12235.76246813664335),367 static_cast<T>(32290.17809718320818),368 static_cast<T>(85713.07608195964685),369 static_cast<T>(228672.1890493117096),370 static_cast<T>(612757.2711915852774)371 };372 return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.65));373 }374 case 14:375 case 15:376 //else if (m < static_cast<T>(0.8))377 {378 constexpr T coef[] =379 {380 static_cast<T>(2.156515647499643235),381 static_cast<T>(1.791805641849463243),382 static_cast<T>(3.826751287465713147),383 static_cast<T>(10.38672468363797208),384 static_cast<T>(31.40331405468070290),385 static_cast<T>(100.9237039498695416),386 static_cast<T>(337.3268282632272897),387 static_cast<T>(1158.707930567827917),388 static_cast<T>(4060.990742193632092),389 static_cast<T>(14454.00184034344795),390 static_cast<T>(52076.66107599404803),391 static_cast<T>(189493.6591462156887),392 static_cast<T>(695184.5762413896145),393 static_cast<T>(2567994.048255284686),394 static_cast<T>(9541921.966748386322),395 static_cast<T>(35634927.44218076174),396 static_cast<T>(133669298.4612040871),397 static_cast<T>(503352186.6866284541),398 static_cast<T>(1901975729.538660119),399 static_cast<T>(7208915015.330103756)400 };401 return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.75));402 }403 case 16:404 //else if (m < static_cast<T>(0.85))405 {406 constexpr T coef[] =407 {408 static_cast<T>(2.318122621712510589),409 static_cast<T>(2.616920150291232841),410 static_cast<T>(7.897935075731355823),411 static_cast<T>(30.50239715446672327),412 static_cast<T>(131.4869365523528456),413 static_cast<T>(602.9847637356491617),414 static_cast<T>(2877.024617809972641),415 static_cast<T>(14110.51991915180325),416 static_cast<T>(70621.44088156540229),417 static_cast<T>(358977.2665825309926),418 static_cast<T>(1847238.263723971684),419 static_cast<T>(9600515.416049214109),420 static_cast<T>(50307677.08502366879),421 static_cast<T>(265444188.6527127967),422 static_cast<T>(1408862325.028702687),423 static_cast<T>(7515687935.373774627)424 };425 return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.825));426 }427 case 17:428 //else if (m < static_cast<T>(0.90))429 {430 constexpr T coef[] =431 {432 static_cast<T>(2.473596173751343912),433 static_cast<T>(3.727624244118099310),434 static_cast<T>(15.60739303554930496),435 static_cast<T>(84.12850842805887747),436 static_cast<T>(506.9818197040613935),437 static_cast<T>(3252.277058145123644),438 static_cast<T>(21713.24241957434256),439 static_cast<T>(149037.0451890932766),440 static_cast<T>(1043999.331089990839),441 static_cast<T>(7427974.817042038995),442 static_cast<T>(53503839.67558661151),443 static_cast<T>(389249886.9948708474),444 static_cast<T>(2855288351.100810619),445 static_cast<T>(21090077038.76684053),446 static_cast<T>(156699833947.7902014),447 static_cast<T>(1170222242422.439893),448 static_cast<T>(8777948323668.937971),449 static_cast<T>(66101242752484.95041),450 static_cast<T>(499488053713388.7989),451 static_cast<T>(37859743397240299.20)452 };453 return boost::math::tools::evaluate_polynomial(coef, m - static_cast<T>(0.875));454 }455 default:456 //457 // This handles all cases where m > 0.9, 458 // including all error handling:459 //460 return ellint_k_imp(k, pol, boost::math::integral_constant<int, 2>());461#if 0462 else463 {464 T lambda_prime = (1 - sqrt(k)) / (2 * (1 + sqrt(k)));465 T k_prime = ellint_k(sqrt((1 - k) * (1 + k))); // K(m')466 T lambda_prime_4th = boost::math::pow<4>(lambda_prime);467 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;468 /*T q_prime_2 = lambda_prime469 + 2 * boost::math::pow<5>(lambda_prime)470 + 15 * boost::math::pow<9>(lambda_prime)471 + 150 * boost::math::pow<13>(lambda_prime)472 + 1707 * boost::math::pow<17>(lambda_prime)473 + 20910 * boost::math::pow<21>(lambda_prime);*/474 return -log(q_prime) * k_prime / boost::math::constants::pi<T>();475 }476#endif477 }478}479template <typename T, typename Policy>480BOOST_MATH_GPU_ENABLED BOOST_MATH_FORCEINLINE T ellint_k_imp(T k, const Policy& pol, boost::math::integral_constant<int, 1> const&)481{482 BOOST_MATH_STD_USING483 using namespace boost::math::tools;484 485 T m = k * k;486 switch (static_cast<int>(m * 20))487 {488 case 0:489 case 1:490 {491 constexpr T coef[] =492 {493 1.5910034537907921801L,494 0.41600074399178691174L,495 0.24579151426410341536L,496 0.17948148291490616181L,497 0.14455605708755514976L,498 0.12320099331242771115L,499 0.10893881157429353105L,500 0.098853409871592910399L,501 0.091439629201749751268L,502 0.085842591595413899672L,503 0.081541118718303214749L,504 0.078199656811256481910L,505 0.075592617535422415648L,506 0.073562939365441925050L507 };508 return boost::math::tools::evaluate_polynomial(coef, m - 0.05L);509 }510 case 2:511 case 3:512 {513 constexpr T coef[] =514 {515 1.6352567322645799924L,516 0.47119062614873229055L,517 0.30972841083149958708L,518 0.25220831177313569923L,519 0.22672562321968464974L,520 0.21577444672958597588L,521 0.21310877187734890963L,522 0.21602912460518828154L,523 0.22325583163305789567L,524 0.23418050129420992492L,525 0.24855768297226407136L,526 0.26636380989261752077L,527 0.28772845215611466775L,528 0.31290024539780334906L,529 0.34223105446381299902L530 };531 return boost::math::tools::evaluate_polynomial(coef, m - 0.15L);532 }533 case 4:534 case 5:535 {536 constexpr T coef[] =537 {538 1.6857503548125960429L,539 0.54173184861328032882L,540 0.40152443839069025682L,541 0.36964247342088908995L,542 0.37606071535458364462L,543 0.40523588708512591863L,544 0.45329438175399907924L,545 0.52051894765118420473L,546 0.60942603920499505544L,547 0.72426352228290886975L,548 0.87101384770981235737L,549 1.0576528727535470365L,550 1.2945970872087764321L,551 1.5953368253888783747L,552 1.9772844873556364793L,553 2.4628890581910021287L554 };555 return boost::math::tools::evaluate_polynomial(coef, m - 0.25L);556 }557 case 6:558 case 7:559 {560 constexpr T coef[] =561 {562 1.7443505972256132429L,563 0.63486427537193530383L,564 0.53984256416444553751L,565 0.57189270519378739093L,566 0.67029513626540610034L,567 0.83258659001097719939L,568 1.0738574482479332654L,569 1.4220914606754977514L,570 1.9203871834023048288L,571 2.6325525483316542006L,572 3.6521097473190391602L,573 5.1158671355588658061L,574 7.2240800073638774108L,575 10.270306349944787227L,576 14.685616935355757348L,577 21.104114212004582734L,578 30.460132808575799413L,579 };580 return boost::math::tools::evaluate_polynomial(coef, m - 0.35L);581 }582 case 8:583 case 9:584 {585 constexpr T coef[] =586 {587 1.8138839368169826437L,588 0.76316324570055724607L,589 0.76192860532159583095L,590 0.95107465366842792679L,591 1.3151806717031612153L,592 1.9285606934774109412L,593 2.9375093425313787550L,594 4.5948944054428780618L,595 7.3300712218817207718L,596 11.871512597425301798L,597 19.458513748229377383L,598 32.206386572464268628L,599 53.737491987005546559L,600 90.273886029409988491L,601 152.53312130253275268L,602 259.02388747148299086L,603 441.78537518096201946L,604 756.39903981567380952L605 };606 return boost::math::tools::evaluate_polynomial(coef, m - 0.45L);607 }608 case 10:609 case 11:610 {611 constexpr T coef[] =612 {613 1.8989249102715535257L,614 0.95052179461824443490L,615 1.1510775899590158079L,616 1.7502391069863005399L,617 2.9526768126368751802L,618 5.2858003961214508892L,619 9.8324857166599797471L,620 18.787148683275595622L,621 36.614686152736981447L,622 72.452923951277718013L,623 145.10795773470691023L,624 293.47863963084970259L,625 598.38518150550101790L,626 1228.4200130758634505L,627 2536.5297553827644880L,628 5263.9832725075189576L,629 10972.138126273491753L,630 22958.388550988306870L,631 48203.103373625406989L632 };633 return boost::math::tools::evaluate_polynomial(coef, m - 0.55L);634 }635 case 12:636 case 13:637 {638 constexpr T coef[] =639 {640 2.0075983984243763017L,641 1.2484572312123473371L,642 1.9262346570764797287L,643 3.7512896400875876798L,644 8.1199445549320458022L,645 18.665721308735553611L,646 44.603924842914370633L,647 109.50920543094983774L,648 274.27795482324134804L,649 697.55980086063261629L,650 1795.7160145002471293L,651 4668.3817167903899100L,652 12235.762468136643348L,653 32290.178097183208178L,654 85713.076081959646847L,655 228672.18904931170958L,656 612757.27119158527740L,657 1.6483233976504668314e6L,658 4.4492251046211960936e6L,659 1.2046317340783185238e7L,660 3.2705187507963254185e7L661 };662 return boost::math::tools::evaluate_polynomial(coef, m - 0.65L);663 }664 case 14:665 case 15:666 {667 constexpr T coef[] =668 {669 2.1565156474996432354L,670 1.7918056418494632425L,671 3.8267512874657131470L,672 10.386724683637972080L,673 31.403314054680702901L,674 100.92370394986954165L,675 337.32682826322728966L,676 1158.7079305678279173L,677 4060.9907421936320917L,678 14454.001840343447947L,679 52076.661075994048028L,680 189493.65914621568866L,681 695184.57624138961450L,682 2.5679940482552846861e6L,683 9.5419219667483863221e6L,684 3.5634927442180761743e7L,685 1.3366929846120408712e8L,686 5.0335218668662845411e8L,687 1.9019757295386601192e9L,688 7.2089150153301037563e9L,689 2.7398741806339510931e10L,690 1.0439286724885300495e11L,691 3.9864875581513728207e11L,692 1.5254661585564745591e12L,693 5.8483259088850315936e12694 };695 return boost::math::tools::evaluate_polynomial(coef, m - 0.75L);696 }697 case 16:698 {699 constexpr T coef[] =700 {701 2.3181226217125105894L,702 2.6169201502912328409L,703 7.8979350757313558232L,704 30.502397154466723270L,705 131.48693655235284561L,706 602.98476373564916170L,707 2877.0246178099726410L,708 14110.519919151803247L,709 70621.440881565402289L,710 358977.26658253099258L,711 1.8472382637239716844e6L,712 9.6005154160492141090e6L,713 5.0307677085023668786e7L,714 2.6544418865271279673e8L,715 1.4088623250287026866e9L,716 7.5156879353737746270e9L,717 4.0270783964955246149e10L,718 2.1662089325801126339e11L,719 1.1692489201929996116e12L,720 6.3306543358985679881e12721 };722 return boost::math::tools::evaluate_polynomial(coef, m - 0.825L);723 }724 case 17:725 {726 constexpr T coef[] =727 {728 2.4735961737513439120L,729 3.7276242441180993105L,730 15.607393035549304964L,731 84.128508428058877470L,732 506.98181970406139349L,733 3252.2770581451236438L,734 21713.242419574342564L,735 149037.04518909327662L,736 1.0439993310899908390e6L,737 7.4279748170420389947e6L,738 5.3503839675586611510e7L,739 3.8924988699487084738e8L,740 2.8552883511008106195e9L,741 2.1090077038766840525e10L,742 1.5669983394779020136e11L,743 1.1702222424224398927e12L,744 8.7779483236689379709e12L,745 6.6101242752484950408e13L,746 4.9948805371338879891e14L,747 3.7859743397240299201e15L,748 2.8775996123036112296e16L,749 2.1926346839925760143e17L,750 1.6744985438468349361e18L,751 1.2814410112866546052e19L,752 9.8249807041031260167e19753 };754 return boost::math::tools::evaluate_polynomial(coef, m - 0.875L);755 }756 default:757 //758 // All cases where m > 0.9759 // including all error handling:760 //761 return ellint_k_imp(k, pol, boost::math::integral_constant<int, 2>());762 }763}764 765template <typename T, typename Policy>766BOOST_MATH_GPU_ENABLED typename tools::promote_args<T>::type ellint_1(T k, const Policy& pol, const boost::math::true_type&)767{768 typedef typename tools::promote_args<T>::type result_type;769 typedef typename policies::evaluation<result_type, Policy>::type value_type;770 typedef boost::math::integral_constant<int, 771#if defined(__clang_major__) && (__clang_major__ == 7)772 2773#else774 boost::math::is_floating_point<T>::value && boost::math::numeric_limits<T>::digits && (boost::math::numeric_limits<T>::digits <= 54) ? 0 :775 boost::math::is_floating_point<T>::value && boost::math::numeric_limits<T>::digits && (boost::math::numeric_limits<T>::digits <= 64) ? 1 : 2776#endif777 > precision_tag_type;778 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%)");779}780 781template <class T1, class T2>782BOOST_MATH_GPU_ENABLED typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi, const boost::math::false_type&)783{784 return boost::math::ellint_1(k, phi, policies::policy<>());785}786 787} // namespace detail788 789// Elliptic integral (Legendre form) of the first kind790template <class T1, class T2, class Policy>791BOOST_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!792{793 typedef typename tools::promote_args<T1, T2>::type result_type;794 typedef typename policies::evaluation<result_type, Policy>::type value_type;795 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%)");796}797 798template <class T1, class T2>799BOOST_MATH_GPU_ENABLED typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi)800{801 typedef typename policies::is_policy<T2>::type tag_type;802 return detail::ellint_1(k, phi, tag_type());803}804 805// Complete elliptic integral (Legendre form) of the first kind806template <typename T>807BOOST_MATH_GPU_ENABLED typename tools::promote_args<T>::type ellint_1(T k)808{809 return ellint_1(k, policies::policy<>());810}811 812}} // namespaces813 814#endif // BOOST_MATH_ELLINT_1_HPP815 816