SharpLibm

The C math library (libm) in C#: sin, pow, fmod, floorf and the rest, with C names, C signatures and C99/C23 semantics. Written to compile under NativeAOT with NoStdLib — no CoreLib math, no native libm, no static constructors, no heap — and usable from any other .NET host as well.

Status: early.

• Exact operations (rounding to integral values, fmod/remainder/remquo, frexp/ldexp/scalbn/ilogb/logb/nextafter, copysign/fabs/fdim/ fmax/fmin, llrint/llround), double and float: bit-exact against musl's libc-test vectors.
• Transcendental functions — 47 for double, 48 for float, from CORE-MATH through CoreMathSharp (correctly rounded): exp/log families, pow, trig and inverse trig, the π-scaled forms, hyperbolic, erf/erfc, tgamma/lgamma, hypot, cbrt, rsqrt, compoundf, and a correctly rounded software fma. Double: bit-exact against MPFR on all of CORE-MATH's worst cases (≈30 M inputs). Float: bit-exact against CORE-MATH's C on every one of the 2^32 inputs.
• CoreMathSharp 1.0.0 is not correctly rounded on hard inputs: about 80 translation errors, several of them gross (pow, tgamma, log1p, expm1, rsqrtf, exp10f). They are fixed here and marked // SharpLibm fix:, together with CORE-MATH fixes made after the port (hypot, cospi, pow, error bounds).
• Bessel functions j0, j1, jn, y0, y1, yn and the inverse error functions erfinv, erfcinv (not ISO C; CUDA's and Go's names), double and float, ported from Go's math (fdlibm's algorithms for the Bessel functions, Wichura's AS 241 plus one Newton step on erf/erfc for erfinv). Not correctly rounded: measured against MPFR, erfinv within 1.33 ulp, erfcinv within 4.6 (below 1e-300, where the step is skipped, otherwise ≈1); j0/j1/y0/y1 within 4 ulp of the function's amplitude (absolute error — near a zero the relative error is unbounded, as in every fdlibm-derived libm), jn/yn within ≈33 for orders up to 50. The float forms go through double and were correctly rounded on every input tested.
• C23 compoundn(x, n) / compoundnf(x, n) = (1 + x)^n for an integer n, with IEEE 754-2019's special cases. compoundnf is correctly rounded where n is exact in float (CORE-MATH's compoundf underneath); compoundn is within one ulp (glibc's vectors). compoundf(x, y) with a real y is CORE-MATH's extension, correctly rounded.
• Builds and links under NativeAOT 8 NoStdLib with no native math left.

Functions

Every function, double and float on one row. .NET: in SharpLibm.Libm; Windows / Linux: exported by the C library (native/). Generated by SharpLibm-tests (tools/readme_table.py) from the builds; accuracy is what its runs measure (round-to-nearest).

function	semantics	.NET	Windows	Linux	accuracy
floor, floorf	largest integer ≤ x	✅	✅	✅	exact
ceil, ceilf	smallest integer ≥ x	✅	✅	✅	exact
trunc, truncf	integer part, towards zero	✅	✅	✅	exact
round, roundf	nearest integer, halves away from zero	✅	✅	✅	exact
roundeven, roundevenf	nearest integer, halves to even	✅	✅	✅	exact
rint, rintf	nearest integer in the rounding mode (to even)	✅	✅	✅	exact
nearbyint, nearbyintf	as rint, raising no inexact	✅	✅	✅	exact
modf, modff	fraction, integer part through a pointer	✅	✅	✅	exact
lrint, lrintf	rint as C long (32-bit on Windows)	✅	✅	✅	exact
lround, lroundf	round as C long (32-bit on Windows)	✅	✅	✅	exact
llrint, llrintf	rint as long long	✅	✅	✅	exact
llround, llroundf	round as long long	✅	✅	✅	exact
fmod, fmodf	x − n·y, n truncated	✅	✅	✅	exact
remainder, remainderf	x − n·y, n to nearest even	✅	✅	✅	exact
remquo, remquof	remainder, low quotient bits through a pointer	✅	✅	✅	exact
frexp, frexpf	mantissa in [½, 1), exponent through a pointer	✅	✅	✅	exact
ldexp, ldexpf	x·2ⁿ	✅	✅	✅	CR
scalbn, scalbnf	x·2ⁿ, n an int	✅	✅	✅	CR
scalbln, scalblnf	x·2ⁿ, n a C long	✅	✅	✅	CR
ilogb, ilogbf	binary exponent as int	✅	✅	✅	exact
logb, logbf	binary exponent as floating point	✅	✅	✅	exact
nextafter, nextafterf	next representable value towards y	✅	✅	✅	exact
nexttoward, nexttowardf	nextafter with a long double y ⁵	✅	❌	❌	exact
copysign, copysignf	|x| with the sign of y	✅	✅	✅	exact
fabs, fabsf	|x|	✅	✅	✅	exact
fdim, fdimf	max(x − y, 0)	✅	✅	✅	CR
fmax, fmaxf	maximum, ignoring a NaN operand	✅	✅	✅	exact
fmin, fminf	minimum, ignoring a NaN operand	✅	✅	✅	exact
fma, fmaf	x·y + z rounded once	✅	✅	✅	CR
sqrt, sqrtf	√x (one instruction)	✅	✅	✅	CR
rsqrt, rsqrtf	1/√x	✅	✅	✅	CR
cbrt, cbrtf	∛x	✅	✅	✅	CR
hypot, hypotf	√(x² + y²) without overflow	✅	✅	✅	CR
exp, expf	eˣ	✅	✅	✅	CR
exp2, exp2f	2ˣ	✅	✅	✅	CR
exp10, exp10f	10ˣ	✅	✅	✅	CR
expm1, expm1f	eˣ − 1	✅	✅	✅	CR
exp2m1, exp2m1f	2ˣ − 1	✅	✅	✅	CR
exp10m1, exp10m1f	10ˣ − 1	✅	✅	✅	CR
log, logf	ln x	✅	✅	✅	CR
log2, log2f	log₂ x	✅	✅	✅	CR
log10, log10f	log₁₀ x	✅	✅	✅	CR
log1p, log1pf	ln(1 + x)	✅	✅	✅	CR
log2p1, log2p1f	log₂(1 + x)	✅	✅	✅	CR
log10p1, log10p1f	log₁₀(1 + x)	✅	✅	✅	CR
pow, powf	xʸ	✅	✅	✅	CR
compoundn, compoundnf	(1 + x)ⁿ, n a long long (C23)	✅	✅	✅	≤ 1 ulp / CR ³
compoundf	(1 + x)ʸ, y real (CORE-MATH extension)	✅	✅	✅	CR
sin, sinf	sin x	✅	✅	✅	CR
cos, cosf	cos x	✅	✅	✅	CR
tan, tanf	tan x	✅	✅	✅	CR
sincos, sincosf	sin x and cos x through pointers	✅	✅	✅	CR
asin, asinf	arcsin x	✅	✅	✅	CR
acos, acosf	arccos x	✅	✅	✅	CR
atan, atanf	arctan x	✅	✅	✅	CR
atan2, atan2f	arctan(y/x) in the quadrant of (x, y)	✅	✅	✅	CR
sinpi, sinpif	sin πx	✅	✅	✅	CR
cospi, cospif	cos πx	✅	✅	✅	CR
tanpi, tanpif	tan πx	✅	✅	✅	CR
asinpi, asinpif	arcsin(x)/π	✅	✅	✅	CR
acospi, acospif	arccos(x)/π	✅	✅	✅	CR
atanpi, atanpif	arctan(x)/π	✅	✅	✅	CR
atan2pi, atan2pif	atan2(y, x)/π	✅	✅	✅	CR
sinh, sinhf	sinh x	✅	✅	✅	CR
cosh, coshf	cosh x	✅	✅	✅	CR
tanh, tanhf	tanh x	✅	✅	✅	CR
asinh, asinhf	arsinh x	✅	✅	✅	CR
acosh, acoshf	arcosh x	✅	✅	✅	CR
atanh, atanhf	artanh x	✅	✅	✅	CR
erf, erff	error function	✅	✅	✅	CR
erfc, erfcf	1 − erf x	✅	✅	✅	CR
tgamma, tgammaf	Γ(x)	✅	✅	✅	CR
lgamma, lgammaf	ln |Γ(x)|	✅	✅	✅	CR
lgamma_r, lgammaf_r	ln |Γ(x)|, sign of Γ through a pointer	✅	✅	✅	CR
erfinv, erfinvf	inverse of erf (not ISO C)	✅	✅	✅	≤ 1.33 ulp / CR ²
erfcinv, erfcinvf	inverse of erfc (not ISO C)	✅	✅	✅	≤ 4.6 ulp ⁴ / CR ²
j0, j0f	Bessel J₀	✅	✅	✅	≤ 4 ulp ¹ / CR ²
j1, j1f	Bessel J₁	✅	✅	✅	≤ 4 ulp ¹ / CR ²
jn, jnf	Bessel Jₙ, n an int	✅	✅	✅	≤ 33 ulp ¹ / CR ²
y0, y0f	Bessel Y₀	✅	✅	✅	≤ 4 ulp ¹ / CR ²
y1, y1f	Bessel Y₁	✅	✅	✅	≤ 4 ulp ¹ / CR ²
yn, ynf	Bessel Yₙ, n an int	✅	✅	✅	≤ 27 ulp ¹ / CR ²

exact — the result is representable, bit for bit. CR — correctly rounded: double against MPFR on CORE-MATH's worst cases (≈30 M inputs), float against CORE-MATH's C on all 2³² inputs (2²⁸ random pairs for two arguments); fma, sqrt, remainder and roundings also against TestFloat; ldexp/scalbn are exact unless the result is subnormal. Two values: double / float.

¹ Near a zero of the function the error is absolute, in ulp of its local amplitude; relative error there is unbounded, as in every fdlibm-derived libm. ² On every input tested (≈40–50 k against MPFR), not proven. ³ Where n is exact in float (|n| ≤ 2²⁴, and larger n with trailing zeros); otherwise through double, rounded once more. ⁴ ≈1 ulp above 1e-300; below, the Newton step is skipped. ⁵ Not exported to C: its long double is double with MSVC and 80-bit with gcc and mingw, so no one ABI fits; use nextafter.

Design

• Base bodies use SSE2 only. The library does not detect the CPU and does not depend on it beyond x86-64's baseline.
• Replaceable entry points. Functions that one instruction does on newer CPUs — floor, ceil, trunc, rint, nearbyint (SSE4.1 roundsd/roundss) and fma, fmaf (FMA3), with their float forms — are kept out of line ([MethodImpl(NoInlining)]). A host that knows the CPU may overwrite their machine code with the instruction; every caller, inside the library as well, goes through the entry point and picks that up (modf calls trunc, so it speeds up with it). Without that, the library works the same, only slower. SharpLibm itself does none of this.
• Semantics are C's, not System.Math's: fmax ignores a NaN operand, modf of ±inf gives a ±0 fraction, NaN signs are IEEE signs.
• The rounding mode is always round-to-nearest (the CLI fixes it) and no floating-point exception flags are raised.
• Future work, for C callers through the exports (<fenv.h>): the other three rounding modes (CORE-MATH has the code paths, CoreMathSharp dropped them), exception flags without spurious ones, and errno. .NET code cannot observe any of these, so System.Math does not need them.

Using it

Import SharpLibm.props into the consuming project; the sources compile into that assembly (in a NativeAOT NoStdLib image the application is its own system module, so a separate compiled assembly would have no CoreLib to bind against).

<Import Project="path/to/SharpLibm/SharpLibm.props" />

Then call SharpLibm.Libm.fmod(x, y) and so on.

Properties:

• SharpLibmExports=true — export every function under its C name ([RuntimeExport("fmod")]) so C code linked into the image can call it. Off by default: an exported method is a root, and ILC would then keep every function and its tables in the image.
• SharpLibmDefineRuntimeExport=true — compile a stand-in System.Runtime.RuntimeExportAttribute for hosts whose CoreLib has none (an ordinary .NET project). A NoStdLib host declares its own.
• SharpLibmNoExceptions=true — CORE-MATH's three "unexpected worst case" guards (pow, atan2, compoundf) end the process instead of throwing InvalidOperationException, for a host with no exception types. The C original prints the input and calls exit(1) there; none of the worst-case tests reaches them.

The three ways it is used:

• An ordinary .NET assembly — import the props (with SharpLibmDefineRuntimeExport=true); SharpLibm-tests is built this way.

• Inside a NativeAOT NoStdLib system — SharpOS imports the props into its std; System.Math forwards to Libm.

• A C library — sharplibm.dll / sharplibm.so (≈700 KB): every function exported under its C name, no imports at all (no libc, no CRT, no OS calls), no .NET runtime, no GC, no allocation. Its core library is native/MinimalRuntime.cs plus std-no-runtime-sharplibm/ — cut copies of SharpOS's apps_native/sdk/MinimalRuntime.cs and std-no-runtime/, same paths, each file says what it keeps. Errors that C would report by crashing (bounds, the guards above) end the process like abort().

  dotnet publish native -c Release
  

  builds it for the host: the DLL on Windows (MSVC link.exe — Visual Studio's C++ workload, vswhere on PATH, as in a Developer Command Prompt), the .so on Linux (the system's clang, as the .NET SDK uses it). Next to it goes sharplibm.h, generated from the same [RuntimeExport]s as the export table, so the two always agree; it compiles alongside <math.h>.

  C ABI details. C's long is 32 bits on Windows and 64 on Linux: lrint, lround, lrintf, lroundf, scalbln, scalblnf are exported at the target's width (CAbiExports.cs). nexttoward/nexttowardf are not exported — their long double is double with MSVC and 80-bit with gcc/mingw; use nextafter.

Host contract

What the consuming CoreLib must provide:

• unsafe code (pointers), MethodImplAttribute;
• System.Math.Sqrt marked [Intrinsic], so the compiler emits sqrtsd (an ordinary .NET CoreLib has this);
• Span<T> / ReadOnlySpan<T> with Slice, and RuntimeHelpers.CreateSpan<T> — constant tables are ReadOnlySpan<T> t = [...], which then live in the image's data instead of being allocated;
• MemoryMarshal.Cast<TFrom, TTo>(ReadOnlySpan<TFrom>), stackalloc;
• ValueTuple<…> (2 to 5 elements) and deconstruction;
• the constants of the primitives (double.NaN, int.MaxValue, …);
• InvalidOperationException, unless SharpLibmNoExceptions=true.

No System.Range, no records, no BitConverter, no double.IsNaN and the like (classification is on the bits, Bits.cs), no interfaces, no delegates, no string formatting, no static constructors, no heap allocation. Builds with C# 12.

Tests

The tests live in a separate repository, SharpLibm-tests: the test suites they use come under many licences (MIT, BSD, Sun, GPL, LGPL), so they are kept apart from this MIT code. They run on an ordinary .NET 10 host and check, bit for bit, in round-to-nearest:

• exact operations (floor … fmod, remquo, frexp, nextafter …) against musl's libc-test;
• every transcendental function within 0.5 ulp on libc-test's vectors;
• the double functions on CORE-MATH's worst cases (≈30 million inputs) against MPFR;
• the float functions on all 2^32 inputs against CORE-MATH's C;
• the Bessel functions and erfinv/erfcinv, error in ulp against MPFR (≈1.9 M inputs), within per-function bounds, and their special values.

Licence

MIT (LICENSE). Code derived from other projects keeps their notices; see THIRD-PARTY-NOTICES.md and the header of each source file.