ellint_2.hpp 25 KB

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