math::beta(a,b) now supports configured binary128 through a Uni20 native implementation, also used when no standard Beta provider is advertised. The implementation is tested directly at float32/64/80/128. Arguments must have the same native real type and be positive. No native path uses MPFR or Boost, and enabling MPFR does not change native dispatch.

Advertised float/double/long-double standard providers remain selected. Their accuracy, range and signgam limitations remain separate from the fallback contract below; see native special functions. The existing MPFR scalar overload retains its own domain and precision rules.

Acceptance established before prototyping

For positive finite inputs the Beta function is positive and symmetric. Zero and negative parameters are outside the native domain; the new algorithm returns NaN for those and for NaN inputs. A positive infinite parameter gives the positive limiting value zero. These explicit conventions describe the new native algorithm; retained providers keep their documented provider behavior. The existing MPFR overloads and their wider domain remain unchanged.

Let eps be spacing at one, B the MPFR reference at the original represented inputs, and L=log(B). The intended error bound for finite native results is 32*eps*max(1,abs(L))*B + denorm_min/2. The logarithmic factor accounts for rounding an exponent before exponentiation. It is an explicit accuracy limitation, not a uniform few-ULP or correctly-rounded promise. The additive term accounts for final subnormal rounding; evaluate the bound with a wider reference so neither factor underflows. MPFR512 is the initial oracle, with precision increased to resolve large-parameter cancellation.

Region / all native formats	Measure	Intended criterion	Guarantee level	Observations
Ordinary positive inputs and region joins	Absolute error divided by the bound above	<=1	Empirical target	Wide scan below passes
Highly unequal inputs, large log-Gamma cancellation, tiny parameters	Same original-input measure	<=1	Empirical target; no cancellation of enormous independently rounded log-Gamma values	Wide scan below passes
Subnormal, zero and overflow boundaries	Same finite-value bound, classification checked separately	<=1; no premature intermediate overflow/underflow	Empirical accuracy target and range contract	896 extra boundary pairs per format pass
Symmetry and B(a,1)=1/a	Exact symmetry; native reciprocal comparison	Symmetric evaluation and exact reciprocal path	API contract	Native tests pass
Independent integer/half-integer cases	Exact rational or pi identities; precision ladder	At least sixteen-fold accuracy improvement on controlled cases	Demonstrated example	Shared precision tests pass

All tests use represented inputs. Do not infer accuracy at an unrounded intended argument. Separate intended criteria from measured maxima and keep retained provider measurements separate from the new algorithm's acceptance. Assume nearest rounding and gradual underflow.

Algorithm and error control

Use the Beta/Gamma identity and Stirling expansions with remainder bounds. Sort the inputs so that a<=b, and handle a parameter equal to one by direct reciprocation. If a is subnormal, the correction to 1/a+1/b is negligible compared with a result spacing throughout the supported native ranges. Sum the half-reciprocals with FMA division residuals and a compensated addition, then restore the factor of two. Plain addition of rounded reciprocals can overflow prematurely: a targeted 300000-pair float32 boundary search found 24104 such classification errors before compensation, and zero afterwards. The failing represented pair is retained as a regression, with smaller boundary searches at every native precision.

The precision-dependent transition c and number of Stirling terms are:

Significand bits	c	Terms	Absolute first-omitted-term bound for S(c)
24	8	4	<2^-37
53	16	8	<2^-70
64	16	10	<2^-80
113	32	16	<2^-136

Here S(x) is the correction in Stirling's log-Gamma expansion:

$$ \log\Gamma(x)=(x-\tfrac12)\log x-x+\tfrac12\log(2\pi)+S(x), \qquad S(x)\sim\sum_{k\ge1}\frac{B_{2k}}{2k(2k-1)x^{2k-1}}. $$

For positive arguments, the first omitted term bounds the asymptotic remainder. The table bounds truncation only; it is not a bound on total floating-point error. Coefficients are exact rational expressions rounded in the result type, generated by scripts/generate-beta-coefficients.py using Python's standard-library rational arithmetic. No provider source or coefficient table was copied. Horner evaluation uses native FMA.

When b<c, lower each parameter larger than two with B(a,b)=B(a,b-1)/(1+a/(b-1)), accumulating the factors. For the remaining small parameters evaluate

$$ B(a,b)=\left(\frac1a+\frac1b\right) \frac{\Gamma(1+a)\Gamma(1+b)}{\Gamma(1+a+b)}. $$

This keeps individual Gamma arguments small and avoids their overflow near zero. Applying the shifted identity directly near the transition failed the binary80 accuracy target because rounding 1+a+b is amplified by Gamma's derivative. The bounded recurrence removes that failure; the corresponding represented inputs are permanent regressions.

For a<c<=b, set r=a/b, l=log1p(r) and form

$$ H=\log\Gamma(1+a)-a\log b+(a-bl)-(a-\tfrac12)l+S(b)-S(a+b). $$

The central cancellation a-bl uses FMA. If r<=epsilon, its omitted value is O(a*r) and is negligible at the stated error scale. Compute 1/(a+b) as (1/b)/(1+r), never forming an overflowing sum. For a<1, return exp(H)/a; otherwise return exp(H-log(a)). Keeping the reciprocal outside the exponential in the first case prevents the observed premature overflow near 1/a=max. In that region exp(H) itself remains representable.

For c<=a<=b, use the entropy expression

$$ \log B(a,b)=a\log r-al-bl+ \tfrac12\bigl(\log(2\pi)-\log a+l\bigr)+S(a)+S(b)-S(a+b). $$

The large terms have the same sign. No difference of enormous independently rounded log-Gamma values is formed. If the negative entropy exceeds the native range, the final Beta value necessarily underflows. Only the final exponential determines the result in the ordinary finite range.

Reference calculations and observed accuracy

Tests use exact represented inputs, initially with 512-bit MPFR references. MPFR's direct Beta routine explicitly has internal overflow/underflow and potentially excessive internal-precision limitations; see its function reference. For parameters above 65536 or with binary exponent below -256, the test oracle evaluates the independent log-Gamma identity instead. Each log-Gamma term receives 512 + max(0,input_exponent) + 64 working bits; addition and subtraction use the largest of these precisions. This permits cancellation without requiring an expensive small-argument log-Gamma evaluation at the largest parameter's precision. Exponentiate before rounding back to 512 bits. Worst-case references are checked again at 1024 bits. This is test-only arithmetic; the native implementation always stays native.

The source-header scan uses a 127 by 127 rational grid, neighbors of the region transitions and reciprocal-overflow boundary, powers of two across the native exponent range, fixed unequal-parameter pairs, and 1000 deterministically generated full-range pairs plus 1000 mixed-scale pairs. On GCC 13/glibc on the local Linux host:

Format	Input pairs	Largest observed fraction of the error bound	Worst pair (decimal rendering)
float32	24899	0.700853	(8.272000312805176, 899471.4375)
float64	25119	0.906467	(1.098000000000000087, 2.229597262146031558e285)
binary80	28749	0.351563	(7, 5.266472321174442579e707)
binary128	28749	0.328125	(16, 4.602094425247528723787021288419903e310)

There were no observed accuracy or range failures. A further 896 pairs per format cover neighborhoods of independently located normal/subnormal and zero-rounding boundaries, and reciprocal-overflow neighborhoods. Their bound fractions approach one as the final result approaches a half-subnormal rounding tie; they also pass. These are sample observations, not a proof of the empirical bound throughout the domain.

The permanent numerical suite retains smaller deterministic scans, boundary searches, independent finite-product values for integer parameters and B(1/2,1/2)=pi precision improvement. It distinguishes the portable path, public fallback and existing provider probes in the coverage registry.

Cost, provider selection and validation

An illustrative GCC 13 -O2 timing loop, with 10000 calls per region, measured about 0.029–0.085 microseconds per portable float32/64 call, 0.32–6.9 for binary80 and 3.7–30.4 for binary128. Standard-provider calls at the same float32/64/80 inputs took about 0.019–0.16 microseconds. The test varied the first parameter around 1/2, 3, 15.5 or 32 and covered both comparable and million-to-one parameter ratios. A background build was running; these are orientation measurements, not a stable benchmark or performance guarantee.

The fallback is appreciably more expensive in some extended-precision regions. Existing providers remain selected; their accepted upstream signgam limitation is therefore unchanged. Their broader accuracy and range behavior must be assessed separately from the fallback's target.

For example, the existing provider loses the Gamma-ratio contribution at B(2,b) with sufficiently large finite b, whereas the independent reference is the simple reciprocal product 1/(b*(b+1)). On the scan above, the retained providers do not satisfy the new fallback's range/accuracy target. This comparison motivates a possible later provider-selection decision; it does not change the current provider contract or silently replace that selection.

All 38 targeted MPFR-enabled tests passed, including the permanent boundary searches, retained-provider checks and existing MPFR special-function regressions. The separate MPFR-disabled Release check passed 13 native tests and ran beta_example.

Full unfiltered integration validation also passed, with examples, CLI and Python bindings enabled:

Configuration	Passing CTest cases	Runtime skips
GCC13 Debug, MPFR/MPC/MPLAPACK, native128	3273	1 NUMA topology test
GCC13 Release, MPFR disabled, native128	2484	1 NUMA topology test
Clang19 Debug, MPFR/MPC/MPLAPACK, native80, native128 disabled	2936	1 NUMA topology test

Each configuration also passed the complete numerical precision executable and coverage report. CUDA runtime and macOS were not exercised. The NUMA skips reflect this host's topology and are separate from the passing counts.