pulp::signal::special Namespace Reference

Functions
double	elliptic_K (double m)
	Complete elliptic integral of the first kind, K(m), m = k^2 in [0, 1).
double	elliptic_E (double m)
	Complete elliptic integral of the second kind, E(m), m = k^2 in [0, 1].
double	jacobi_nome (double m)
	Jacobi nome q(m) = exp(-pi * K(1 - m) / K(m)).
void	jacobi_sncndn (double u, double m, double *sn_out, double *cn_out, double *dn_out)
	Jacobi elliptic functions sn(u, m), cn(u, m), dn(u, m) — all three at once because the AGM-based scheme produces them together.
double	jacobi_sn (double u, double m)
	Convenience: sn(u, m).
double	jacobi_cn (double u, double m)
	Convenience: cn(u, m).
double	jacobi_dn (double u, double m)
	Convenience: dn(u, m).
double	carlson_RF (double x, double y, double z)
	Carlson's symmetric elliptic integral of the first kind, RF(x, y, z).
double	jacobi_asn (double x, double m)
	Inverse Jacobi sn: returns u such that sn(u, m) = x, for x in [-1, 1].

Variables
constexpr double	kPi = 3.141592653589793238462643383279502884
	Pi as a double-precision constant. Used by every routine below.

Function Documentation

carlson_RF()

double pulp::signal::special::carlson_RF ( double x, double y, double z )

inline

Carlson's symmetric elliptic integral of the first kind, RF(x, y, z).

Used internally by jacobi_asn(). Carlson's duplication theorem reduces the integral to a power series around the average of the three arguments, converging to machine precision in ~5 iterations for well-conditioned inputs. Reference: B.C. Carlson, "Numerical Computation of Real or Complex Elliptic Integrals" (1995), and Numerical Recipes 3rd ed. §6.12 (function rf).

Domain: x, y, z >= 0 with at most one of them zero. Caller is responsible for staying inside this domain; we simply clamp negatives to zero so a stray rounding error in jacobi_asn() doesn't escape as NaN into a coefficient stream.

elliptic_E()

double pulp::signal::special::elliptic_E ( double m )

inline

Complete elliptic integral of the second kind, E(m), m = k^2 in [0, 1].

AGM-based: track c_n = (a_n - b_n)/2 alongside the AGM iteration; then E(m) = K(m) * (1 - sum_{n=0}^{inf} 2^{n-1} c_n^2) where c_0 = sqrt(m). The c_n sequence shrinks quadratically so a handful of terms suffices. The sum coefficient pattern is 2^{-1}, 2^0, 2^1, ... so we seed pow2 = 1/2 before doubling each step.

E(0) = pi/2 and E(1) = 1 exactly (returned as edge cases).

elliptic_K()

double pulp::signal::special::elliptic_K ( double m )

inline

Complete elliptic integral of the first kind, K(m), m = k^2 in [0, 1).

Computed via the arithmetic-geometric mean (AGM): start with (a0, b0) = (1, sqrt(1 - m)); iterate (a_{n+1}, b_{n+1}) = ((a + b)/2, sqrt(a * b)); K(m) = pi / (2 * AGM). Converges quadratically; 8 to 12 iterations reach machine precision for m up to ~1 - 1e-12.

At m = 1 the integral diverges; the implementation returns +infinity. At m = 0 it returns pi/2 exactly.

jacobi_asn()

double pulp::signal::special::jacobi_asn ( double x, double m )

inline

Inverse Jacobi sn: returns u such that sn(u, m) = x, for x in [-1, 1].

Implementation: asn(x, m) = x * RF(1 - x^2, 1 - m * x^2, 1). This is the standard Carlson form (Numerical Recipes §6.12), accurate to machine precision across the full m ∈ [0, 1] range. The sign of x carries through naturally since RF is even in its arguments and we only multiply by x.

For m = 0 this reduces to asin(x); for m = 1 it reduces to atanh(x). Inputs outside [-1, 1] are clamped to keep filter-design code numerically safe rather than throwing or NaN-poisoning the cascade.

jacobi_cn()

double pulp::signal::special::jacobi_cn ( double u, double m )

inline

Convenience: cn(u, m).

jacobi_dn()

double pulp::signal::special::jacobi_dn ( double u, double m )

inline

Convenience: dn(u, m).

jacobi_nome()

double pulp::signal::special::jacobi_nome ( double m )

inline

Jacobi nome q(m) = exp(-pi * K(1 - m) / K(m)).

q is the natural argument for theta-function series and the design of elliptic rational functions. For small m the nome is approximately m/16; for m close to 1 it approaches 1 from below.

jacobi_sn()

double pulp::signal::special::jacobi_sn ( double u, double m )

inline

Convenience: sn(u, m).

jacobi_sncndn()

void pulp::signal::special::jacobi_sncndn ( double u, double m, double * sn_out, double * cn_out, double * dn_out )

inline

Jacobi elliptic functions sn(u, m), cn(u, m), dn(u, m) — all three at once because the AGM-based scheme produces them together.

Algorithm: Numerical Recipes 3rd ed. §6.12 (Press et al.). We run the AGM until the two sequences a_n and emc_n agree, then unwind the recursion using the recurrence dn_{n-1} = (em_n + dn_n * a) / (b + a), a := c / b to recover sn / cn / dn at the original modulus. This converges to machine precision in ~6 iterations across the entire m ∈ (0, 1) range without needing per-step asin() calls.

Outputs are written through pointer arguments so callers can request any subset cheaply (pass nullptr for outputs they don't need).

Variable Documentation

kPi

double pulp::signal::special::kPi = 3.141592653589793238462643383279502884

inlineconstexpr

Pi as a double-precision constant. Used by every routine below.