Skip to content

Improve IEEE decimal transcendental accuracy #133921

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
tannergooding merged 3 commits into from
Sep 15, 2026
Merged

Improve IEEE decimal transcendental accuracy #133921

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
e1f576d
Improve IEEE decimal transcendental accuracy
tannergooding Sep 14, 2026
2cb2d0c
Cover decimal Pi cohorts and transcendental boundaries
tannergooding Sep 15, 2026
ba06161
Reduce redundant decimal transcendental work and improve log-base acc…
tannergooding Sep 15, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
64 changes: 24 additions & 40 deletions src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128Conversions.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,10 +17,9 @@ internal static partial class Number
// `UInt{64,128}Powers10` tables already present for parsing/formatting and build the required power
// of ten on the fly by chunked multiply/divide in the engine. A coefficient (< 2^113) loads into a
// binary128 significand exactly, and every 10^k with k below the format precision is exact in the
// 128-bit `ux` fraction (5^34 < 2^114), so the only rounding is the final round-to-nearest-even
// extraction of the P-digit result. That keeps the transcendental cores bit-faithful to Intel while
// the conversion stays within the <= 1 ulp faithful target; the extended-precision table path is a
// documented later refinement.
// 128-bit `ux` fraction (5^34 < 2^114). The scaled value can still round in each multiply/divide.
// Cancellation-sensitive reductions must therefore preserve small residuals in decimal before
// conversion. The result is rounded once at its final decimal quantum, including for subnormals.

/// <summary>
/// Builds a normalized <see cref="DiyFp128"/> holding the exact value of the non-zero magnitude
Expand Down Expand Up @@ -50,11 +49,17 @@ private static DiyFp128 DiyFp128ScaleByPow10<TDecimal, TValue>(DiyFp128 value, i
{
int remaining = int.Abs(power);
int maxChunk = TDecimal.Precision - 1;
int previousChunk = 0;
DiyFp128 pow = default;

while (remaining > 0)
{
int chunk = int.Min(remaining, maxChunk);
DiyFp128 pow = DiyFp128FromUInt128(UInt128.CreateTruncating(TDecimal.Power10(chunk)), 0);
if (chunk != previousChunk)
{
pow = DiyFp128FromUInt128(UInt128.CreateTruncating(TDecimal.Power10(chunk)), 0);
previousChunk = chunk;
}

if (power > 0)
{
Expand Down Expand Up @@ -112,7 +117,7 @@ private static TValue DiyFp128ToDecimal<TDecimal, TValue>(DiyFp128 value)
int binaryExponent = value._exponent - 1;
const double Log10Of2 = 0.30102999566398119521;
int d = (int)double.Floor(binaryExponent * Log10Of2);
int q = d - (precision - 1);
int q = int.Max(d - (precision - 1), TDecimal.MinAdjustedExponent);

UInt128 pow10P = UInt128.CreateTruncating(TDecimal.MaxSignificand);
pow10P++; // 10^P
Expand All @@ -132,7 +137,7 @@ private static TValue DiyFp128ToDecimal<TDecimal, TValue>(DiyFp128 value)
continue;
}

if ((coefficient != UInt128.Zero) && (coefficient < pow10Pm1))
if ((coefficient != UInt128.Zero) && (coefficient < pow10Pm1) && (q > TDecimal.MinAdjustedExponent))
{
// Under-shot (estimate was one high); pull in another decimal place.
q--;
Expand All @@ -152,6 +157,14 @@ private static TValue DiyFp128ToDecimal<TDecimal, TValue>(DiyFp128 value)
/// </summary>
private static UInt128 DiyFp128RoundToUInt128(DiyFp128 value)
{
if (value._exponent <= 0)
{
// Values below 1/2, including the tie itself, round to the even integer zero.
return ((value._exponent == 0) && ((value._hi != UxMsb) || (value._lo != 0)))
? UInt128.One
: UInt128.Zero;
}

int shift = 128 - value._exponent;
Debug.Assert(shift is > 0 and < 128);

Expand All @@ -169,14 +182,15 @@ private static UInt128 DiyFp128RoundToUInt128(DiyFp128 value)
}

/// <summary>
/// Encodes <c>(sign, coefficient, exponent)</c> into the BID bit pattern, reducing the coefficient to
/// a representable subnormal (ties-to-even) when the exponent is below the minimum and returning the
/// format's infinity when it is above the maximum.
/// Encodes an already rounded <c>(sign, coefficient, exponent)</c> into the BID bit pattern,
/// returning the format's infinity when the exponent is above the maximum.
/// </summary>
private static TValue EncodeDecimalFromUInt128<TDecimal, TValue>(bool signed, UInt128 coefficient, int exponent)
where TDecimal : unmanaged, IDecimalIeee754ParseAndFormatInfo<TDecimal, TValue>
where TValue : unmanaged, IBinaryInteger<TValue>
{
Debug.Assert(exponent >= TDecimal.MinAdjustedExponent);

if (coefficient == UInt128.Zero)
{
return DecimalIeee754FiniteNumberBinaryEncoding<TDecimal, TValue>(signed, TValue.Zero, TDecimal.MinAdjustedExponent);
Expand All @@ -187,36 +201,6 @@ private static TValue EncodeDecimalFromUInt128<TDecimal, TValue>(bool signed, UI
return signed ? TDecimal.NegativeInfinity : TDecimal.PositiveInfinity;
}

if (exponent < TDecimal.MinAdjustedExponent)
{
// Fold the extra magnitude into the coefficient as a subnormal, rounding ties-to-even.
int deficit = TDecimal.MinAdjustedExponent - exponent;

if (deficit >= UInt128.PowersOf10.Length)
{
// The coefficient has at most 34 digits, so a larger divisor rounds it entirely to zero.
return DecimalIeee754FiniteNumberBinaryEncoding<TDecimal, TValue>(signed, TValue.Zero, TDecimal.MinAdjustedExponent);
}

UInt128 power = UInt128.PowersOf10[deficit];
UInt128 quotient = coefficient / power;
UInt128 remainder = coefficient - (quotient * power);
UInt128 half = power >> 1; // 10^deficit is even, so this is an exact half

if ((remainder > half) || ((remainder == half) && UInt128.IsOddInteger(quotient)))
{
quotient++;
}

coefficient = quotient;
exponent = TDecimal.MinAdjustedExponent;

if (coefficient == UInt128.Zero)
{
return DecimalIeee754FiniteNumberBinaryEncoding<TDecimal, TValue>(signed, TValue.Zero, exponent);
}
}

return DecimalIeee754FiniteNumberBinaryEncoding<TDecimal, TValue>(signed, TValue.CreateTruncating(coefficient), exponent);
}
}
15 changes: 2 additions & 13 deletions src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128Exp.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -122,17 +122,6 @@ private static ReadOnlySpan<DiyFp128FixedCoefficient> DiyFp128FixedCoefficients(
// 1.0 as an unpacked value (Intel's UX_ONE).
private static DiyFp128 DiyFp128One => new DiyFp128(0, 1, 0x8000000000000000, 0);

// ln2 as a full unpacked value, built from the exp table's high and low pieces.
private static DiyFp128 DiyFp128Ln2
{
get
{
DiyFp128 single = default;
DiyFp128AddSub(new DiyFp128(0, 0, ExpLn2High, 0), ExpLn2Low, UxSub, new Span<DiyFp128>(ref single));
return single;
}
}

/// <summary>
/// Reduces <paramref name="orig"/> as <c>lnb*x = scale*ln2 + reduced</c> with <c>|reduced| &lt;=
/// ln2/2</c> (Intel's <c>UX_EXP_REDUCE</c>), returning <c>scale</c>. For <c>|x| &gt; 2^17</c> it
Expand Down Expand Up @@ -618,7 +607,7 @@ private static DiyFp128 DiyFp128Exp10M1(scoped in DiyFp128 argument) =>
private static DiyFp128 DiyFp128Exp2(scoped in DiyFp128 argument)
{
DiyFp128 argumentLocal = argument;
DiyFp128 ln2 = DiyFp128Ln2;
DiyFp128 ln2 = LogLn2;
DiyFp128Multiply(ref argumentLocal, ref ln2, out DiyFp128 scaled);
return DiyFp128Exp(scaled);
}
Expand All @@ -627,7 +616,7 @@ private static DiyFp128 DiyFp128Exp2(scoped in DiyFp128 argument)
private static DiyFp128 DiyFp128Exp2M1(scoped in DiyFp128 argument)
{
DiyFp128 argumentLocal = argument;
DiyFp128 ln2 = DiyFp128Ln2;
DiyFp128 ln2 = LogLn2;
DiyFp128Multiply(ref argumentLocal, ref ln2, out DiyFp128 scaled);
return DiyFp128ExpM1(scaled);
}
Expand Down
22 changes: 16 additions & 6 deletions src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128InvHyper.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,8 +17,8 @@ internal static partial class Number
// acosh(x) = log(x + sqrt(x^2 - 1)), atanh(x) = (1/2) * log((1 + x) / (1 - x)). Near the point where
// the reduced argument is 1 the naive ratio loses significance, so a small-argument path forms the
// reduced ratio directly and evaluates it with `DiyFp128LogPoly`; otherwise the big path forms the
// full argument and calls `DiyFp128Ln`. The evaluation runs entirely in the software binary128
// engine, so Decimal64/Decimal128 obtain the full ~34-digit accuracy Intel's reference does.
// full argument and calls `DiyFp128Ln`. The decimal dispatch supplies the small endpoint residuals
// for acosh and atanh before conversion to the software binary128 engine.

// Loss-of-significance thresholds (dpml_inv_hyper_x.h): the MSD boundaries selecting the small path.
private const ulong InvHyperSqrt2Over4 = 0xB504F333F9DE6484; // sqrt(2) / 4
Expand Down Expand Up @@ -64,13 +64,18 @@ private static DiyFp128 DiyFp128Asinh(DiyFp128 x)
}

/// <summary>Computes <c>acosh(x)</c> for a finite <paramref name="x"/> &gt;= 1 (Intel's <c>F_ACOSH</c>).</summary>
private static DiyFp128 DiyFp128Acosh(DiyFp128 x)
private static DiyFp128 DiyFp128Acosh(DiyFp128 x, DiyFp128 magnitudeMinusOne)
{
int exponent = x._exponent;
ulong fHi = x._hi;

Span<DiyFp128> parts = [default, default];
DiyFp128AddSub(x, DiyFp128One, UxAddSub, parts); // parts[0] = x + 1, parts[1] = x - 1
bool hasResidual = !DiyFp128IsZero(magnitudeMinusOne);
DiyFp128AddSub(x, DiyFp128One, hasResidual ? UxAdd : UxAddSub, parts);
if (hasResidual)
{
parts[1] = magnitudeMinusOne;
}

if ((exponent == 1) && (fHi <= InvHyperThreeSqrt2Over4))
{
Expand All @@ -88,7 +93,7 @@ private static DiyFp128 DiyFp128Acosh(DiyFp128 x)
}

/// <summary>Computes <c>atanh(x)</c> for a finite <paramref name="x"/> with <c>|x| &lt; 1</c> (Intel's <c>F_ATANH</c>).</summary>
private static DiyFp128 DiyFp128Atanh(DiyFp128 x)
private static DiyFp128 DiyFp128Atanh(DiyFp128 x, DiyFp128 magnitudeMinusOne)
{
uint sign = x._sign;
x._sign = 0; // |x|
Expand All @@ -105,7 +110,12 @@ private static DiyFp128 DiyFp128Atanh(DiyFp128 x)
else
{
Span<DiyFp128> parts = [default, default];
DiyFp128AddSub(x, DiyFp128One, UxAddSub, parts); // parts[0] = |x| + 1, parts[1] = |x| - 1
bool hasResidual = !DiyFp128IsZero(magnitudeMinusOne);
DiyFp128AddSub(x, DiyFp128One, hasResidual ? UxAdd : UxAddSub, parts);
if (hasResidual)
{
parts[1] = magnitudeMinusOne;
}
DiyFp128Divide(parts[1], parts[0], DiyFp128FullPrecision, out DiyFp128 ratio); // (|x| - 1) / (|x| + 1)
DiyFp128Normalize(ref ratio);
result = DiyFp128Ln(ratio); // magnitude only: log((1 - |x|) / (1 + |x|))
Expand Down
13 changes: 7 additions & 6 deletions src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128InvTrig.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -207,7 +207,7 @@ private static DiyFp128 DiyFp128Atan2(DiyFp128 y, DiyFp128 x, bool haveX)
private static DiyFp128 DiyFp128Atan(scoped in DiyFp128 arg) => DiyFp128Atan2(arg, default, false);

// UX_ASIN_ACOS with the asin/acos interval maps precomputed. Callers guarantee |arg| <= 1.
private static DiyFp128 DiyFp128AsinAcos(DiyFp128 arg, bool isAcos)
private static DiyFp128 DiyFp128AsinAcos(DiyFp128 arg, DiyFp128 magnitudeMinusOne, bool isAcos)
{
int indexMap = isAcos ? InvTrigAcosMap : InvTrigAsinMap;

Expand All @@ -223,9 +223,10 @@ private static DiyFp128 DiyFp128AsinAcos(DiyFp128 arg, bool isAcos)
{
// 1/2 <= |x| < 1: compute sqrt((1-x)/2).
exponentIncrement = 1;
DiyFp128 t = default;
DiyFp128AddSub(new DiyFp128(0, 1, UxMsb, 0), arg, UxSub | UxMagnitudeOnly, new Span<DiyFp128>(ref t));
arg = t;
arg = DiyFp128IsZero(magnitudeMinusOne)
? DiyFp128Difference(arg, DiyFp128One)
: magnitudeMinusOne;
arg._sign = 0;
arg._exponent -= 1;
arg = DiyFp128Sqrt(arg);
}
Expand Down Expand Up @@ -255,9 +256,9 @@ private static DiyFp128 DiyFp128AsinAcos(DiyFp128 arg, bool isAcos)
return value;
}

private static DiyFp128 DiyFp128Asin(scoped in DiyFp128 arg) => DiyFp128AsinAcos(arg, false);
private static DiyFp128 DiyFp128Asin(scoped in DiyFp128 arg, scoped in DiyFp128 magnitudeMinusOne) => DiyFp128AsinAcos(arg, magnitudeMinusOne, false);

private static DiyFp128 DiyFp128Acos(scoped in DiyFp128 arg) => DiyFp128AsinAcos(arg, true);
private static DiyFp128 DiyFp128Acos(scoped in DiyFp128 arg, scoped in DiyFp128 magnitudeMinusOne) => DiyFp128AsinAcos(arg, magnitudeMinusOne, true);

// True when a normalized, non-zero |arg| is strictly greater than 1 (outside the asin/acos domain).
private static bool DiyFp128MagnitudeExceedsOne(in DiyFp128 arg)
Expand Down
Loading
Loading