1[/ 2Copyright (c) 2012 John Maddock 3Use, modification and distribution are subject to the 4Boost Software License, Version 1.0. (See accompanying file 5LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) 6] 7 8[section:jacobi Jacobi Elliptic Functions] 9 10[section:jac_over Overview of the Jacobi Elliptic Functions] 11 12There are twelve Jacobi Elliptic functions, of which the three copolar functions ['sn], ['cn] and ['dn] are the most important 13as the other nine can be computed from these three 14[footnote [@http://en.wikipedia.org/wiki/Jacobi_elliptic_functions Wikipedia: Jacobi elliptic functions]] 15[footnote [@http://mathworld.wolfram.com/JacobiEllipticFunctions.html Weisstein, Eric W. "Jacobi Elliptic Functions." From MathWorld - A Wolfram Web Resource.]] 16[footnote [@http://dlmf.nist.gov/22 Digital Library of Mathematical Functions: Jacobian Elliptic Functions, Reinhardt, W. P., Walker, O. L.]]. 17 18These functions each take two arguments: a parameter, and a variable as described below. 19 20Like all elliptic functions these can be parameterised in a number of ways: 21 22* In terms of a parameter ['m]. 23* In terms of the elliptic modulus ['k] where ['m = k[super 2]]. 24* In terms of the modular angle [alpha], where ['m = sin[super 2][thin][alpha]]. 25 26In our implementation, these functions all take the elliptic modulus /k/ as the parameter. 27 28In addition the variable /u/ is sometimes expressed as an amplitude [phi], in our implementation we always use /u/. 29 30Finally note that our functions all take the elliptic modulus /k/ as the *first* argument - this is for alignment with 31the Elliptic Integrals (but is different from other implementations, for example Mathworks). 32 33A simple example comparing use of __WolframAlpha with Boost.Math (including much higher precision using Boost.Multiprecision) 34is [@../../example/jacobi_zeta_example.cpp jacobi_zeta_example.cpp]. 35 36There are twelve functions for computing the twelve individual Jacobi elliptic functions: __jacobi_cd, __jacobi_cn, __jacobi_cs, 37__jacobi_dc, __jacobi_dn, __jacobi_ds, __jacobi_nc, __jacobi_nd, __jacobi_ns, __jacobi_sc, __jacobi_sd and __jacobi_sn. 38 39They are all called as for example: 40 41 jacobi_cs(k, u); 42 43Note however that these individual functions are all really thin wrappers around the function __jacobi_elliptic which calculates 44the three copolar functions ['sn], ['cn] and ['dn] in a single function call. 45 46[tip If you need more than one of these functions 47for a given set of arguments, it's most efficient to use __jacobi_elliptic.] 48 49[endsect] [/section:jac_over Overview of the Jacobi Elliptic Functions] 50 51[section:jacobi_elliptic Jacobi Elliptic SN, CN and DN] 52 53[heading Synopsis] 54 55`` 56 #include <boost/math/special_functions/jacobi_elliptic.hpp> 57`` 58 59 namespace boost { namespace math { 60 61 template <class T, class U, class V> 62 ``__sf_result`` jacobi_elliptic(T k, U u, V* pcn, V* pdn); 63 64 template <class T, class U, class V, class Policy> 65 ``__sf_result`` jacobi_elliptic(T k, U u, V* pcn, V* pdn, const Policy&); 66 67 }} // namespaces 68 69[heading Description] 70 71The function __jacobi_elliptic calculates the three copolar Jacobi elliptic functions 72['sn(u, k)], ['cn(u, k)] and ['dn(u, k)]. The returned value is 73['sn(u, k)], and if provided, `*pcn` is 74set to ['cn(u, k)], and `*pdn` is set to ['dn(u, k)]. 75 76The functions are defined as follows, given: 77 78[equation jacobi1] 79 80The the angle ['[phi]] is called the ['amplitude] and: 81 82[equation jacobi2] 83 84[note 85 ['[phi]] is called the amplitude. 86 ['k] is called the elliptic modulus. 87] 88 89[caution Rather like other elliptic functions, the Jacobi functions 90are expressed in a variety of different ways. In particular, 91the parameter /k/ (the modulus) may also be expressed using a modular 92angle [alpha], or a parameter /m/. These are related by: 93 94[expression k = sin [alpha]] 95 96[expression m = k[super 2] = sin[super 2][alpha]] 97 98So that the function ['sn] (for example) may be expressed as 99either: 100 101[expression sn(u, k)] 102 103[expression sn(u \\ [alpha])] 104 105[expression sn(u | m)] 106 107To further complicate matters, some texts refer to the ['complement 108of the parameter m], or 1 - m, where: 109 110[expression 1 - m = 1 - k[super 2] = cos[super 2][alpha]] 111 112This implementation uses /k/ throughout, and makes this the first argument 113to the functions: this is for alignment with the elliptic integrals which match the requirements 114of the [@http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2005/n1836.pdf 115Technical Report on C++ Library Extensions]. However, you should 116be extra careful when using these functions!] 117 118[optional_policy] 119 120The following graphs illustrate how these functions change as /k/ changes: for small /k/ 121these are sine waves, while as /k/ tends to 1 they become hyperbolic functions: 122 123[graph jacobi_sn] 124 125[graph jacobi_cn] 126 127[graph jacobi_dn] 128 129[heading Accuracy] 130 131These functions are computed using only basic arithmetic operations and trigonometric functions, so 132there isn't much variation in accuracy over differing platforms. 133Typically errors are trivially small for small angles, and as is typical for cyclic 134functions, grow as the angle increases. 135Note that only results for the widest floating-point type on the 136system are given as narrower types have __zero_error. All values 137are relative errors in units of epsilon. 138 139[table_jacobi_cn] 140 141[table_jacobi_dn] 142 143[table_jacobi_sn] 144 145[heading Testing] 146 147The tests use a mixture of spot test values calculated using the online 148calculator at [@http://functions.wolfram.com/ functions.wolfram.com], 149and random test data generated using 150MPFR at 1000-bit precision and this implementation. 151 152[heading Implementation] 153 154For ['k > 1] we apply the relations: 155 156[equation jacobi3] 157 158Then filter off the special cases: 159 160[expression ['sn(0, k) = 0] and ['cn(0, k) = dn(0, k) = 1]] 161 162[expression ['sn(u, 0) = sin(u), cn(u, 0) = cos(u) and dn(u, 0) = 1]] 163 164[expression ['sn(u, 1) = tanh(u), cn(u, 1) = dn(u, 1) = 1 / cosh(u)]] 165 166And for ['k[super 4] < [epsilon]] we have: 167 168[equation jacobi4] 169 170Otherwise the values are calculated using the method of [@http://dlmf.nist.gov/22.20#SS2 arithmetic geometric means]. 171 172[endsect] [/section:jacobi_elliptic Jacobi Elliptic SN, CN and DN] 173 174 175[section:jacobi_cd Jacobi Elliptic Function cd] 176 177[heading Synopsis] 178 179`` 180 #include <boost/math/special_functions/jacobi_elliptic.hpp> 181`` 182 183 namespace boost { namespace math { 184 185 template <class T, class U> 186 ``__sf_result`` jacobi_cd(T k, U u); 187 188 template <class T, class U, class Policy> 189 ``__sf_result`` jacobi_cd(T k, U u, const Policy& pol); 190 191 }} // namespaces 192 193[heading Description] 194 195This function returns the Jacobi elliptic function ['cd]. 196 197[optional_policy] 198 199This function is a trivial wrapper around __jacobi_elliptic, with: 200 201[expression ['cd(u, k) = cn(u, k) / dn(u, k)]] 202 203[graph jacobi_cd] 204 205[endsect] [/section:jacobi_cd Jacobi Elliptic Function cd] 206 207 208[section:jacobi_cn Jacobi Elliptic Function cn] 209 210[heading Synopsis] 211 212`` 213 #include <boost/math/special_functions/jacobi_elliptic.hpp> 214`` 215 216 namespace boost { namespace math { 217 218 template <class T, class U> 219 ``__sf_result`` jacobi_cn(T k, U u); 220 221 template <class T, class U, class Policy> 222 ``__sf_result`` jacobi_cn(T k, U u, const Policy& pol); 223 224 }} // namespaces 225 226[heading Description] 227 228This function returns the Jacobi elliptic function ['cn]. 229 230[optional_policy] 231 232This function is a trivial wrapper around __jacobi_elliptic. 233 234[graph jacobi_cn] 235 236[endsect] [/section:jacobi_cn Jacobi Elliptic Function cn] 237 238 239[section:jacobi_cs Jacobi Elliptic Function cs] 240 241[heading Synopsis] 242 243`` 244 #include <boost/math/special_functions/jacobi_elliptic.hpp> 245`` 246 247 namespace boost { namespace math { 248 249 template <class T, class U> 250 ``__sf_result`` jacobi_cs(T k, U u); 251 252 template <class T, class U, class Policy> 253 ``__sf_result`` jacobi_cs(T k, U u, const Policy& pol); 254 255 }} // namespaces 256 257[heading Description] 258 259This function returns the Jacobi elliptic function ['cs]. 260 261[optional_policy] 262 263This function is a trivial wrapper around __jacobi_elliptic, with: 264 265[expression ['cs(u, k) = cn(u, k) / sn(u, k)]] 266 267[graph jacobi_cs] 268 269[endsect] [/section:jacobi_cs Jacobi Elliptic Function cs] 270 271 272[section:jacobi_dc Jacobi Elliptic Function dc] 273 274[heading Synopsis] 275 276`` 277 #include <boost/math/special_functions/jacobi_elliptic.hpp> 278`` 279 280 namespace boost { namespace math { 281 282 template <class T, class U> 283 ``__sf_result`` jacobi_dc(T k, U u); 284 285 template <class T, class U, class Policy> 286 ``__sf_result`` jacobi_dc(T k, U u, const Policy& pol); 287 288 }} // namespaces 289 290[heading Description] 291 292This function returns the Jacobi elliptic function ['dc]. 293 294[optional_policy] 295 296This function is a trivial wrapper around __jacobi_elliptic, with: 297 298[expression ['dc(u, k) = dn(u, k) / cn(u, k)]] 299 300[graph jacobi_dc] 301 302[endsect] [/section:jacobi_dc Jacobi Elliptic Function dc] 303 304[section:jacobi_dn Jacobi Elliptic Function dn] 305 306[heading Synopsis] 307 308`` 309 #include <boost/math/special_functions/jacobi_elliptic.hpp> 310`` 311 312 namespace boost { namespace math { 313 314 template <class T, class U> 315 ``__sf_result`` jacobi_dn(T k, U u); 316 317 template <class T, class U, class Policy> 318 ``__sf_result`` jacobi_dn(T k, U u, const Policy& pol); 319 320 }} // namespaces 321 322[heading Description] 323 324This function returns the Jacobi elliptic function ['dn]. 325 326[optional_policy] 327 328This function is a trivial wrapper around __jacobi_elliptic. 329 330[graph jacobi_dn] 331 332[endsect] 333 334[section:jacobi_ds Jacobi Elliptic Function ds] 335 336[heading Synopsis] 337 338`` 339 #include <boost/math/special_functions/jacobi_elliptic.hpp> 340`` 341 342 namespace boost { namespace math { 343 344 template <class T, class U> 345 ``__sf_result`` jacobi_ds(T k, U u); 346 347 template <class T, class U, class Policy> 348 ``__sf_result`` jacobi_ds(T k, U u, const Policy& pol); 349 350 }} // namespaces 351 352[heading Description] 353 354This function returns the Jacobi elliptic function ['ds]. 355 356[optional_policy] 357 358This function is a trivial wrapper around __jacobi_elliptic, with: 359 360[expression ['ds(u, k) = dn(u, k) / sn(u, k)]] 361 362[graph jacobi_ds] 363 364[endsect] [/section:jacobi_dn Jacobi Elliptic Function dn] 365 366[section:jacobi_nc Jacobi Elliptic Function nc] 367 368[heading Synopsis] 369 370`` 371 #include <boost/math/special_functions/jacobi_elliptic.hpp> 372`` 373 374 namespace boost { namespace math { 375 376 template <class T, class U> 377 ``__sf_result`` jacobi_nc(T k, U u); 378 379 template <class T, class U, class Policy> 380 ``__sf_result`` jacobi_nc(T k, U u, const Policy& pol); 381 382 }} // namespaces 383 384[heading Description] 385 386This function returns the Jacobi elliptic function ['nc]. 387 388[optional_policy] 389 390This function is a trivial wrapper around __jacobi_elliptic, with: 391 392[expression ['nc(u, k) = 1 / cn(u, k)]] 393 394[graph jacobi_nc] 395 396[endsect] [/section:jacobi_nc Jacobi Elliptic Function nc] 397 398[section:jacobi_nd Jacobi Elliptic Function nd] 399 400[heading Synopsis] 401 402`` 403 #include <boost/math/special_functions/jacobi_elliptic.hpp> 404`` 405 406 namespace boost { namespace math { 407 408 template <class T, class U> 409 ``__sf_result`` jacobi_nd(T k, U u); 410 411 template <class T, class U, class Policy> 412 ``__sf_result`` jacobi_nd(T k, U u, const Policy& pol); 413 414 }} // namespaces 415 416[heading Description] 417 418This function returns the Jacobi elliptic function ['nd]. 419 420[optional_policy] 421 422This function is a trivial wrapper around __jacobi_elliptic, with: 423 424[expression ['nd(u, k) = 1 / dn(u, k)]] 425 426[graph jacobi_nd] 427 428[endsect] [/section:jacobi_nd Jacobi Elliptic Function nd] 429 430[section:jacobi_ns Jacobi Elliptic Function ns] 431 432[heading Synopsis] 433 434`` 435 #include <boost/math/special_functions/jacobi_elliptic.hpp> 436`` 437 438 namespace boost { namespace math { 439 440 template <class T, class U> 441 ``__sf_result`` jacobi_ns(T k, U u); 442 443 template <class T, class U, class Policy> 444 ``__sf_result`` jacobi_ns(T k, U u, const Policy& pol); 445 446 }} // namespaces 447 448[heading Description] 449 450This function returns the Jacobi elliptic function ['ns]. 451 452[optional_policy] 453 454This function is a trivial wrapper around __jacobi_elliptic, with: 455 456[expression ['ns(u, k) = 1 / sn(u, k)]] 457 458[graph jacobi_ns] 459 460[endsect] [/section:jacobi_ns Jacobi Elliptic Function ns] 461 462[section:jacobi_sc Jacobi Elliptic Function sc] 463 464[heading Synopsis] 465 466`` 467 #include <boost/math/special_functions/jacobi_elliptic.hpp> 468`` 469 470 namespace boost { namespace math { 471 472 template <class T, class U> 473 ``__sf_result`` jacobi_sc(T k, U u); 474 475 template <class T, class U, class Policy> 476 ``__sf_result`` jacobi_sc(T k, U u, const Policy& pol); 477 478 }} // namespaces 479 480[heading Description] 481 482This function returns the Jacobi elliptic function ['sc]. 483 484[optional_policy] 485 486This function is a trivial wrapper around __jacobi_elliptic, with: 487 488[expression ['sc(u, k) = sn(u, k) / cn(u, k)]] 489 490[graph jacobi_sc] 491 492[endsect] [/section:jacobi_sc Jacobi Elliptic Function sc] 493 494[section:jacobi_sd Jacobi Elliptic Function sd] 495 496[heading Synopsis] 497 498`` 499 #include <boost/math/special_functions/jacobi_elliptic.hpp> 500`` 501 502 namespace boost { namespace math { 503 504 template <class T, class U> 505 ``__sf_result`` jacobi_sd(T k, U u); 506 507 template <class T, class U, class Policy> 508 ``__sf_result`` jacobi_sd(T k, U u, const Policy& pol); 509 510 }} // namespaces 511 512[heading Description] 513 514This function returns the Jacobi elliptic function ['sd]. 515 516[optional_policy] 517 518This function is a trivial wrapper around __jacobi_elliptic, with: 519 520[expression ['sd(u, k) = sn(u, k) / dn(u, k)]] 521 522[graph jacobi_sd] 523 524[endsect] [/section:jacobi_sd Jacobi Elliptic Function sd] 525 526[section:jacobi_sn Jacobi Elliptic Function sn] 527 528[heading Synopsis] 529 530`` 531 #include <boost/math/special_functions/jacobi_elliptic.hpp> 532`` 533 534 namespace boost { namespace math { 535 536 template <class T, class U> 537 ``__sf_result`` jacobi_sn(T k, U u); 538 539 template <class T, class U, class Policy> 540 ``__sf_result`` jacobi_sn(T k, U u, const Policy& pol); 541 542 }} // namespaces 543 544[heading Description] 545 546This function returns the Jacobi elliptic function ['sn]. 547 548[optional_policy] 549 550This function is a trivial wrapper around __jacobi_elliptic. 551 552[graph jacobi_sn] 553 554[endsect] [/section:jacobi_sn Jacobi Elliptic Function sn] 555 556[endsect] [/section:jacobi Jacobi Elliptic Functions] 557 558 559