ellpj

Jacobian Elliptic Functions.

SYNOPSIS:

double u, m, sn, cn, dn, phi;
int ellpj();
ellpj(u, m, _&sn, _&cn, _&dn, _&phi);

DESCRIPTION:
Evaluates the Jacobian elliptic functions sn(u|m), cn(u|m), and dn(u|m) of parameter m between 0 and 1, and real argument u.
These functions are periodic, with quarter-period on the real axis equal to the complete elliptic integral ellpk(1.0-m).

Relation to incomplete elliptic integral:
If u = ellik(phi,m), then sn(u|m) = sin(phi), and cn(u|m) = cos(phi).
Phi is called the amplitude of u.

Computation is by means of the arithmetic-geometric mean algorithm, except when m is within 1e-9 of 0 or 1.
In the latter case with m close to 1, the approximation applies only for phi < pi/2.

ACCURACY:

Tested at random points with u between 0 and 10, m between 0 and 1.

            Absolute error (* = relative error):
arithmetic   function   # trials      peak         rms
    DEC       sn           1800       4.5e-16     8.7e-17
    IEEE      phi         10000       9.2e-16*    1.4e-16*
    IEEE      sn          50000       4.1e-15     4.6e-16
    IEEE      cn          40000       3.6e-15     4.4e-16
    IEEE      dn         100000       3.9e-15     1.7e-16

Larger errors occur for m near 1.
Peak error observed in consistency check using addition theorem for sn(u+v) was 4e-16 (absolute). Also tested by the earlier relation to the incomplete elliptic integral. Accuracy deteriorates when u is large.