Module jacobi_elliptic 

The Jacobi elliptic functions sn, cn and dn.

What they are

All three are built from one quantity, the amplitude $\varphi = \mathrm{am}(u, k)$, defined as the angle whose incomplete elliptic integral of the first kind is u ($F(\varphi, k) = u$, see ellint_impl). Then

\mathrm{sn}(u, k) = \sin\varphi, \qquad
\mathrm{cn}(u, k) = \cos\varphi, \qquad
\mathrm{dn}(u, k) = \sqrt{1 - k^2\sin^2\varphi}

Hence the names: sine amplitude, cosine amplitude, delta amplitude. At k = 0 the amplitude is u itself and they degenerate to sin u, cos u and 1. At k = 1 they become tanh u, sech u and sech u.

They are returned together because they are a closed system, not merely because it is cheaper: differentiating any one of them produces a product of the other two ($\mathrm{sn}' = \mathrm{cn}\,\mathrm{dn}$, $\mathrm{cn}' = -\mathrm{sn}\,\mathrm{dn}$, $\mathrm{dn}' = -k^2\mathrm{sn}\,\mathrm{cn}$), exactly the way sin and cos close under differentiation. The other nine Jacobi functions in Glaisher’s notation (ns, nc, nd, sc, sd, cs, cd, ds, dc) are reciprocals and ratios of these three, so a caller holding the triple holds all twelve.

Algorithm

The descending Landen transformation, in the arithmetic-only form due to Bulirsch (1965) rather than the textbook one. Both walk the same AGM ladder down from k to modulus zero and then climb back up, but they differ in what the climb costs:

• The textbook descent (A&S 16.4, and Boost’s jacobi_elliptic) carries an angle back up, $\varphi_{n-1} = \tfrac12(\varphi_n + \arcsin(\tfrac{c_n}{a_n}\sin\varphi_n))$. That is one sin and one asin per level, on a ladder several levels deep, the worst possible shape for a vector unit, where every lane pays for both.
• Bulirsch carries the tangent of the angle instead. The half-angle step becomes rational, so the entire climb is multiplies and divides, and the whole function needs exactly one sin_cos, at the bottom of the ladder where the modulus is zero and the amplitude is just the argument.

Measured against mpmath at 40 digits over k in [0, 1) and |u| <= 8, worst absolute error 8.3 eps for sn, 4.1 for cn, 3.8 for dn. Absolute is the honest metric here: all three are bounded by 1 and all three have zeros, so relative error at a zero is governed by how well the zero’s location is known, exactly as for sin.

Accuracy degrades with |u| the way sin’s does and for the same reason: the one trig call takes u scaled by the AGM limit, so a large |u| is a large argument to reduce. The error above was measured to |u| = 8, and grows slowly beyond that.

The ladder is bounded, and short

$k' = \sqrt{1 - k^2}$ is what the AGM starts from, and for any k strictly below 1 in binary64 the cancellation-free one_minus_sq bottoms out at $2^{-52}$, so k' never falls below about 1.5e-8 and the ladder is never deeper than 8 rungs (measured, 4 to 6 is typical). NMAX carries two rungs of margin on top of that.

Lanes converge at different depths, so the loop runs until every lane has converged and the climb then runs the full depth for all of them. That is safe: past convergence $a_n = b_n$, so the extra rungs are identity transformations, and running them for every lane unconditionally was measured to give bit-identical results to stopping each lane at its own depth.

Constants

NMAX

Maximum depth of the AGM ladder.

Functions

jacobi_elliptic

(sn, cn, dn) at argument u and modulus k.