1[/ 2Copyright (c) 2006 Xiaogang Zhang 3Copyright (c) 2006 John Maddock 4Use, modification and distribution are subject to the 5Boost Software License, Version 1.0. (See accompanying file 6LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) 7] 8 9[section:ellint_1 Elliptic Integrals of the First Kind - Legendre Form] 10 11[heading Synopsis] 12 13`` 14 #include <boost/math/special_functions/ellint_1.hpp> 15`` 16 17 namespace boost { namespace math { 18 19 template <class T1, class T2> 20 ``__sf_result`` ellint_1(T1 k, T2 phi); 21 22 template <class T1, class T2, class ``__Policy``> 23 ``__sf_result`` ellint_1(T1 k, T2 phi, const ``__Policy``&); 24 25 template <class T> 26 ``__sf_result`` ellint_1(T k); 27 28 template <class T, class ``__Policy``> 29 ``__sf_result`` ellint_1(T k, const ``__Policy``&); 30 31 }} // namespaces 32 33[heading Description] 34 35These two functions evaluate the incomplete elliptic integral of the first kind 36['F([phi], k)] and its complete counterpart ['K(k) = F([pi]/2, k)]. 37 38[graph ellint_1] 39 40The return type of these functions is computed using the __arg_promotion_rules 41when T1 and T2 are different types: when they are the same type then the result 42is the same type as the arguments. 43 44 template <class T1, class T2> 45 ``__sf_result`` ellint_1(T1 k, T2 phi); 46 47 template <class T1, class T2, class ``__Policy``> 48 ``__sf_result`` ellint_1(T1 k, T2 phi, const ``__Policy``&); 49 50Returns the incomplete elliptic integral of the first kind ['F([phi], k)]: 51 52[equation ellint2] 53 54Requires k[super 2]sin[super 2](phi) < 1, otherwise returns the result of __domain_error. 55 56[optional_policy] 57 58 template <class T> 59 ``__sf_result`` ellint_1(T k); 60 61 template <class T> 62 ``__sf_result`` ellint_1(T k, const ``__Policy``&); 63 64Returns the complete elliptic integral of the first kind ['K(k)]: 65 66[equation ellint6] 67 68Requires |k| < 1, otherwise returns the result of __domain_error. 69 70[optional_policy] 71 72[heading Accuracy] 73 74These functions are computed using only basic arithmetic operations, so 75there isn't much variation in accuracy over differing platforms. 76Note that only results for the widest floating point type on the 77system are given as narrower types have __zero_error. All values 78are relative errors in units of epsilon. 79 80[table_ellint_1] 81 82The following error plot are based on an exhaustive search of the functions domain, MSVC-15.5 at `double` precision, 83and GCC-7.1/Ubuntu for `long double` and `__float128`. 84 85[graph elliptic_integral_k__double] 86 87[graph elliptic_integral_k__80_bit_long_double] 88 89[graph elliptic_integral_k____float128] 90 91[heading Testing] 92 93The tests use a mixture of spot test values calculated using the online 94calculator at [@http://functions.wolfram.com/ functions.wolfram.com], 95and random test data generated using 96NTL::RR at 1000-bit precision and this implementation. 97 98[heading Implementation] 99 100These functions are implemented in terms of Carlson's integrals using the relations: 101 102[equation ellint19] 103 104and 105 106[equation ellint20] 107 108[endsect] [/section:ellint_1 Elliptic Integrals of the First Kind - Legendre Form] 109 110[section:ellint_2 Elliptic Integrals of the Second Kind - Legendre Form] 111 112[heading Synopsis] 113 114`` 115 #include <boost/math/special_functions/ellint_2.hpp> 116`` 117 118 namespace boost { namespace math { 119 120 template <class T1, class T2> 121 ``__sf_result`` ellint_2(T1 k, T2 phi); 122 123 template <class T1, class T2, class ``__Policy``> 124 ``__sf_result`` ellint_2(T1 k, T2 phi, const ``__Policy``&); 125 126 template <class T> 127 ``__sf_result`` ellint_2(T k); 128 129 template <class T, class ``__Policy``> 130 ``__sf_result`` ellint_2(T k, const ``__Policy``&); 131 132 }} // namespaces 133 134[heading Description] 135 136These two functions evaluate the incomplete elliptic integral of the second kind 137['E([phi], k)] and its complete counterpart ['E(k) = E([pi]/2, k)]. 138 139[graph ellint_2] 140 141The return type of these functions is computed using the __arg_promotion_rules 142when T1 and T2 are different types: when they are the same type then the result 143is the same type as the arguments. 144 145 template <class T1, class T2> 146 ``__sf_result`` ellint_2(T1 k, T2 phi); 147 148 template <class T1, class T2, class ``__Policy``> 149 ``__sf_result`` ellint_2(T1 k, T2 phi, const ``__Policy``&); 150 151Returns the incomplete elliptic integral of the second kind ['E([phi], k)]: 152 153[equation ellint3] 154 155Requires k[super 2]sin[super 2](phi) < 1, otherwise returns the result of __domain_error. 156 157[optional_policy] 158 159 template <class T> 160 ``__sf_result`` ellint_2(T k); 161 162 template <class T> 163 ``__sf_result`` ellint_2(T k, const ``__Policy``&); 164 165Returns the complete elliptic integral of the second kind ['E(k)]: 166 167[equation ellint7] 168 169Requires |k| < 1, otherwise returns the result of __domain_error. 170 171[optional_policy] 172 173[heading Accuracy] 174 175These functions are computed using only basic arithmetic operations, so 176there isn't much variation in accuracy over differing platforms. 177Note that only results for the widest floating point type on the 178system are given as narrower types have __zero_error. All values 179are relative errors in units of epsilon. 180 181[table_ellint_2] 182 183The following error plot are based on an exhaustive search of the functions domain, MSVC-15.5 at `double` precision, 184and GCC-7.1/Ubuntu for `long double` and `__float128`. 185 186[graph elliptic_integral_e__double] 187 188[graph elliptic_integral_e__80_bit_long_double] 189 190[graph elliptic_integral_e____float128] 191 192[heading Testing] 193 194The tests use a mixture of spot test values calculated using the online 195calculator at [@http://functions.wolfram.com 196functions.wolfram.com], and random test data generated using 197NTL::RR at 1000-bit precision and this implementation. 198 199[heading Implementation] 200 201These functions are implemented in terms of Carlson's integrals 202using the relations: 203 204[equation ellint21] 205 206and 207 208[equation ellint22] 209 210[endsect] [/section:ellint_2 Elliptic Integrals of the Second Kind - Legendre Form] 211 212[section:ellint_3 Elliptic Integrals of the Third Kind - Legendre Form] 213 214[heading Synopsis] 215 216`` 217 #include <boost/math/special_functions/ellint_3.hpp> 218`` 219 220 namespace boost { namespace math { 221 222 template <class T1, class T2, class T3> 223 ``__sf_result`` ellint_3(T1 k, T2 n, T3 phi); 224 225 template <class T1, class T2, class T3, class ``__Policy``> 226 ``__sf_result`` ellint_3(T1 k, T2 n, T3 phi, const ``__Policy``&); 227 228 template <class T1, class T2> 229 ``__sf_result`` ellint_3(T1 k, T2 n); 230 231 template <class T1, class T2, class ``__Policy``> 232 ``__sf_result`` ellint_3(T1 k, T2 n, const ``__Policy``&); 233 234 }} // namespaces 235 236[heading Description] 237 238These two functions evaluate the incomplete elliptic integral of the third kind 239['[Pi](n, [phi], k)] and its complete counterpart ['[Pi](n, k) = E(n, [pi]/2, k)]. 240 241[graph ellint_3] 242 243The return type of these functions is computed using the __arg_promotion_rules 244when the arguments are of different types: when they are the same type then the result 245is the same type as the arguments. 246 247 template <class T1, class T2, class T3> 248 ``__sf_result`` ellint_3(T1 k, T2 n, T3 phi); 249 250 template <class T1, class T2, class T3, class ``__Policy``> 251 ``__sf_result`` ellint_3(T1 k, T2 n, T3 phi, const ``__Policy``&); 252 253Returns the incomplete elliptic integral of the third kind ['[Pi](n, [phi], k)]: 254 255[equation ellint4] 256 257Requires ['k[super 2]sin[super 2](phi) < 1] and ['n < 1/sin[super 2]([phi])], otherwise 258returns the result of __domain_error (outside this range the result 259would be complex). 260 261[optional_policy] 262 263 template <class T1, class T2> 264 ``__sf_result`` ellint_3(T1 k, T2 n); 265 266 template <class T1, class T2, class ``__Policy``> 267 ``__sf_result`` ellint_3(T1 k, T2 n, const ``__Policy``&); 268 269Returns the complete elliptic integral of the first kind ['[Pi](n, k)]: 270 271[equation ellint8] 272 273Requires ['|k| < 1] and ['n < 1], otherwise returns the 274result of __domain_error (outside this range the result would be complex). 275 276[optional_policy] 277 278[heading Accuracy] 279 280These functions are computed using only basic arithmetic operations, so 281there isn't much variation in accuracy over differing platforms. 282Note that only results for the widest floating point type on the 283system are given as narrower types have __zero_error. All values 284are relative errors in units of epsilon. 285 286[table_ellint_3] 287 288[heading Testing] 289 290The tests use a mixture of spot test values calculated using the online 291calculator at [@http://functions.wolfram.com 292functions.wolfram.com], and random test data generated using 293NTL::RR at 1000-bit precision and this implementation. 294 295[heading Implementation] 296 297The implementation for [Pi](n, [phi], k) first siphons off the special cases: 298 299[expression ['[Pi](0, [phi], k) = F([phi], k)]] 300 301[expression ['[Pi](n, [pi]/2, k) = [Pi](n, k)]] 302 303and 304 305[equation ellint23] 306 307Then if n < 0 the relations (A&S 17.7.15/16): 308 309[equation ellint24] 310 311are used to shift /n/ to the range \[0, 1\]. 312 313Then the relations: 314 315[expression ['[Pi](n, -[phi], k) = -[Pi](n, [phi], k)]] 316 317[expression ['[Pi](n, [phi]+m[pi], k) = [Pi](n, [phi], k) + 2m[Pi](n, k) ; n <= 1]] 318 319[expression ['[Pi](n, [phi]+m[pi], k) = [Pi](n, [phi], k) ; n > 1] [indent] [indent] 320[footnote I haven't been able to find a literature reference for this 321relation, but it appears to be the convention used by Mathematica. 322Intuitively the first ['2 * m * [Pi](n, k)] terms cancel out as the 323derivative alternates between +[infin] and -[infin].]] 324 325are used to move [phi] to the range \[0, [pi]\/2\]. 326 327The functions are then implemented in terms of Carlson's integrals using the relations: 328 329[equation ellint25] 330 331and 332 333[equation ellint26] 334 335[endsect] [/section:ellint_3 Elliptic Integrals of the Third Kind - Legendre Form] 336 337[section:ellint_d Elliptic Integral D - Legendre Form] 338 339[heading Synopsis] 340 341`` 342 #include <boost/math/special_functions/ellint_d.hpp> 343`` 344 345 namespace boost { namespace math { 346 347 template <class T1, class T2> 348 ``__sf_result`` ellint_d(T1 k, T2 phi); 349 350 template <class T1, class T2, class ``__Policy``> 351 ``__sf_result`` ellint_d(T1 k, T2 phi, const ``__Policy``&); 352 353 template <class T1> 354 ``__sf_result`` ellint_d(T1 k); 355 356 template <class T1, class ``__Policy``> 357 ``__sf_result`` ellint_d(T1 k, const ``__Policy``&); 358 359 }} // namespaces 360 361[heading Description] 362 363These two functions evaluate the incomplete elliptic integral 364['D([phi], k)] and its complete counterpart ['D(k) = D([pi]/2, k)]. 365 366The return type of these functions is computed using the __arg_promotion_rules 367when the arguments are of different types: when they are the same type then the result 368is the same type as the arguments. 369 370 template <class T1, class T2> 371 ``__sf_result`` ellint_d(T1 k, T2 phi); 372 373 template <class T1, class T2, class ``__Policy``> 374 ``__sf_result`` ellint_3(T1 k, T2 phi, const ``__Policy``&); 375 376Returns the incomplete elliptic integral: 377 378[equation ellint_d] 379 380Requires ['k[super 2]sin[super 2](phi) < 1], otherwise 381returns the result of __domain_error (outside this range the result 382would be complex). 383 384[optional_policy] 385 386 template <class T1> 387 ``__sf_result`` ellint_d(T1 k); 388 389 template <class T1, class ``__Policy``> 390 ``__sf_result`` ellint_d(T1 k, const ``__Policy``&); 391 392Returns the complete elliptic integral ['D(k) = D([pi]/2, k)] 393 394Requires ['-1 <= k <= 1] otherwise returns the 395result of __domain_error (outside this range the result would be complex). 396 397[optional_policy] 398 399[heading Accuracy] 400 401These functions are trivially computed in terms of other elliptic integrals 402and generally have very low error rates (a few epsilon) unless parameter [phi] 403is very large, in which case the usual trigonometric function argument-reduction issues apply. 404 405[table_ellint_d_complete_] 406 407[table_ellint_d] 408 409The following error plot are based on an exhaustive search of the functions domain, MSVC-15.5 at `double` precision, 410and GCC-7.1/Ubuntu for `long double` and `__float128`. 411 412[graph elliptic_integral_d__double] 413 414[graph elliptic_integral_d__80_bit_long_double] 415 416[graph elliptic_integral_d____float128] 417 418 419[heading Testing] 420 421The tests use a mixture of spot test values calculated using 422values calculated at __WolframAlpha, and random test data generated using 423MPFR at 1000-bit precision and a deliberately naive implementation in terms of 424the Legendre integrals. 425 426[heading Implementation] 427 428The implementation for D([phi], k) first performs argument reduction using the relations: 429 430[expression ['D(-[phi], k) = -D([phi], k)]] 431 432and 433 434[expression ['D(n[pi]+[phi], k) = 2nD(k) + D([phi], k)]] 435 436to move [phi] to the range \[0, [pi]\/2\]. 437 438The functions are then implemented in terms of Carlson's integral R[sub D] 439using the relation: 440 441[equation ellint_d] 442 443[endsect] [/section:ellint_d Elliptic Integral D - Legendre Form] 444 445[section:jacobi_zeta Jacobi Zeta Function] 446 447[heading Synopsis] 448 449`` 450 #include <boost/math/special_functions/jacobi_zeta.hpp> 451`` 452 453 namespace boost { namespace math { 454 455 template <class T1, class T2> 456 ``__sf_result`` jacobi_zeta(T1 k, T2 phi); 457 458 template <class T1, class T2, class ``__Policy``> 459 ``__sf_result`` jacobi_zeta(T1 k, T2 phi, const ``__Policy``&); 460 461 }} // namespaces 462 463[heading Description] 464 465This function evaluates the Jacobi Zeta Function ['Z([phi], k)] 466 467[equation jacobi_zeta] 468 469Please note the use of [phi], and /k/ as the parameters, the function is often defined as ['Z([phi], m)] 470with ['m = k[super 2]], see for example [@http://mathworld.wolfram.com/JacobiZetaFunction.html Weisstein, Eric W. "Jacobi Zeta Function." From MathWorld--A Wolfram Web Resource.] 471Or else as [@https://dlmf.nist.gov/22.16#E32 ['Z(x, k)]] with ['[phi] = am(x, k)], 472where ['am] is the [@https://dlmf.nist.gov/22.16#E1 Jacobi amplitude function] 473which is equivalent to ['asin(jacobi_elliptic(k, x))]. 474 475The return type of this function is computed using the __arg_promotion_rules 476when the arguments are of different types: when they are the same type then the result 477is the same type as the arguments. 478 479Requires ['-1 <= k <= 1], otherwise 480returns the result of __domain_error (outside this range the result would be complex). 481 482[optional_policy] 483 484Note that there is no complete analogue of this function (where [phi] = [pi] / 2) 485as this takes the value 0 for all ['k]. 486 487[heading Accuracy] 488 489These functions are trivially computed in terms of other elliptic integrals 490and generally have very low error rates (a few epsilon) unless parameter [phi] 491is very large, in which case the usual trigonometric function argument-reduction issues apply. 492 493[table_jacobi_zeta] 494 495[heading Testing] 496 497The tests use a mixture of spot test values calculated using 498values calculated at __WolframAlpha, and random test data generated using 499MPFR at 1000-bit precision and a deliberately naive implementation in terms of 500the Legendre integrals. 501 502[heading Implementation] 503 504The implementation for Z([phi], k) first makes the argument [phi] positive using: 505 506[expression ['Z(-[phi], k) = -Z([phi], k)]] 507 508The function is then implemented in terms of Carlson's integral R[sub J] 509using the relation: 510 511[equation jacobi_zeta] 512 513There is one special case where the above relation fails: when ['k = 1], in that case 514the function simplifies to 515 516[expression ['Z([phi], 1) = sign(cos([phi])) sin([phi])]] 517 518[h5:jacobi_zeta_example Example] 519 520A simple example comparing use of __WolframAlpha with Boost.Math (including much higher precision using Boost.Multiprecision) 521is [@../../example/jacobi_zeta_example.cpp jacobi_zeta_example.cpp]. 522 523[endsect] [/section:jacobi_zeta Jacobi Zeta Function] 524 525[section:heuman_lambda Heuman Lambda Function] 526 527[heading Synopsis] 528 529`` 530 #include <boost/math/special_functions/heuman_lambda.hpp> 531`` 532 533 namespace boost { namespace math { 534 535 template <class T1, class T2> 536 ``__sf_result`` heuman_lambda(T1 k, T2 phi); 537 538 template <class T1, class T2, class ``__Policy``> 539 ``__sf_result`` heuman_lambda(T1 k, T2 phi, const ``__Policy``&); 540 541 }} // namespaces 542 543[heading Description] 544 545This function evaluates the Heuman Lambda Function ['[Lambda][sub 0]([phi], k)] 546 547[equation heuman_lambda] 548 549The return type of this function is computed using the __arg_promotion_rules 550when the arguments are of different types: when they are the same type then the result 551is the same type as the arguments. 552 553Requires ['-1 <= k <= 1], otherwise 554returns the result of __domain_error (outside this range the result would be complex). 555 556[optional_policy] 557 558Note that there is no complete analogue of this function (where [phi] = [pi] / 2) 559as this takes the value 1 for all ['k]. 560 561[heading Accuracy] 562 563These functions are trivially computed in terms of other elliptic integrals 564and generally have very low error rates (a few epsilon) unless parameter [phi] 565is very large, in which case the usual trigonometric function argument-reduction issues apply. 566 567[table_heuman_lambda] 568 569[heading Testing] 570 571The tests use a mixture of spot test values calculated using 572values calculated at __WolframAlpha, and random test data generated using 573MPFR at 1000-bit precision and a deliberately naive implementation in terms of 574the Legendre integrals. 575 576[heading Implementation] 577 578The function is then implemented in terms of Carlson's integrals R[sub J] and R[sub F] 579using the relation: 580 581[equation heuman_lambda] 582 583This relation fails for ['|[phi]| >= [pi]/2] in which case the definition in terms of the 584Jacobi Zeta is used. 585 586[endsect] [/section:heuman_lambda Heuman Lambda Function] 587 588 589