Class NDComplexMath
Complex-math helpers backing NumSharp's complex (complex128) unary ufuncs. Each entry
point reproduces NumPy 2.4.2 bit-for-bit (or within 3 ULP on the finite interior), verified by
a 504-point bit-exact sweep and a layout sweep (contiguous / F-contiguous / strided /
transposed / reversed / sliced / broadcast / 0-d). Most are direct ports of NumPy's own
routines in npy_math_complex.c.src (the FreeBSD msun implementations) because
System.Numerics.Complex diverges on large magnitudes, the unit circle, tiny/subnormal
values, branch cuts, and signed zeros.
Ported NumPy algorithms.
Log(Complex) = npy_clog (four-regime rescale incl. the near-|z|=1 log1p
path; drives Log10(Complex) and the engine's log2); Sinh(Complex)/
Cosh(Complex) = textbook sinh/cosh(x)·trig(y) with a y==0 guard (so a huge
real part doesn't become inf·0 = NaN) and the C99 Annex G non-finite tables;
Tanh(Complex) = Kahan's npy_ctanh (markedly more accurate than the BCL near
±π/2) + the |x|≥22 overflow-safe branch; Sin(Complex)/Cos(Complex)/
Tan(Complex) route through those exactly as NumPy defines csin/ccos/ctan;
Atan(Complex) = the full npy_catanh (real atanh/atan on the axes, the
log1p interior, and an exponent-classified real_part_reciprocal);
Exp(Complex) = npy_cexp; Sqrt(Complex) = npy_csqrt;
Expm1(Complex) = nc_expm1 with a Goldberg real expm1;
Square(Complex) = FMA-contracted z·z (matches NumPy's complex multiply
overflow/cancellation); Reciprocal(Complex) = Smith's nc_recip;
Exp2(Complex)/Log1p(Complex) compose the above; Abs(Complex) = npy_cabs
(C99 hypot: an infinite component yields +inf even alongside a NaN — the .NET 8
Complex.Abs returns NaN there).
Still delegating to the BCL (at parity): Asin(Complex) and Acos(Complex) use Asin(Complex)/Acos(Complex) on the finite interior with signed-zero / branch-cut fixups, and the C99 non-finite tables otherwise.
Accepted residuals (pathological inputs only, beyond 3 ULP): cos/sin with a
NaN imaginary part pick the C99-unspecified sign for the resulting zero; arccos
with a sub-DBL_MIN imaginary part flushes the denormal real part to 0 where NumPy's
cacos hard-work kernel keeps it (~5.8e-309); sinh/cosh at the |x|∈[710,710.13]
overflow edge differ because Windows' CRT sinh overflows where .NET's stays finite.
Perf: each public entry point is a tiny finite-path wrapper marked
AggressiveInlining so the JIT folds it into the IL-emitted
unary kernel (no per-element call frame); the rare non-finite / special-value tables live in
cold helpers marked AggressiveOptimization (kept out-of-line so
the hot wrapper stays inlineable, fully optimized when hit). A benchmark of an IL-inlined
variant vs this call-based form showed the per-element cost is dominated by the
transcendental, so hand-emitting the formulas is not worth the duplication.
public static class NDComplexMath
- Inheritance
-
NDComplexMath
- Inherited Members
Methods
Abs(Complex)
|z| = hypot(re, im) with NumPy/C99 infinity semantics: a ±infinite real or
imaginary part returns +inf regardless of the other part (including NaN).
All other inputs defer to Abs(Complex).
public static double Abs(Complex z)
Parameters
zComplex
Returns
Acos(Complex)
Complex arccosine matching NumPy. Finite inputs defer to Acos(Complex) with a
branch-cut fixup; non-finite inputs follow C99 (ported from npy_cacos).
public static Complex Acos(Complex z)
Parameters
zComplex
Returns
Asin(Complex)
Complex arcsine matching NumPy. Finite inputs defer to Asin(Complex) with
signed-zero/branch-cut fixups; non-finite inputs use the identity
asin(z) = iconj(casinh(iconj z)) where i*conj(z) = (y, x).
public static Complex Asin(Complex z)
Parameters
zComplex
Returns
Atan(Complex)
Complex arctangent matching NumPy. NumPy defines npy_catan(z) = iconj(catanh(iconj z))
verbatim (i*conj(z) = (y, x), then (w.Im, w.Re)), so routing through the full
Catanh(Complex) port reproduces NumPy bit-for-bit — including the tiny-imaginary and
subnormal cases where Atan(Complex) cancels (its internal log) or underflows to 0.
public static Complex Atan(Complex z)
Parameters
zComplex
Returns
Cos(Complex)
Complex cosine matching NumPy. NumPy defines npy_ccos(z) = cosh(iz) verbatim
(iz = (-y, x)), so routing through the C99-correct Cosh(Complex) reproduces
NumPy bit-for-bit, including the large-imaginary and signed-zero cases that
Cos(Complex) mishandles.
public static Complex Cos(Complex z)
Parameters
zComplex
Returns
Cosh(Complex)
Complex hyperbolic cosine matching NumPy. The y == 0 guard returns (cosh(x), x*y)
so a large x does not produce the BCL's sinh(x)sin(0) = inf0 = NaN imaginary
part (this is what lets Cos(Complex) handle a huge imaginary input); the general finite
case is the textbook cosh(x)cos(y) + i sinh(x)sin(y), which overflows to ±inf exactly
where NumPy's libm does (probed: |x| >= ~710.5). Non-finite inputs follow C99 Annex G.
public static Complex Cosh(Complex z)
Parameters
zComplex
Returns
Exp(Complex)
Complex exponential matching NumPy (npy_cexp). A finite real part defers to
Exp(Complex) (ULP-identical to NumPy); a non-finite real part follows C99.
public static Complex Exp(Complex z)
Parameters
zComplex
Returns
Exp2(Complex)
Complex base-2 exponential matching NumPy: exp2(z) = exp(z*ln2). Routing through the
C99-correct Exp(Complex) reproduces NumPy's non-finite results (e.g. exp2(+Inf+Inf i) =
+Inf + I NaN) that Pow(Complex, double) turned into (NaN, NaN).
public static Complex Exp2(Complex z)
Parameters
zComplex
Returns
Expm1(Complex)
Complex exp(z)-1 matching NumPy (nc_expm1):
real = expm1(x)cos(y) - 2sin^2(y/2), imag = exp(x)*sin(y). This reproduces
NumPy's structure (e.g. expm1(+Inf+0i).imag = exp(+Inf)*sin(0) = NaN) and signed
zeros. The BCL has no real expm1, so RealExpm1(double) falls back to
exp(x)-1 (with a signed-zero guard) — bit-identical to NumPy except for inputs
extremely close to the origin, where libm's expm1 is more accurate (a documented
divergence, like arctan's interior).
public static Complex Expm1(Complex z)
Parameters
zComplex
Returns
Log(Complex)
Complex natural logarithm matching NumPy (full npy_clog port, FreeBSD msun). The real
part is log|z| computed by rescaling to dodge the four problem regimes the naive
log(hypot(re,im)) hits: |z| huge (rescale by ½), subnormal (rescale by
2^53), near 1 (0.71 <= |z| <= 1.73 uses ½·log1p((m-1)(m+1)+n²) so the
real part doesn't cancel to 0), and 0 (handled by -1/re). The imaginary part is
atan2(im, re). Log(Complex) reproduces only the common regime.
public static Complex Log(Complex z)
Parameters
zComplex
Returns
Log10(Complex)
Complex base-10 logarithm matching NumPy (clog(z) * LOG10E), built from
Log(Complex) scaled by 1/ln(10). Log10(Complex) uses a different
scaling that drifts well past 1 ULP from NumPy.
public static Complex Log10(Complex z)
Parameters
zComplex
Returns
Log1p(Complex)
Complex log(1+z) matching NumPy's nc_log1p (funcs.inc.src). Computed as
real = log(hypot(re+1, im)), imag = atan2(im, re+1) — the naive magnitude
log, NOT Log(Complex)'s near-|z|=1 refinement (NumPy's np.log1p(1e-10j).real == 0,
not the 5e-21 a refined clog(1+z) would give). Routing through the C99
Hypot(double, double) (instead of the BCL Log(Complex), whose Complex.Abs
turns hypot(±inf, NaN) into NaN) reproduces NumPy's non-finite results — most
importantly log1p(nan ± inf i) = (+inf, nan). Building re+1 separately also
keeps a negative-zero imaginary part that the naive 1 + z would flip on the cut.
public static Complex Log1p(Complex z)
Parameters
zComplex
Returns
Reciprocal(Complex)
Complex reciprocal matching NumPy (nc_recip): Smith's algorithm specialised to a unit
numerator. This is overflow-safe (so 1/(huge) doesn't prematurely flush to 0) AND
reproduces NumPy's signed zeros (the imaginary part is -rat*scl), neither of which
operator /(double, Complex) gets right. Bit-identical to NumPy across finite, zero,
and infinite inputs.
public static Complex Reciprocal(Complex z)
Parameters
zComplex
Returns
Sin(Complex)
Complex sine matching NumPy. NumPy defines npy_csin(z) = -isinh(iz) verbatim
(i*z = (-y, x), then (w.Im, -w.Re)), so routing through the C99-correct
Sinh(Complex) reproduces NumPy bit-for-bit — including the large-imaginary case where
Sin(Complex) returns NaN from cosh(huge)*0.
public static Complex Sin(Complex z)
Parameters
zComplex
Returns
Sinh(Complex)
Complex hyperbolic sine matching NumPy. The y == 0 guard returns (sinh(x), y)
(so sinh(huge + I0) stays (±inf, ±0) instead of the BCL's cosh(x)sin(0) =
inf0 = NaN); the general finite case is the textbook sinh(x)cos(y) + i cosh(x)sin(y),
overflowing to ±inf exactly where NumPy's libm does. Non-finite inputs follow C99 Annex G.
public static Complex Sinh(Complex z)
Parameters
zComplex
Returns
Sqrt(Complex)
Complex square root matching NumPy (npy_csqrt, FreeBSD msun). Ported in full
(exact arithmetic, not transcendental) so the branch-cut signs and signed zeros that
Sqrt(Complex) drops are reproduced bit-for-bit.
public static Complex Sqrt(Complex z)
Parameters
zComplex
Returns
Square(Complex)
Complex square matching NumPy (np.square(z) == zz). NumPy's complex multiply is the
textbook (arbr - aibi, arbi + aibr), but compiled with FMA contraction; for
zz this is (fma(re, re, -(imim)), fma(re, im, imre)). The FMA path is what
produces NumPy's square(1e-10+1e-10i).real = -2.275e-37 (exact re² minus rounded im²)
and square(1e300+1e300i).real = -inf (not NaN) — both of which operator *(double, Complex)
(no FMA) turns into 0 and NaN respectively.
public static Complex Square(Complex z)
Parameters
zComplex
Returns
Tan(Complex)
Complex tangent matching NumPy. NumPy defines npy_ctan(z) = -itanh(iz) verbatim
(i*z = (-y, x), then (w.Im, -w.Re)), so routing through the Kahan-based
Tanh(Complex) reproduces NumPy bit-for-bit (the BCL Tan(Complex) drifts
tens of ULP near ±π/2).
public static Complex Tan(Complex z)
Parameters
zComplex
Returns
Tanh(Complex)
Complex hyperbolic tangent matching NumPy (full npy_ctanh port, FreeBSD msun). The
finite case uses Kahan's algorithm — t=tan(y); beta=1+t^2; s=sinh(x); rho=sqrt(1+s^2);
tanh = (betarhos + i t) / (1 + beta*s^2) — which is materially more accurate than
Tanh(Complex) (e.g. tan(1.5) drifts ~33 ULP through the BCL). The
|x| >= 22 branch avoids spurious overflow, and non-finite inputs follow C99.
public static Complex Tanh(Complex z)
Parameters
zComplex