Class NDFloatMath
Transcendental helpers backing NumSharp's float unary ufuncs. Unlike the BCL's MathF/Math (which route to the platform libm and are only ~correctly rounded), each entry point here is a port of the kernel NumPy 2.4.2 actually runs, so the result is bit-identical to NumPy rather than merely within a couple of ULP.
Most of these are float32-only, because at float64 the platform libm already agrees with NumPy bit-for-bit on this host and there is nothing to port. Tanh(double) is the exception: NumPy ships its own kernel at both widths, so both diverged from the BCL and both are ported.
Exp(float) = simd_exp_FLOAT, Log(float) = simd_log_FLOAT
(numpy/_core/src/umath/loops_exponent_log.dispatch.c.src, the SIMD_AVX2_FMA3
instantiation). NumPy's own algorithm: clamp/flag the overflow and underflow ends, Cody-Waite
range reduction y = x - k·ln2 with k = rint(x·log2(e)), evaluate exp(y) as
the ratio of a 5th-order over a 2nd-order Remez minimax polynomial, then scale by 2^k.
It is NOT correctly rounded — NumPy documents a max error of 2.52 ULP (at
x = 0xc2781e37) — so reproducing NumPy means reproducing that error, operation for
operation, in the same order.
Why this is lane-width independent. Every step is elementwise (no cross-lane
reduction, shuffle or horizontal op), and the whole-vector if inside NumPy's
fma_scalef_ps is followed by a per-lane blend whose non-denormal lanes are unchanged by
the taken branch. This scalar entry point therefore yields the same bits the 8-lane
__m256 kernel yields — as do the Vector{128,256,512} overloads in
NDFloatMath.Simd.cs, which the IL kernels prefer when the host has FMA.
The FMA contraction is part of the answer, not an optimization. NumPy writes the
quadrant as _mm256_mul_ps(x, log2e) followed by _mm256_add_ps(quadrant, magic),
but the wheel's compiler (MSVC 19.44, the numpy==2.4.2 win-amd64 build) contracts that pair
into a single vfmadd. The difference is observable: at x = 0xc26d0e6c the
un-contracted product lands on the exact tie -85.5 and the magic-constant rint takes it
to -86 (half-to-even), while the fused form keeps the extra bits of
-85.4999987495703 and rounds to -85 — a 1-ULP difference in the result. Hence
FusedMultiplyAdd(float, float, float) here, matching the binary NumPy ships.
Specials (probed against 2.4.2, and a consequence of NumPy zeroing the NaN lanes
before the range comparisons, which makes the three masks mutually exclusive): any
NaN — quiet or signalling, either sign, any payload — returns the canonical
0x7fc00000; x ≥ 88.72283935546875 (incl. +inf) returns +inf;
x ≤ -103.97208404541015625 (incl. -inf) returns +0; ±0 returns
1. Results between those ends denormalize gracefully through the
ScalefDenormal(float, float) path. NumPy additionally raises the FP overflow/underflow status
flags (surfacing as a RuntimeWarning); NumSharp models no FP status word, so that
signalling is absent — a pre-existing, engine-wide difference, not one this port introduces.
Verified exhaustively: all 2^32 float32 bit patterns agree with NumPy 2.4.2 (chunked checksum sweep), not a sample. Perf: the finite path is AggressiveInlining so the JIT folds it into the IL-emitted unary kernel with no per-element call frame; the NaN/overflow/underflow ends and the denormal scale-back live in NoInlining cold helpers so the hot path stays inlineable.
public static class NDFloatMath- Inheritance
-
NDFloatMath
- Inherited Members
Methods
Cos(Vector128<float>)
NumPy 2.4.2's float32 cosine, 4 lanes at a time. See Cos(float).
public static Vector128<float> Cos(Vector128<float> x)Parameters
Returns
Cos(Vector256<float>)
NumPy 2.4.2's float32 cosine, 8 lanes at a time. See Cos(float).
public static Vector256<float> Cos(Vector256<float> x)Parameters
Returns
Cos(Vector512<float>)
NumPy 2.4.2's float32 cosine, 16 lanes at a time. See Cos(float).
public static Vector512<float> Cos(Vector512<float> x)Parameters
Returns
Cos(float)
NumPy 2.4.2's float32 cosine, bit-for-bit (simd_sincos_f32 with
SIMD_COMPUTE_COS - the same kernel, biasing the quadrant by one). See
Sin(float); note NumPy uses a SMALLER Cody-Waite limit for cosine.
public static float Cos(float x)Parameters
xfloat
Returns
Exp(Vector128<float>)
NumPy 2.4.2's float32 exponential, 4 lanes at a time. See Exp(float).
public static Vector128<float> Exp(Vector128<float> x)Parameters
Returns
Exp(Vector256<float>)
NumPy 2.4.2's float32 exponential, 8 lanes at a time. See Exp(float).
public static Vector256<float> Exp(Vector256<float> x)Parameters
Returns
Exp(Vector512<float>)
NumPy 2.4.2's float32 exponential, 16 lanes at a time. See Exp(float).
public static Vector512<float> Exp(Vector512<float> x)Parameters
Returns
Exp(float)
NumPy 2.4.2's float32 exponential, bit-for-bit (port of simd_exp_FLOAT).
public static float Exp(float x)Parameters
xfloat
Returns
Exp2(Vector128<float>)
float32 2^x, 4 lanes at a time. See Exp2(float).
public static Vector128<float> Exp2(Vector128<float> x)Parameters
Returns
Exp2(Vector256<float>)
float32 2^x, 8 lanes at a time. See Exp2(float).
public static Vector256<float> Exp2(Vector256<float> x)Parameters
Returns
Exp2(Vector512<float>)
float32 2^x, 16 lanes at a time. See Exp2(float).
public static Vector512<float> Exp2(Vector512<float> x)Parameters
Returns
Exp2(float)
float32 base-2 exponential (2^x), the fast replacement for
(float)Math.Pow(2, (double)x). Unlike Exp(float) this is not a
bit-for-bit port of a NumPy kernel — see the constants block for why (the numpy 2.4.2
win-amd64 wheel runs a SCALAR exp2f, its SVML vector kernel being AVX-512/Linux
only). It agrees with NumPy to ≤ 1 ULP (≈0.2% of inputs differ, all producing
subnormal outputs), the same accuracy class the Math.Pow path had, computed by
evaluating 2^r in double and rounding once.
Algorithm. Split x = n + r with n = rint(x) and
r ∈ [-0.5, 0.5], so 2^x = 2^r · 2^n. 2^r is a degree-8 minimax
polynomial evaluated in double (its 2^-37.7 approximation error is far below float
precision, so the narrowing to float is correctly rounded), then scaled by the integer
power 2^n straight into the exponent field — exact, so it introduces no second
rounding for a normal result. The n ≤ -125 denormal branch reuses the exp kernel's
split (ScalefDenormal(float, float)): clamp the field to 2^-125 and divide.
Lane-width independent, same as the exp/log ports: every step is elementwise,
so this scalar entry point yields the same bits as the Vector{128,256,512} overloads
in NDFloatMath.Simd.cs — verified 0-diff against the vector form over the full
~7M-input accuracy corpus.
Specials (probed against 2.4.2): any NaN returns the canonical
0x7fc00000; x ≥ 128 (incl. +inf) returns +inf; x ≤ -150
(incl. -inf) returns +0; ±0 returns 1.
public static float Exp2(float x)Parameters
xfloat
Returns
IsExpVectorAccelerated(int)
Whether Exp(Vector{ lowers to hardware
FMA on this machine. NumPy gates its own vector exp the same way (its kernel is compiled
only for AVX2+FMA3 / AVX-512F targets, and hosts without them run the scalar libm).vectorBits}<float>)
public static bool IsExpVectorAccelerated(int vectorBits)Parameters
vectorBitsint
Returns
IsTanhVectorAccelerated(int)
Whether Tanh(Vector{ can run its table
lookup in hardware on this machine. Distinct from IsExpVectorAccelerated(int)
because tanh needs a capability the pure-arithmetic kernels do not: its coefficients are
chosen PER LANE by the exponent of vectorBits}<float>)|x|, so they cannot be held in registers and the
vector form must read the table — AVX2, on top of the FMA they all need. A host without it
takes the scalar entry point, which computes the SAME bits one lane at a time; this gate
is about speed, never results.
The 128 arm is NOT dead code — do not "simplify" it away. It is tempting to
argue that the emitter picks 128 only when Vector256.IsHardwareAccelerated is
false, that this means AVX2 is absent, and therefore that 128 => Avx2.IsSupported
can never be true. That is WRONG, and the counter-example is a supported runtime knob
rather than exotic hardware: with DOTNET_PreferredVectorBitWidth=128 an ordinary
AVX2 host reports VectorBits = 128 while Avx2.IsSupported stays true
(probed). Dropping the arm demotes tanh to the scalar loop for everyone who caps the
width — a pure perf regression, invisible to any correctness gate because the bits are
identical either way. A host that genuinely lacks AVX2 is still handled: the arm is false
there.
public static bool IsTanhVectorAccelerated(int vectorBits)Parameters
vectorBitsint
Returns
Log(Vector128<float>)
NumPy 2.4.2's float32 natural logarithm, 4 lanes at a time. See Log(float).
public static Vector128<float> Log(Vector128<float> x)Parameters
Returns
Log(Vector256<float>)
NumPy 2.4.2's float32 natural logarithm, 8 lanes at a time. See Log(float).
public static Vector256<float> Log(Vector256<float> x)Parameters
Returns
Log(Vector512<float>)
NumPy 2.4.2's float32 natural logarithm, 16 lanes at a time. See Log(float).
public static Vector512<float> Log(Vector512<float> x)Parameters
Returns
Log(float)
NumPy 2.4.2's float32 natural logarithm, bit-for-bit (port of simd_log_FLOAT,
the sibling of Exp(float) in the same NumPy source file). Splits x into
exponent and mantissa, folds the mantissa into [1/sqrt 2, sqrt 2), evaluates
log(1+m) as a 5th-over-5th-order Remez minimax ratio, and adds exponent*ln2
with a final fused multiply-add. NumPy documents a max error of 3.83 ULP here (at
x = 0x3f486945) — reproducing NumPy means reproducing that.
public static float Log(float x)Parameters
xfloat
Returns
Sin(Vector128<float>)
NumPy 2.4.2's float32 sine, 4 lanes at a time. See Sin(float).
public static Vector128<float> Sin(Vector128<float> x)Parameters
Returns
Sin(Vector256<float>)
NumPy 2.4.2's float32 sine, 8 lanes at a time. See Sin(float).
public static Vector256<float> Sin(Vector256<float> x)Parameters
Returns
Sin(Vector512<float>)
NumPy 2.4.2's float32 sine, 16 lanes at a time. See Sin(float).
public static Vector512<float> Sin(Vector512<float> x)Parameters
Returns
Sin(float)
NumPy 2.4.2's float32 sine, bit-for-bit (port of simd_sincos_f32 with
SIMD_COMPUTE_SIN). Cody-Waite reduction against pi/2 in three fused steps, then a
sine/cosine polynomial pair selected by the quadrant.
Large arguments deliberately fall back to the platform libm, exactly as
NumPy does: beyond |x| = 117435.992 (sine) its own comment says Cody-Waite "becomes
inaccurate and we will call libc". Those inputs therefore inherit whatever difference
exists between .NET's Sin(float) and the CRT's sinf - the port cannot
close that, and extending the polynomial past NumPy's own limit would move us AWAY from
NumPy rather than toward it.
public static float Sin(float x)Parameters
xfloat
Returns
Tanh(double)
simd_tanh_f64 — the double-precision half of the same kernel: 16 subintervals and a
degree-16 Horner evaluation. See Tanh(float) for the algorithm, the
lane-width argument and the specials; the only structural difference is the index
extraction. NumPy computes it with a 32-bit max/min over the 64-bit lanes, noting it
"is fine … since we're not crossing 32-bit edge": the mask and the bias both have a zero
low half, so the subtraction cannot borrow and the whole index lives in the high 32 bits.
This port clamps that high half directly, which is the same arithmetic without the
bitcast.
That re-spelling is the one line here a reader cannot diff against NumPy's, so
it is PROVEN rather than argued — and proven exhaustively, not sampled.
ndnan = bits & 0x7ff8000000000000 keeps only 11 exponent bits plus mantissa bit
51 and clears everything else, and both formulas depend on nothing else, so the input
domain is exactly 2^12 = 4096 values. TanhFloat64IndexExtraction_MatchesNumPySpelling
(Math/np.Tanh.IndexExtraction.Test.cs) enumerates all 4096 and asserts NumPy's spelling and
this one agree on every single one — and separately confirms that real double bit patterns
actually reach all 4096, so the domain is exact rather than an over-approximation.
Verified against NumPy 2.4.2 over a 300K-value structured sweep (specials, subnormals, both tails, interior) and the committed corpus tier; float64 has no exhaustive analogue.
public static double Tanh(double x)Parameters
xdouble
Returns
Tanh(Vector128<float>)
NumPy 2.4.2's float32 tanh, 4 lanes at a time. See Tanh(float).
This width gathers where its 256-bit sibling transposes, and that is deliberate —
the verdict flips with the lane count. A subinterval is 8 consecutive floats, so the
transpose reads whole rows and is unbeatable when 8 lanes each need one (see
Tanh(Vector256<float>)). At 4 lanes it is not: expressing this width as the
low half of an 8-lane call — which would give the file a single table read — computes four
lanes it throws away, and measured 0.171 ms / 17.56 ms against the gather's
0.139 ms / 13.58 ms at 100K / 10M under DOTNET_PreferredVectorBitWidth=128,
i.e. the gather is 1.23-1.29x faster and the composed form is barely ahead of the scalar
loop (1.06x) where the gather is 1.31-1.39x. Two lookup strategies is the right answer
here; both are verified bit-exact against the same table.
public static Vector128<float> Tanh(Vector128<float> x)Parameters
Returns
Tanh(Vector256<float>)
NumPy 2.4.2's float32 tanh, 8 lanes at a time. See Tanh(float).
The coefficient fetch is a transpose, not a gather — the one place this port
deliberately departs from the shape of NumPy's own AVX2 code, and it is worth the words.
A subinterval's 8 coefficients are 8 consecutive floats, i.e. exactly one
Vector256, so the whole table read is 8 contiguous row loads followed by an 8x8
transpose that lands each coefficient across the lanes — where a gather-per-coefficient
would issue 8 vgatherdps, 64 separate element loads for the same 8 rows. Both were
implemented and measured here, bit-identical and against the same table: transpose
0.071 ms vs gather 0.124 ms at 100K, and 7.41 ms vs 12.87 ms at 10M — 1.7x,
and the difference between losing to NumPy and beating it. (NumPy reaches the same place
differently, permuting a pre-transposed table; its layout is chosen to suit AVX-512's
two-table lookup, which this host does not have.)
public static Vector256<float> Tanh(Vector256<float> x)Parameters
Returns
Tanh(Vector512<float>)
NumPy 2.4.2's float32 tanh, 16 lanes at a time — composed from two 8-lane halves rather than written against AVX-512 gathers. Every lane is independent here, so splitting the register is an identity, and it means the 512-bit width runs exactly the code path that was verified exhaustively rather than an untested transcription of it. The AVX-512 gather would save one instruction per coefficient and is not worth an unverifiable kernel. See Tanh(float).
public static Vector512<float> Tanh(Vector512<float> x)Parameters
Returns
Tanh(float)
simd_tanh_f32 — NumPy's table-driven tanh, converted by NumPy from Intel SVML's
svml_z0_tanh_s_la. tanh is odd, so the kernel works on |x| and re-applies the
sign bit at the end. [0, SATURATION_THRESHOLD) is split into 32 subintervals; the
exponent field of |x| indexes straight into a table of per-subinterval minimax
coefficients and a recentring offset b, and the result is a degree-6 Horner
evaluation of P(|x| - b). The top subintervals hold 1 + 0·y + 0·y² …, which is
how the saturated range returns exactly 1 through the same arithmetic.
Why this is lane-width independent. Every step is elementwise; the only
width-dependent code in NumPy is how the table is fetched (a 4×4 lane transpose on
SSE4, TwoTablesLookupLanes on AVX2/AVX-512, a gather elsewhere), and all three fetch
the same coefficients. So this scalar entry point, the
Vector{128,256,512} overloads in NDFloatMath.Simd.cs, and NumPy's own kernel
at any width all produce identical bits.
Unlike the exp/log/sincos ports, the FMA here is not a host pin. Those had to
reproduce MSVC's contraction of a separate multiply and add; this kernel spells its Horner
steps as hn::MulAdd — an explicit fused multiply-add — so
FusedMultiplyAdd(float, float, float) is a literal transcription rather than a guess about
the compiler.
Specials (probed against 2.4.2): any NaN — quiet or signalling, either sign,
any payload — returns the canonical 0x7fc00000; ±inf and every |x| past
the saturation threshold return ±1; ±0 returns ±0 (subinterval 0 has
b = c0 = 0 and c1 = 1, so the polynomial collapses to |x| and the sign
bit is OR-ed back). NumPy also clears the FP status word around this loop; NumSharp models
no FP status, a pre-existing engine-wide difference this port does not change.
Verified exhaustively: all 2^32 float32 bit patterns agree with NumPy 2.4.2, through both this entry point and the SIMD kernel.
public static float Tanh(float x)Parameters
xfloat