Native special functions and constants

<uni20/core/math.hpp> exposes math::beta, math::zeta, math::expint, math::upper_gamma, math::airy_ai, math::lgamma, math::lgamma_sign, math::bessel_j, and math::bessel_y for native real types when a suitable provider or Uni20 implementation is available. These paths add no library dependency and do not use Boost or MPFR to evaluate native arguments.

#include <uni20/core/math.hpp>
 
// On a platform with the double-precision providers:
auto gamma = uni20::math::lgamma_sign(-0.5);
// gamma.value = log(abs(Gamma(-0.5))), gamma.sign = -1
auto logarithm = uni20::math::lgamma(-0.5); // Same value, without the sign.
auto j = uni20::math::bessel_j(3, 1.0);     // J_3(1)
auto y = uni20::math::bessel_y(-3, 1.0);    // Y_-3(1) = -Y_3(1)

Standard-library special functions

When <cmath> advertises __cpp_lib_math_special_functions >= 201603L, Uni20 selects this adapter for float, double and long double:

Uni20 callable	Standard provider	Domain and interpretation
math::beta(x,y)	std::beta(x,y)	Both arguments positive; same native real type

auto b = uni20::math::beta(0.5L, 0.5L); // long double, approximately pi

The standard feature macro is sufficient to select this interface; there are no redundant per-function CMake probes. The return type is the input type. Integer arguments, mixed native types and conversion-only classes are not accepted by these adapters. Convert explicitly when needed.

The standard library determines Beta's numerical accuracy, floating-point error reporting and range behavior. Native Beta requires positive arguments and does not promise MPFR's extended real domain or correct rounding. The MPFR overload retains its own domains and precision rules.

If the macro is absent, Beta uses Uni20's native fallback, which also supplies configured binary128. Native Ei and zeta use Uni20 algorithms at every configured precision. Their former standard paths were replaced after accuracy and performance comparisons; neither requires the standard special-function feature macro. A standard template declaration accepting an extension type alone does not establish precision-preserving evaluation.

Shared numerical tests cover all configured native formats and MPFR128/256, with independent identities, wider references and measured accuracy gains. Identity tests also run without MPFR.

Global Gamma sign state

Application code using these adapters must not read, write, or rely on the global signgam variable, including reading it after a non-reentrant platform lgamma call. Use math::lgamma_sign(x) when the Gamma sign is needed; its returned sign belongs to that invocation.

This restriction accommodates a known upstream limitation. The inspected libstdc++ implementation of native Beta calls ordinary lgamma, which writes global signgam on glibc. Its numerical results do not consume that sign output. Uni20 retains the standard-library Beta adapters, but concurrent calls on affected providers still have an internal write/write data race even when application code never accesses signgam. The restriction does not make those providers formally race-free, and race detectors can report the upstream defect. Uni20 adds no serialization around the retained Beta providers. The same underlying issue is tracked for std::poisson_distribution in GCC bug 111726.

This limitation concerns the affected native Beta providers. Uni20's Beta fallback is reentrant and does not access signgam; its presence does not remove the limitation from the retained standard providers. Native zeta now uses reentrant Uni20 evaluation without signgam. Uni20's own log-Gamma adapters use the reentrant routines described below.

Pi-scaled trigonometric functions

math::sinpi, cospi, tanpi, asinpi, acospi, atanpi and atan2pi support native float32/64/80 and configured binary128, as well as the existing MPFR scalar path. The forward functions evaluate at pi times the argument; the inverse functions return angles divided by pi. Native unary results retain the argument type. atan2pi(y,x) uses the common native floating type, subject to the usual facade check that it preserves both operands' precision and exponent range. Integer and conversion-only inputs are not accepted. As with the existing elementary math interface, a trailing Precision is an MPFR operation, not a native overload.

auto phase = uni20::math::sinpi(0.5L); // exactly 1, without forming pi/2
auto turns = uni20::math::atan2pi(1.0, -1.0); // exactly 0.75

CMake checks the linked C23 forward family (sinpi, cospi, tanpi) and inverse family (asinpi, acospi, atanpi, atan2pi) separately for each native type. It requires exact signatures, including _Float128 signatures for binary128. The capability macros are UNI20_HAS_PI_TRIG_{FLOAT,DOUBLE,LONG_DOUBLE,FLOAT128} and UNI20_HAS_PI_INVERSE_{FLOAT,DOUBLE,LONG_DOUBLE,FLOAT128}. A complete detected family uses its libc provider. A missing or incomplete family uses a portable implementation at the original precision. No newer libc is required, and MPFR is never used to evaluate native arguments.

The portable forward implementation uses remquo(x, 0.5) to recover the quadrant and an exact binary remainder in [-1/4,1/4]. It multiplies only that remainder by the typed pi constant before calling native sin, cos or tan. This preserves phase for large arguments, avoids overflow in 2*x or pi*x, and retains relative accuracy near integer zeros and half-integer poles. The inverse implementation scales the native inverse function, with explicit handling of exact endpoints, axes and diagonals. atan2pi uses atan2, avoiding an overflowing or underflowing intermediate y/x.

The native value conventions follow IEEE pi-scaled arithmetic:

• At every integer, sinpi and tanpi return zero with the input's sign.
• cospi returns positive zero at every half-integer.
• tanpi(n+1/2) returns positive infinity for even n and negative infinity for odd n; at quarter-integers its finite values are exactly +1 or -1.
• Forward functions return NaN for infinite or NaN inputs. asinpi and acospi have domain [-1,1], with NaN outside it.
• atanpi maps infinities to +/-1/2. atan2pi preserves signed-zero quadrants, including both zero arguments, and handles infinite axes and diagonals. NaN in either argument yields NaN.

The MPFR path retains its provider conventions. In particular, MPFR 4.2.1 uses an alternating zero sign for tanpi at odd integers; native results follow the input sign, as described in the current MPFR reference.

See the MPFR function reference and the GNU libc forward and inverse interfaces. Native implementations do not promise correct rounding or a portable errno/floating-point-exception contract. Their accuracy depends on the native elementary provider and rounding of the typed pi constant. Tests use round-to-nearest with gradual underflow, require errors within four output grid spacings on a deterministic adversarial sample, and separately check improvement with precision using analytic identities and a convergent series. This is numerical validation, not a proven global error bound.

The shared precision matrix tests both the selected public path and the portable implementation explicitly. Thus a successful fallback does not count as validation of an unavailable libc provider. Core boundary tests run without MPFR; the wider 512-bit oracle checks are reported unavailable in that build.

Numerical acceptance

The following acceptance criteria apply to configured native float32/64/80/128 under round-to-nearest with gradual underflow. The reference is the function at the exact represented input, not at an intended real input before rounding. PiScaledAccuracy and PiScaledPortableAccuracy cover the public and portable paths separately in tests/numerics/test_pi_math_precision.cpp.

Region	Measure and intended criterion	Guarantee or evidence
Integer/half-integer grids; quarter-integer tangent; inverse endpoints, axes and equal-magnitude diagonals	Exact value and specified zero/infinity sign	Value contract; core tests plus precision-matrix identities
Ordinary non-dyadic sample (i/37, -128 <= i <= 128); both neighbors of quarter integers from -4 to 4	Absolute error <= four output grid spacings against 512-bit MPFR	Empirical acceptance target, not a domain-wide bound
Large arguments with representable quarter offsets; native minimum normal/subnormal multiples and largest finite values	Same absolute grid-spacing target, with exact checks for zeros and nonfinite results	Tests phase retention and exponent range separately from significand precision
Selected nonzero analytic values at 1/8, 1/2 and 1/4	Four-spacing target and error at least sixteen times smaller than the preceding precision's rounded reference	Precision-improvement evidence from radical/rational identities and an independent atan series
Out-of-domain and nonfinite arguments	NaN classification, or the specified finite/infinite value and sign	Value contract, independently checked in core tests

For a finite nonzero result r, an output grid spacing is next_up(abs(r)) - abs(r), evaluated in the wider reference type. This also gives an absolute measure in the subnormal region. No relative error is used at zeros; exact-zero signs have separate assertions. Pole-neighbor tests assess evaluation at the represented input, not the conditioning of an unrounded mathematical argument.

The recorded GCC 13 and Clang 19 runs satisfy these thresholds for their configured portable paths. The following measurements use the existing PiScaledAccuracy/PiScaledPortableAccuracy deterministic grid, with GCC 13.3, glibc 2.39 and MPFR 4.2.1 at 512 bits. Public and direct portable paths have identical measurements on this host; all eight typed libc pi-family probes are unavailable. These sample maxima are not domain-wide guarantees. The intended bound remains four output-grid spacings.

Function	float32	float64	binary80	binary128
sinpi	0.983119	0.920274	0.839963	0.810794
cospi	0.997654	0.860209	0.974611	0.927585
tanpi	1.668568	1.867259	1.452005	1.960705
asinpi	1.304351	1.327294	0.897331	1.203063
acospi	1.256780	1.238159	1.318098	1.255125
atanpi	1.562144	1.401984	1.323161	1.233538
atan2pi	1.716090	1.401984	1.330636	1.832722

Maxima are rounded upward to six decimal places. Each format checks 402/416/ 378/112/112/438/1756 finite nonzero references for the seven functions in table order. Zeros and nonfinite values have separate classification/sign checks and are excluded from these maxima. Worst inputs are below; decimal values retain enough digits to round back to the tested native value. atan2pi entries give (y,x) in call order.

Function	float32 input	float64 input	binary80 input	binary128 input
sinpi	-0.243243247	-3.1891891891891890	-3.10810810810810810804	-3.13513513513513513513513513513513510
cospi	-1.54054058	-1.5945945945945945	-1.51351351351351351350	-0.594594594594594594594594594594594592
tanpi	-3.35135126	-3.7500000000000004	-3.40540540540540540536	-3.45945945945945945945945945945945958
asinpi	-0.999999940	-0.62162162162162160	-0.108108108108108108106	-0.513513513513513513513513513513513472
acospi	-0.837837815	0.43243243243243246	0.75	0.0810810810810810810810810810810810840
atanpi	-2.37837839	67108863.25	-0.837837837837837837835	-2.32432432432432432432432432432432448
atan2pi	(-1.89189184, -1)	(67108863.25, 1)	(-0.249999999999999999986, -1)	(-2.86486486486486486486486486486486490, 3)

At all 28 worst points, the 1024-bit oracle rounds back to the same 512-bit reference. Measurement adds maximum/input bookkeeping to the existing sample checker without changing its inputs or acceptance assertions. The permanent tests continue to enforce the four-spacing target and independent precision improvement checks. A provider unavailable on a test host has no numerical evidence from that run, even if its wiring is checked synthetically. Accuracy on a new libc provider must be measured on that provider before calling it validated.

Compound powers

math::compound(x,n) supports native float32/64/80/128 and MPFR. It computes (1+x)^n from the original input, including increments lost by first rounding 1+x. The integer exponent has the shared checked signed64 domain; the real input must be at least -1. Native paths use per-type C23 compoundn detection or a bounded expansion of native values, with separate binary scaling. See native compound for exceptional values, numerical acceptance, range handling and measured performance.

Exponential and logarithmic extensions

exp10, exp2m1, exp10m1, log2p1 and log10p1 support native float32/64/80/128 through detected libc routines or precision-preserving portable implementations. See their value contracts, algorithms and numericalacceptance criteria". @section autotoc_md1622 Reciprocal functions <tt>math::sec</tt>, <tt>csc</tt>, <tt>cot</tt>, <tt>sech</tt>, <tt>csch</tt> and <tt>coth</tt> support all configured native formats and MPFR. Native hyperbolic reciprocals retain representable tails beyond the overflow of cosh/sinh. See @ref "/home/runner/work/uni20/uni20/docs/tensor/native_reciprocal_math.md" "reciprocal functions and Catalan" for formulas, signed-zero behavior and numerical acceptance. @section autotoc_md1623 Digamma <tt>math::digamma(x)</tt> supports all configured native formats through compensated recurrence, Bernoulli asymptotics and pi-scaled reflection. Negative nonintegers are supported; negative integer poles give NaN. Its accuracy criterion is absolute near zeros and relative away from them. See @ref "/home/runner/work/uni20/uni20/docs/tensor/native_digamma.md" "native digamma" for regions, coefficients, tail bounds and independent validation. @section autotoc_md1624 Riemann zeta <tt>math::zeta(x)</tt> supports the whole real line at each configured native precision, including binary128. It uses Euler–Maclaurin evaluation, pole cancellation near zero, and reflection with exact pi-scaled phase reduction. Native negative even integers give positive zero; the pole at one gives positive infinity. NaN and negative infinity give NaN, and positive infinity gives one. See @ref "/home/runner/work/uni20/uni20/docs/tensor/native_zeta.md" "native zeta" for accuracy targets, truncation bounds, provider comparisons and scaled range handling. MPFR semantics are unchanged. @section autotoc_md1625 Real exponential integral <tt>math::expint(x)</tt> uses Uni20 algorithms at every configured native precision, including binary128. It combines a compensated series, a positive E1 continued fraction and positive-Ei asymptotics. Exponential scaling preserves finite results beyond exp(x) overflow and negative subnormal tails. Accuracy and performance comparisons found substantial errors in GCC13's standard Ei, justifying the selection change. See @ref "/home/runner/work/uni20/uni20/docs/tensor/native_expint.md" "native Ei" for contracts, acceptance, error control and the measured cost. MPFR scalar behavior is unchanged. @section autotoc_md1626 Real dilogarithm <tt>math::dilog_real(x)</tt> supports all configured native real types through a bounded series after real domain transformations. For x>1 it returns the real part of Li2; this is not a complex branch API. Signed zeros are retained, both infinite arguments give negative infinity, and finite arguments give finite results. See @ref "/home/runner/work/uni20/uni20/docs/tensor/native_dilogarithm.md" "native real dilogarithm" for accuracy near zero and the nonzero root, reduction formulas and validation. @section autotoc_md1627 Arithmetic-geometric mean <tt>math::agm(a,b)</tt> supports matching native real types and MPFR. Native iteration keeps separate operand exponents until a common scale is safe, including minimum-subnormal/maximum input pairs. See @ref "/home/runner/work/uni20/uni20/docs/tensor/native_agm.md" "native AGM" for domains, convergence, acceptance criteria and independent references. @section autotoc_md1628 Typed factorials <tt>uni20::factorial\<Real\>(n)</tt> explicitly selects a native floating result type. For generic precision-aware code, <tt>uni20::factorial\<Real\>(n,p)</tt> checks the native precision or selects finite/exact MPFR arithmetic. The original <tt>math::factorial(n,p)</tt> MPFR convenience remains available. The standalone header is <tt>\<uni20/core/factorial.hpp\></tt>, also included by <tt>math.hpp</tt>. See @ref "/home/runner/work/uni20/uni20/docs/tensor/native_factorial.md" "factorial contracts and accuracy". @section autotoc_md1629 Typed constants <tt>\<uni20/core/math_constants.hpp\></tt> supplies <tt>uni20::pi\<Real\></tt>, <tt>uni20::log_two\<Real\></tt>, <tt>uni20::euler_gamma\<Real\></tt> and <tt>uni20::catalan\<Real\></tt>. It is also included by <tt>math.hpp</tt>. For native types these are <tt>constexpr</tt> values of type <tt>Real</tt>, using <tt>std::numbers::pi_v\<Real\></tt>, <tt>ln2_v\<Real\></tt> and <tt>egamma_v\<Real\></tt> for the first three. Catalan is assembled from a binary prefix and independently checked; see @ref "/home/runner/work/uni20/uni20/docs/tensor/native_reciprocal_math.md" "native Catalan". This includes configured binary128 with the project's supported typed standard constants. No value is constructed by widening a double constant. @icode{cpp} constexpr auto angle = uni20::pi<long double> / 3; constexpr auto ln2 = uni20::log_two<double>; // With MPFR enabled and mpreal.hpp included: auto p = uni20::Precision::bits(256); auto wide_pi = uni20::pi<uni20::mpreal>.at(p); auto twice_pi = uni20::mpreal(2, p) * uni20::pi<uni20::mpreal>; @endicode The MPFR specializations remain precision-aware descriptors, including <tt>catalan\<mpreal\></tt>. They require finite precision, explicitly or from an approximate operand. Integral and complex constant types are not provided. Constants are validated against wider MPFR values in the shared precision matrix, including binary128 where configured. @section autotoc_md1630 Log-Gamma and Bessel providers CMake compiles and links probes independently for <tt>float</tt>, <tt>double</tt>, <tt>long double</tt>, and configured binary128. The native result retains the input's type. Unsupported signatures are constrained out; a missing extended-precision routine never falls back to double. Generic callers can use <tt>std::invocable</tt> on the corresponding <tt>uni20::math</tt> callable. These APIs accept native real arguments, not conversion-only classes or an added <tt>Precision</tt> argument. <table class="markdownTable"> <tr class="markdownTableHead"> <th class="markdownTableHeadNone"> Scalar

Log-Gamma

Bessel J/Y

float

lgammaf_r

jnf, ynf

double

lgamma_r

jn, yn

long double (including native float80)

lgammal_r

jnl, ynl

Configured float128

lgammaf128_r

jnf128, ynf128

The last row uses the platform's binary128 C math routines where available. It does not require those symbols on every platform with binary128 storage. The generated capability macros are UNI20_HAS_LGAMMA_R_{FLOAT,DOUBLE,LONG_DOUBLE,FLOAT128} and UNI20_HAS_BESSEL_{FLOAT,DOUBLE,LONG_DOUBLE,FLOAT128}. Each Bessel capability requires both J and Y. Capability detection also applies to macOS; this is not a claim that each macOS SDK supplies every routine in the table.

Log-Gamma uses only reentrant routines returning the sign through an output parameter. Neither interface reads or modifies the global signgam variable. The ordinary non-reentrant lgamma routine is not a fallback.

Numerical contracts

lgamma(x) means the natural logarithm of the absolute Gamma value. lgamma_sign(x) returns lgamma_result<R>{value, sign}. The sign is -1 or +1 where defined, zero for NaN, negative infinity and negative integer poles. At positive/negative zero the sign is +1/-1, respectively. The value follows the native provider, including infinities at poles and NaN propagation.

Bessel order comes first and must be an integer source other than bool. The common signed-64-bit argument check is followed by provider lowering: the absolute order must fit int, otherwise std::out_of_range is thrown before calling the provider. Negative orders use J_-n(x) = (-1)^n J_n(x) and Y_-n(x) = (-1)^n Y_n(x). Native Bessel Y has the provider's real domain: negative real arguments give NaN; the zero endpoint is a pole. The provider determines floating-point exceptions, range behavior and signed zeros. The MPFR overloads retain their existing signed-64-bit order domain and exact-value contracts.

Native results do not promise MPFR's correct rounding. Validation includes Gamma recurrence and half-integer identities, Bessel parity and recurrence, special values, checked order conversion, and preservation of global sign state. The shared precision matrix measures log-Gamma and Bessel accuracy at float32/64/80/128 and MPFR 128/256 bits where configured, against independent identities/series and wider MPFR references. Unavailable provider combinations are recorded in the coverage report rather than counted as passing probes.

The native coverage trajectory is tracked in issue #67 and the scalar math coverage plan.

Upper incomplete Gamma

math::upper_gamma(a,x) evaluates the upper, unregularized incomplete Gamma for matching native real types at every configured precision. Finite negative shapes are supported when x>0. The x=0 endpoint follows analytic Gamma(a), including its poles, rather than the divergent-integral interpretation for nonpositive shape. Native evaluation uses scaled series and continued fractions, with two native words where range decisions or exponent cancellation need more precision. It does not require a standard-library special-function provider or use MPFR internally. See native upper Gamma for the complete domain, accuracy targets, truncation bounds and validation evidence.

Airy Ai

math::airy_ai(x) supports all configured native real formats. It combines a four-word central series with decaying/oscillatory asymptotics. Negative-axis phase reduction starts from the exact binary input, including when a native x*sqrt(-x) would overflow; positive subnormal rounding uses a guarded path. NaN propagates and both infinities return positive zero. See native Airy Ai for the accuracy contract, reproducible constants, bounded work and validation evidence. MPFR keeps its own provider limitations; native arithmetic never calls MPFR.