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. Every non-finite input additionally matches NumPy's
exact NaN SIGN bit-for-bit (win-amd64 = MSVC UCRT complex functions): the "produce a NaN" slots
emit the positive NPY_NAN; csqrt/clog/cnc_log1p/cexp PROPAGATE an input NaN's sign;
csinh/ccosh CANONICALISE to +NaN (so csin/ccos follow the transform's negate),
while ctanh propagates; and a genuine 0/0 · inf-inf · inf*0 keeps its x86-negative
sign on both engines. Gated by the complex-unary NaN-contract oracle tier (raw-byte NaN compare).
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) = vfmaddsub z·z (matches NumPy's SIMD complex multiply
overflow/cancellation AND NaN sign); Reciprocal(Complex) = the CDOUBLE_reciprocal
ufunc loop (division-form imaginary term, NaN-sign correct);
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 FINITE inputs only, beyond 3 ULP — NaN sign is NOT
among them, it is byte-exact): 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 (and the sin/cos that route through them) 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 (npy_cabs) infinity/NaN semantics: a
±infinite real or imaginary part returns +inf regardless of the other part (including
NaN); a NaN component with no infinity returns the POSITIVE NaN (NumPy's npy_hypot
yields 0x7ff8…, where Abs(Complex) emits .NET's negative 0xfff8…).
All finite inputs defer to Abs(Complex) (bit-exact with NumPy).
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
Acosh(Complex)
Complex inverse hyperbolic cosine matching NumPy (npy_cacosh): a verbatim port of
the msun formula cacosh(z) = ±I*cacos(z) (sign chosen so Re >= 0) plus its
NaN/Inf special-value block, run on NumSharp's byte-exact Acos(Complex).
public static Complex Acosh(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
Asinh(Complex)
Complex inverse hyperbolic sine matching NumPy (npy_casinh). msun relation
casin(z) = Iconj(casinh(Iconj z)) inverts to casinh(z) = Iconj(casin(Iconj z)),
run on NumSharp's byte-exact Asin(Complex) — I*conj(z) = (z.Im, z.Re), a pure swap.
public static Complex Asinh(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
Atanh(Complex)
Complex inverse hyperbolic tangent matching NumPy (npy_catanh). The full msun port
already lives in Catanh(Complex) (it drives Atan(Complex)); this exposes it.
public static Complex Atanh(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. A FULLY finite input defers to Exp(Complex)
(ULP-identical to NumPy); any non-finite component follows ExpSpecial(double, double), which
reproduces MSVC UCRT cexp (NumPy 2.4.2 win-amd64) bit-for-bit — including the positive
NaN sign that Exp(Complex) emits negative for a finite real + non-finite imag.
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 — a verbatim port of the CDOUBLE_reciprocal ufunc
loop (numpy/_core/src/umath/loops.c.src), the Smith-style branch on the larger
component. This is overflow-safe (so 1/(huge) doesn't prematurely flush to 0) and
reproduces NumPy's signed zeros, which operator /(double, Complex) does not. Critically
the imaginary term is written as NumPy writes it — a DIVISION -1/d / -r/d, NOT
a negate of a reciprocal — so a NaN d keeps its (positive) sign instead of being
flipped by an explicit unary minus. Bit-identical to NumPy 2.4.2 across finite, zero,
infinite AND NaN-component 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 SIMD complex multiply
(simd_cmul in loops_arithm_fp.dispatch.c.src) is a fused multiply-add/subtract
(vfmaddsub): real = fused(rere - imim), imag = fused(reim + imre). The
fused subtract 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). It is spelled as
a TRUE MultiplySubtractScalar(Vector128<double>, Vector128<double>, Vector128<double>) rather than
FusedMultiplyAdd(re, re, -(imim)) so that a NaN imim keeps its (positive) sign
instead of being flipped by the explicit unary minus — matching NumPy 2.4.2 bit-for-bit on
NaN/inf-component inputs too. Off x86 (e.g. AArch64) there is no such intrinsic, so it uses
the portable FusedMultiplyAdd(double, double, double) — still a real fused
op, so the finite AND overflow results stay bit-identical to the x86 path (validated over 2M
random inputs plus the cancellation/overflow edges); only a NaN imim's sign differs
there, which is host-dependent off x86 and tokenized by the oracle (nanExact is x86-gated).
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