diff --git a/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128Conversions.cs b/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128Conversions.cs index 16969e7fb5145c..1396983a25e1d3 100644 --- a/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128Conversions.cs +++ b/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128Conversions.cs @@ -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. /// /// Builds a normalized holding the exact value of the non-zero magnitude @@ -50,11 +49,17 @@ private static DiyFp128 DiyFp128ScaleByPow10(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) { @@ -112,7 +117,7 @@ private static TValue DiyFp128ToDecimal(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 @@ -132,7 +137,7 @@ private static TValue DiyFp128ToDecimal(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--; @@ -152,6 +157,14 @@ private static TValue DiyFp128ToDecimal(DiyFp128 value) /// 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); @@ -169,14 +182,15 @@ private static UInt128 DiyFp128RoundToUInt128(DiyFp128 value) } /// - /// Encodes (sign, coefficient, exponent) 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 (sign, coefficient, exponent) into the BID bit pattern, + /// returning the format's infinity when the exponent is above the maximum. /// private static TValue EncodeDecimalFromUInt128(bool signed, UInt128 coefficient, int exponent) where TDecimal : unmanaged, IDecimalIeee754ParseAndFormatInfo where TValue : unmanaged, IBinaryInteger { + Debug.Assert(exponent >= TDecimal.MinAdjustedExponent); + if (coefficient == UInt128.Zero) { return DecimalIeee754FiniteNumberBinaryEncoding(signed, TValue.Zero, TDecimal.MinAdjustedExponent); @@ -187,36 +201,6 @@ private static TValue EncodeDecimalFromUInt128(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(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(signed, TValue.Zero, exponent); - } - } - return DecimalIeee754FiniteNumberBinaryEncoding(signed, TValue.CreateTruncating(coefficient), exponent); } } diff --git a/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128Exp.cs b/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128Exp.cs index 4ff9eb8410ee11..db4568727734eb 100644 --- a/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128Exp.cs +++ b/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128Exp.cs @@ -122,17 +122,6 @@ private static ReadOnlySpan 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(ref single)); - return single; - } - } - /// /// Reduces as lnb*x = scale*ln2 + reduced with |reduced| <= /// ln2/2 (Intel's UX_EXP_REDUCE), returning scale. For |x| > 2^17 it @@ -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); } @@ -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); } diff --git a/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128InvHyper.cs b/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128InvHyper.cs index 8a666a6844c3d3..4acab8a727c2ea 100644 --- a/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128InvHyper.cs +++ b/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128InvHyper.cs @@ -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 @@ -64,13 +64,18 @@ private static DiyFp128 DiyFp128Asinh(DiyFp128 x) } /// Computes acosh(x) for a finite >= 1 (Intel's F_ACOSH). - private static DiyFp128 DiyFp128Acosh(DiyFp128 x) + private static DiyFp128 DiyFp128Acosh(DiyFp128 x, DiyFp128 magnitudeMinusOne) { int exponent = x._exponent; ulong fHi = x._hi; Span 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)) { @@ -88,7 +93,7 @@ private static DiyFp128 DiyFp128Acosh(DiyFp128 x) } /// Computes atanh(x) for a finite with |x| < 1 (Intel's F_ATANH). - private static DiyFp128 DiyFp128Atanh(DiyFp128 x) + private static DiyFp128 DiyFp128Atanh(DiyFp128 x, DiyFp128 magnitudeMinusOne) { uint sign = x._sign; x._sign = 0; // |x| @@ -105,7 +110,12 @@ private static DiyFp128 DiyFp128Atanh(DiyFp128 x) else { Span 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|)) diff --git a/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128InvTrig.cs b/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128InvTrig.cs index 919290fada96cb..109191e2645b2d 100644 --- a/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128InvTrig.cs +++ b/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128InvTrig.cs @@ -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; @@ -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(ref t)); - arg = t; + arg = DiyFp128IsZero(magnitudeMinusOne) + ? DiyFp128Difference(arg, DiyFp128One) + : magnitudeMinusOne; + arg._sign = 0; arg._exponent -= 1; arg = DiyFp128Sqrt(arg); } @@ -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) diff --git a/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128PiTrig.cs b/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128PiTrig.cs index 36ba2af64b8729..205cda1f08e0da 100644 --- a/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128PiTrig.cs +++ b/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128PiTrig.cs @@ -13,14 +13,8 @@ internal static partial class Number // `sinpi`/`cospi`/`tanpi` from amd/aocl-libm-ose, BSD 3-Clause; see THIRD-PARTY-NOTICES.TXT): the // magnitude is split exactly into an integer and a fractional part in [0, 1), the fraction folds by // quarter turns, and a small ux sin/cos of (reduced * pi) with reduced in [0, 1/4] is evaluated. The - // integer/fractional split is exact in binary128 for every non-integer decimal (its magnitude is - // below 2^113), so the pi-scaled reduction avoids the large-argument cancellation that motivates a - // dedicated *Pi routine. The inverse variants are the radian result divided by pi. - - private static DiyFp128 UxQuarter => new DiyFp128(0, -1, UxMsb, 0); - private static DiyFp128 UxHalf => new DiyFp128(0, 0, UxMsb, 0); - private static DiyFp128 UxThreeQuarter => new DiyFp128(0, 0, 0xC000000000000000, 0); - private static DiyFp128 UxOne => new DiyFp128(0, 1, UxMsb, 0); + // reduction is performed in decimal before conversion so small distances from integers and + // half-integers are preserved. The inverse variants are the radian result divided by pi. // 0, 1/4, 1/2, 3/4, 1 -- InvTrigConstants (0, pi/4, pi/2, 3pi/4, pi) divided by pi, for the exact // signed-zero/infinity quadrant results of the inverse *Pi variants. @@ -34,96 +28,69 @@ private static DiyFp128 GetPiFractionConstant(int index) private static bool DiyFp128IsZero(in DiyFp128 value) => (value._hi | value._lo) == 0; - // Compares the magnitudes of two normalized non-negative DiyFp128 values (returns a <= b). - private static bool DiyFp128MagnitudeLessOrEqual(in DiyFp128 a, in DiyFp128 b) + private static DiyFp128 ReduceDecimalIeee754Pi( + in DecodedDecimalIeee754 decoded, out int octant) + where TDecimal : unmanaged, IDecimalIeee754ParseAndFormatInfo + where TValue : unmanaged, IBinaryInteger { - if (DiyFp128IsZero(a)) - { - return true; - } + int exponent = decoded.UnbiasedExponent; + TValue coefficient = decoded.Significand; + octant = 0; - if (DiyFp128IsZero(b)) + if (exponent >= 0) { - return false; + octant = ((exponent == 0) && TValue.IsOddInteger(coefficient)) ? 4 : 0; + return new DiyFp128(decoded.Signed ? UxSignBit : 0, UxZeroExponent, 0, 0); } - if (a._exponent != b._exponent) + if (-exponent > TDecimal.Precision) { - return a._exponent < b._exponent; + // |x| < 1/10: no reduction is needed, and 10^-exponent need not fit in TValue. + return DecimalToDiyFp128(decoded.Signed, exponent, coefficient); } - if (a._hi != b._hi) - { - return a._hi < b._hi; - } - - return a._lo <= b._lo; - } + int scale = -exponent; + TValue one = (scale == TDecimal.Precision) ? TDecimal.MaxSignificand + TValue.One : TDecimal.Power10(scale); + TValue integer = coefficient / one; + TValue fraction = coefficient - (integer * one); + octant = TValue.IsOddInteger(integer) ? 4 : 0; - // Splits |value| (assumed normalized) into its fractional part in [0, 1); reports whether floor(|value|) - // is odd and whether the value is an exact integer. - private static DiyFp128 DiyFp128SplitInteger(in DiyFp128 value, out bool oddInteger, out bool isInteger) - { - if (DiyFp128IsZero(value)) + TValue fourFraction = fraction << 2; + if (fourFraction <= one) { - oddInteger = false; - isInteger = true; - return default; + coefficient = fraction; } - - int exponent = value._exponent; - - if (exponent <= 0) + else if ((fraction << 1) <= one) { - // |value| < 1, so the whole value is fractional and floor is 0 (even). - oddInteger = false; - isInteger = false; - DiyFp128 fraction = value; - fraction._sign = 0; - return fraction; + octant += 1; + coefficient = (one >> 1) - fraction; } - - if (exponent >= 128) + else if (fourFraction <= (one + (one << 1))) { - // The 128-bit significand has no fractional bits; the value is an even integer (a power-of-two scale). - oddInteger = false; - isInteger = true; - return default; + octant += 2; + coefficient = fraction - (one >> 1); } - - UInt128 significand = new UInt128(value._hi, value._lo); - int shift = 128 - exponent; - UInt128 fractionBits = significand & ((UInt128.One << shift) - UInt128.One); - - // The integer part's low bit is bit `shift` of the significand; read it from the half that - // holds it rather than materializing the full 128-bit shifted integer for one bit. - oddInteger = (shift < 64) ? (((value._lo >> shift) & 1) != 0) - : (((value._hi >> (shift - 64)) & 1) != 0); - - if (fractionBits == UInt128.Zero) + else { - isInteger = true; - return default; + octant += 3; + coefficient = one - fraction; } - isInteger = false; - DiyFp128 result = new DiyFp128(0, exponent, fractionBits.Upper, fractionBits.Lower); - DiyFp128Normalize(ref result); - return result; + return TValue.IsZero(coefficient) + ? new DiyFp128(decoded.Signed ? UxSignBit : 0, UxZeroExponent, 0, 0) + : DecimalToDiyFp128(decoded.Signed, exponent, coefficient); } - private static DiyFp128 DiyFp128Product(in DiyFp128 a, in DiyFp128 b) + // reduced (in [0, 1/4]) * pi -> a small angle in [0, pi/4]. + private static DiyFp128 DiyFp128TimesPi(DiyFp128 reduced) { - DiyFp128 x = a; - DiyFp128 y = b; - DiyFp128Multiply(ref x, ref y, out DiyFp128 z); - DiyFp128Normalize(ref z); - return z; + reduced._sign = 0; + DiyFp128 pi = GetInvTrigConstant(4); + DiyFp128Multiply(ref reduced, ref pi, out DiyFp128 result); + DiyFp128Normalize(ref result); + return result; } - // reduced (in [0, 1/4]) * pi -> a small angle in [0, pi/4]. - private static DiyFp128 DiyFp128TimesPi(in DiyFp128 reduced) => DiyFp128Product(reduced, GetInvTrigConstant(4)); - private static DiyFp128 DiyFp128Difference(in DiyFp128 a, in DiyFp128 b) { DiyFp128 result = default; @@ -131,94 +98,85 @@ private static DiyFp128 DiyFp128Difference(in DiyFp128 a, in DiyFp128 b) return result; } - private static DiyFp128 DiyFp128WithSignFlipped(DiyFp128 value, uint sign) + private static DiyFp128 DiyFp128EvaluatePiTrig(in DiyFp128 reduced, bool cosine) { - value._sign ^= sign; - return value; + // Decimal reduction already bounds the angle to [0, pi/4], so no radian reduction is needed. + DiyFp128 angle = DiyFp128TimesPi(reduced); + Span results = [default, default]; + DiyFp128EvaluateRational(angle, cosine ? default : TrigSinCoefficients, 1, + cosine ? TrigCosCoefficients : default, 1, TrigSinCosDegree, + TrigSkip | (cosine ? TrigCosPolyFlags : TrigSinPolyFlags), results); + return results[0]; } - /// Computes sin(pi * x) for a finite non-zero binary128 argument. - private static DiyFp128 DiyFp128SinPi(in DiyFp128 x) + /// Computes sin(pi * x) from its decimal-reduced argument and octant. + private static DiyFp128 DiyFp128SinPi(in DiyFp128 reduced, int octant) { - DiyFp128 magnitude = x; - magnitude._sign = 0; - DiyFp128 fraction = DiyFp128SplitInteger(magnitude, out bool oddInteger, out bool isInteger); - - if (isInteger) - { - // sin(pi * n) = +/-0, keeping the sign of x. - return new DiyFp128(x._sign, UxZeroExponent, 0, 0); - } - - uint sign = x._sign ^ (oddInteger ? UxSignBit : 0u); + bool useCosine = (octant & 3) is 1 or 2; DiyFp128 result; - if (DiyFp128MagnitudeLessOrEqual(fraction, UxQuarter)) - { - result = DiyFp128Sin(DiyFp128TimesPi(fraction)); - } - else if (DiyFp128MagnitudeLessOrEqual(fraction, UxHalf)) + if (DiyFp128IsZero(reduced)) { - result = DiyFp128Cos(DiyFp128TimesPi(DiyFp128Difference(UxHalf, fraction))); - } - else if (DiyFp128MagnitudeLessOrEqual(fraction, UxThreeQuarter)) - { - result = DiyFp128Cos(DiyFp128TimesPi(DiyFp128Difference(fraction, UxHalf))); + if (!useCosine) + { + // sin(pi * n) = +/-0, keeping the sign of x. + return reduced; + } + + result = DiyFp128One; } else { - result = DiyFp128Sin(DiyFp128TimesPi(DiyFp128Difference(UxOne, fraction))); + result = DiyFp128EvaluatePiTrig(reduced, useCosine); } - return DiyFp128WithSignFlipped(result, sign); + result._sign = reduced._sign ^ (((octant & 4) != 0) ? UxSignBit : 0u); + return result; } - /// Computes cos(pi * x) for a finite non-zero binary128 argument. - private static DiyFp128 DiyFp128CosPi(in DiyFp128 x) + /// Computes cos(pi * x) from its decimal-reduced argument and octant. + private static DiyFp128 DiyFp128CosPi(in DiyFp128 reduced, int octant) { - DiyFp128 magnitude = x; - magnitude._sign = 0; - DiyFp128 fraction = DiyFp128SplitInteger(magnitude, out bool oddInteger, out bool isInteger); + bool useCosine = (octant & 3) is not (1 or 2); + DiyFp128 result = DiyFp128IsZero(reduced) + ? (useCosine ? DiyFp128One : new DiyFp128(0, UxZeroExponent, 0, 0)) + : DiyFp128EvaluatePiTrig(reduced, useCosine); - if (isInteger) + // cos(pi * (n + 1/2)) is exactly +0; the reduced result is +0 and must not take the odd-integer sign. + if (DiyFp128IsZero(result)) { - // cos(pi * n) = (-1)^n. - return DiyFp128WithSignFlipped(UxOne, oddInteger ? UxSignBit : 0u); + return new DiyFp128(0, UxZeroExponent, 0, 0); } - uint sign = oddInteger ? UxSignBit : 0u; - DiyFp128 result; + result._sign = (((octant + 2) & 4) != 0) ? UxSignBit : 0u; + return result; + } - if (DiyFp128MagnitudeLessOrEqual(fraction, UxQuarter)) - { - result = DiyFp128Cos(DiyFp128TimesPi(fraction)); - } - else if (DiyFp128MagnitudeLessOrEqual(fraction, UxHalf)) - { - result = DiyFp128Sin(DiyFp128TimesPi(DiyFp128Difference(UxHalf, fraction))); - } - else if (DiyFp128MagnitudeLessOrEqual(fraction, UxThreeQuarter)) + /// Computes sin(pi * x) and cos(pi * x) from their decimal-reduced argument and octant. + private static void DiyFp128SinCosPi(in DiyFp128 reduced, int octant, out DiyFp128 sin, out DiyFp128 cos) + { + if (DiyFp128IsZero(reduced)) { - result = DiyFp128WithSignFlipped(DiyFp128Sin(DiyFp128TimesPi(DiyFp128Difference(fraction, UxHalf))), UxSignBit); + sin = reduced; + cos = DiyFp128One; } else { - result = DiyFp128WithSignFlipped(DiyFp128Cos(DiyFp128TimesPi(DiyFp128Difference(UxOne, fraction))), UxSignBit); + DiyFp128 angle = DiyFp128TimesPi(reduced); + Span results = [default, default]; + DiyFp128EvaluateRational(angle, TrigSinCoefficients, 1, TrigCosCoefficients, 1, TrigSinCosDegree, + TrigSinPolyFlags | TrigCosPolyFlags | TrigNoDivide, results); + sin = results[0]; + cos = results[1]; } - // cos(pi * (n + 1/2)) is exactly +0; the reduced result is +0 and must not take the odd-integer sign. - if (DiyFp128IsZero(result)) + if ((octant & 3) is 1 or 2) { - return new DiyFp128(0, UxZeroExponent, 0, 0); + (sin, cos) = (cos, sin); } - return DiyFp128WithSignFlipped(result, sign); - } - - /// Computes sin(pi * x) and cos(pi * x) for a finite non-zero binary128 argument. - private static void DiyFp128SinCosPi(in DiyFp128 x, out DiyFp128 sin, out DiyFp128 cos) - { - sin = DiyFp128SinPi(x); - cos = DiyFp128CosPi(x); + // Integer sine zeros retain the input sign; half-integer cosine zeros are always positive. + sin._sign = reduced._sign ^ ((!DiyFp128IsZero(sin) && ((octant & 4) != 0)) ? UxSignBit : 0u); + cos._sign = (!DiyFp128IsZero(cos) && (((octant + 2) & 4) != 0)) ? UxSignBit : 0u; } } diff --git a/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128Pow.cs b/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128Pow.cs index 8dc189059510d1..7157ddeb8dbe12 100644 --- a/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128Pow.cs +++ b/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.DiyFp128Pow.cs @@ -132,7 +132,7 @@ private static void DiyFp128EvaluatePow2Polynomial(scoped in DiyFp128 hIn, out D /// Computes x^y for a positive finite (Intel's UX_POW). The caller /// handles the IEEE special cases and the sign of a negative base raised to an integer power. /// - private static DiyFp128 DiyFp128Pow(DiyFp128 x, DiyFp128 y) + private static DiyFp128 DiyFp128Pow(DiyFp128 x, DiyFp128 y, DiyFp128 xMinusOne) { Span tmp = [default, default, default]; DiyFp128 single = default; @@ -150,7 +150,13 @@ private static DiyFp128 DiyFp128Pow(DiyFp128 x, DiyFp128 y) // z = 2(g - 1) / ((g + 1) * ln2) DiyFp128 one = DiyFp128One; - DiyFp128AddSub(x, one, UxAddSub, pair); // pair[0] = g + 1, pair[1] = g - 1 + bool hasResidual = (exponent == 0) && !DiyFp128IsZero(xMinusOne); + DiyFp128AddSub(x, one, hasResidual ? UxAdd : UxAddSub, pair); + if (hasResidual) + { + // A large y amplifies conversion error in x - 1; retain the decimal residual. + pair[1] = xMinusOne; + } tmp[0] = pair[0]; tmp[1] = pair[1]; diff --git a/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.Transcendental.cs b/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.Transcendental.cs index f1a98c135e94c8..152672bd478522 100644 --- a/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.Transcendental.cs +++ b/src/libraries/System.Private.CoreLib/src/System/Number.DecimalIeee754.Transcendental.cs @@ -20,6 +20,56 @@ internal static partial class Number private static bool DecimalIeee754UsesDouble() => false; + // The log/pow reduction writes x = 2^n * g. Only n == 0 can make log(x) small, so other + // intervals need no decimal residual. Use the engine's existing 1/sqrt(2) boundary. + private static bool DiyFp128LogNeedsResidual(in DiyFp128 argument) + => argument._exponent == ((argument._hi <= LogOneOverSqrt2) ? 1 : 0); + + // Restores exact dyadic operands before domain checks and binary subtraction. Otherwise returns + // the nonzero decimal residual where cancellation is possible; zero means to use the binary path. + private static DiyFp128 DecimalIeee754MagnitudeMinusOne( + in DecodedDecimalIeee754 decoded, ref DiyFp128 argument) + where TDecimal : unmanaged, IDecimalIeee754ParseAndFormatInfo + where TValue : unmanaged, IBinaryInteger + { + if (argument._exponent is 0 or 1) + { + // For 1/2 <= |x| < 2, subtraction is exact in the operand's decimal quantum. + int scale = -decoded.UnbiasedExponent; + + TValue one = (scale == TDecimal.Precision) ? TDecimal.MaxSignificand + TValue.One : TDecimal.Power10(scale); + + // An exact dyadic input has coefficient n * 5^scale and value n * 2^-scale. + // Division can leave guard bits even for such inputs. Round to a 64-bit candidate, + // check its alignment to 2^-scale, then verify its decimal coefficient exactly. + // exponent + scale is in [1, 35]; a carry out of hi represents 2^exponent. + ulong hi = argument._hi + (argument._lo >> 63); + int shift = argument._exponent + scale; + if ((hi << shift) == 0) + { + ulong integer = (hi == 0) ? 1UL << shift : hi >> (64 - shift); + TValue coefficient = TValue.CreateTruncating(integer) * (one >> scale); + if (coefficient == decoded.Significand) + { + argument._hi = (hi == 0) ? UxMsb : hi; + argument._lo = 0; + argument._exponent += (hi == 0) ? 1 : 0; + return default; + } + } + + bool negative = decoded.Significand < one; + TValue difference = negative ? one - decoded.Significand : decoded.Significand - one; + + DiyFp128 numerator = DiyFp128FromUInt128(UInt128.CreateTruncating(difference), negative ? UxSignBit : 0); + DiyFp128 denominator = DiyFp128FromUInt128(UInt128.CreateTruncating(one), 0); + DiyFp128Divide(numerator, denominator, DiyFp128FullPrecision, out DiyFp128 result); + return result; + } + + return default; + } + /// Computes e^x. internal static TValue ExpDecimalIeee754(TValue x) where TDecimal : unmanaged, IDecimalIeee754ParseAndFormatInfo @@ -311,8 +361,24 @@ internal static TValue LogDecimalIeee754(TValue x) return ConvertFloatToDecimalIeee754(double.Log(value)); } + return DiyFp128ToDecimal(LogDecimalIeee754(decoded)); + } + + private static DiyFp128 LogDecimalIeee754(in DecodedDecimalIeee754 decoded) + where TDecimal : unmanaged, IDecimalIeee754ParseAndFormatInfo + where TValue : unmanaged, IBinaryInteger + { DiyFp128 argument = DecimalToDiyFp128(decoded.Signed, decoded.UnbiasedExponent, decoded.Significand); - return DiyFp128ToDecimal(DiyFp128Ln(argument)); + if (DiyFp128LogNeedsResidual(argument)) + { + DiyFp128 residual = DecimalIeee754MagnitudeMinusOne(decoded, ref argument); + if (!DiyFp128IsZero(residual)) + { + return DiyFp128Ln1p(residual); + } + } + + return DiyFp128Ln(argument); } /// Computes log_newBase(x) as log(x) / log(newBase), mirroring the @@ -341,19 +407,33 @@ internal static TValue LogDecimalIeee754(TValue x, TValue newB } DecodedDecimalIeee754 decodedX = UnpackDecimalIeee754(x); - bool xIsOne = !TDecimal.IsInfinity(x) && !TDecimal.IsNegative(x) - && DecimalIeee754MagnitudeIsOne(decodedX.UnbiasedExponent, decodedX.Significand); bool baseIsZero = !TDecimal.IsInfinity(newBase) && TValue.IsZero(decodedBase.Significand); bool baseIsPositiveInfinity = TDecimal.IsInfinity(newBase) && !TDecimal.IsNegative(newBase); - if (!xIsOne && (baseIsZero || baseIsPositiveInfinity)) + if (baseIsZero || baseIsPositiveInfinity) { - return TDecimal.NaNMask; + bool xIsOne = !TDecimal.IsInfinity(x) && !TDecimal.IsNegative(x) + && DecimalIeee754MagnitudeIsOne(decodedX.UnbiasedExponent, decodedX.Significand); + if (!xIsOne) + { + return TDecimal.NaNMask; + } } - TValue logX = LogDecimalIeee754(x); - TValue logBase = LogDecimalIeee754(newBase); - return DivideDecimalIeee754(logX, logBase); + if (DecimalIeee754UsesDouble() || TDecimal.IsInfinity(x) || TDecimal.IsInfinity(newBase) + || decodedX.Signed || decodedBase.Signed || TValue.IsZero(decodedX.Significand) || baseIsZero) + { + TValue logX = LogDecimalIeee754(x); + TValue logBase = LogDecimalIeee754(newBase); + return DivideDecimalIeee754(logX, logBase); + } + + // Keep both logarithms and their quotient wide; rounding them to decimal first loses + // accuracy even when each logarithm is individually correctly rounded. + DiyFp128 numerator = LogDecimalIeee754(decodedX); + DiyFp128 denominator = LogDecimalIeee754(decodedBase); + DiyFp128Divide(numerator, denominator, DiyFp128FullPrecision, out DiyFp128 result); + return DiyFp128ToDecimal(result); } /// Computes log2(x). @@ -393,6 +473,15 @@ internal static TValue Log2DecimalIeee754(TValue x) } DiyFp128 argument = DecimalToDiyFp128(decoded.Signed, decoded.UnbiasedExponent, decoded.Significand); + if (DiyFp128LogNeedsResidual(argument)) + { + DiyFp128 residual = DecimalIeee754MagnitudeMinusOne(decoded, ref argument); + if (!DiyFp128IsZero(residual)) + { + return DiyFp128ToDecimal(DiyFp128Log2P1(residual)); + } + } + return DiyFp128ToDecimal(DiyFp128Log2(argument)); } @@ -433,6 +522,15 @@ internal static TValue Log10DecimalIeee754(TValue x) } DiyFp128 argument = DecimalToDiyFp128(decoded.Signed, decoded.UnbiasedExponent, decoded.Significand); + if (DiyFp128LogNeedsResidual(argument)) + { + DiyFp128 residual = DecimalIeee754MagnitudeMinusOne(decoded, ref argument); + if (!DiyFp128IsZero(residual)) + { + return DiyFp128ToDecimal(DiyFp128Log10P1(residual)); + } + } + return DiyFp128ToDecimal(DiyFp128Log10(argument)); } @@ -514,21 +612,32 @@ private static TValue Log1pDecimalIeee754(TValue x, LogBase lo DiyFp128 argument = DecimalToDiyFp128(decoded.Signed, decoded.UnbiasedExponent, decoded.Significand); - // Guard the 1 + x domain in the binary128 engine (the double path gets this from IEEE): the - // conversion error is always below the decimal granularity near x = -1, so 1 + x is exact here. - DiyFp128 onePlus = default; - DiyFp128AddSub(DiyFp128One, argument, UxAdd, new Span(ref onePlus)); - - if ((onePlus._hi | onePlus._lo) == 0) + if (decoded.Signed && (argument._exponent >= 0)) { - // logP1(-1) = -inf. - return TDecimal.NegativeInfinity; - } + DiyFp128 onePlus = DecimalIeee754MagnitudeMinusOne(decoded, ref argument); - if (onePlus._sign != 0) - { - // logP1(x < -1) is invalid and produces the canonical quiet NaN. - return TDecimal.NaNMask; + if (DiyFp128MagnitudeExceedsOne(argument)) + { + return TDecimal.NaNMask; + } + + if (DiyFp128MagnitudeIsOne(argument)) + { + return TDecimal.NegativeInfinity; + } + + if (!DiyFp128IsZero(onePlus)) + { + onePlus._sign ^= UxSignBit; + + DiyFp128 logarithm = logBase switch + { + LogBase.Two => DiyFp128Log2(onePlus), + LogBase.Ten => DiyFp128Log10(onePlus), + _ => DiyFp128Ln(onePlus), + }; + return DiyFp128ToDecimal(logarithm); + } } DiyFp128 result128 = logBase switch @@ -674,7 +783,7 @@ internal static TValue PowDecimalIeee754(TValue x, TValue y) bool yIsOddInteger = false; bool yIsInteger = false; - if (!yInf) + if (!yInf && TDecimal.IsNegative(x)) { yIsInteger = DecimalIeee754IsInteger(dy.UnbiasedExponent, dy.Significand, out yIsOddInteger); } @@ -762,7 +871,10 @@ internal static TValue PowDecimalIeee754(TValue x, TValue y) // The engine evaluates |x|^y; a negative base with an odd integer exponent carries the sign. DiyFp128 baseValue = DecimalToDiyFp128(signed: false, dx.UnbiasedExponent, dx.Significand); DiyFp128 exponentValue = DecimalToDiyFp128(dy.Signed, dy.UnbiasedExponent, dy.Significand); - DiyFp128 magnitude = DiyFp128Pow(baseValue, exponentValue); + DiyFp128 baseMinusOne = DiyFp128LogNeedsResidual(baseValue) + ? DecimalIeee754MagnitudeMinusOne(dx, ref baseValue) + : default; + DiyFp128 magnitude = DiyFp128Pow(baseValue, exponentValue, baseMinusOne); if (dx.Signed && yIsOddInteger) { @@ -926,6 +1038,14 @@ internal static TValue RootNDecimalIeee754(TValue x, int n) return TDecimal.NaNMask; } + if (n == 1) + { + // Retain the full-precision result cohort without a binary conversion or approximation. + int padding = int.Min(TDecimal.Precision - TDecimal.CountDigits(dx.Significand), dx.UnbiasedExponent - TDecimal.MinAdjustedExponent); + return DecimalIeee754FiniteNumberBinaryEncoding( + dx.Signed, dx.Significand * TDecimal.Power10(padding), dx.UnbiasedExponent - padding); + } + if (DecimalIeee754UsesDouble()) { double value = ConvertDecimalIeee754ToFloat(x); @@ -939,15 +1059,17 @@ internal static TValue RootNDecimalIeee754(TValue x, int n) return ConvertFloatToDecimalIeee754(result); } - // The engine evaluates |x|^(1/n) with the reciprocal formed exactly in the binary128 domain; + // The engine evaluates |x|^(1/n) with the reciprocal formed in binary128 working precision; // a negative base only reaches here with an odd n, so it simply carries the sign. `n` is taken // through `long` so `int.MinValue`'s magnitude does not overflow. - DiyFp128 one = new DiyFp128(0u, 1, 0x8000_0000_0000_0000, 0); - DiyFp128 degree = DecimalToDiyFp128(nNegative, 0, TValue.CreateTruncating(long.Abs(n))); - DiyFp128Divide(one, degree, DiyFp128FullPrecision, out DiyFp128 exponent); + DiyFp128 degree = DiyFp128FromWord(n); + DiyFp128Divide(DiyFp128One, degree, DiyFp128FullPrecision, out DiyFp128 exponent); DiyFp128 baseValue = DecimalToDiyFp128(signed: false, dx.UnbiasedExponent, dx.Significand); - DiyFp128 magnitude = DiyFp128Pow(baseValue, exponent); + DiyFp128 baseMinusOne = DiyFp128LogNeedsResidual(baseValue) + ? DecimalIeee754MagnitudeMinusOne(dx, ref baseValue) + : default; + DiyFp128 magnitude = DiyFp128Pow(baseValue, exponent, baseMinusOne); if (dx.Signed) { @@ -1190,11 +1312,14 @@ internal static TValue AsinDecimalIeee754(TValue x) DiyFp128 argument = DecimalToDiyFp128(decoded.Signed, decoded.UnbiasedExponent, decoded.Significand); + DiyFp128 magnitudeMinusOne = DecimalIeee754MagnitudeMinusOne(decoded, ref argument); + if (DiyFp128MagnitudeExceedsOne(argument)) { return TDecimal.NaNMask; } - return DiyFp128ToDecimal(DiyFp128Asin(argument)); + + return DiyFp128ToDecimal(DiyFp128Asin(argument, magnitudeMinusOne)); } /// Computes acos(x), the result in radians. @@ -1235,11 +1360,14 @@ internal static TValue AcosDecimalIeee754(TValue x) DiyFp128 argument = DecimalToDiyFp128(decoded.Signed, decoded.UnbiasedExponent, decoded.Significand); + DiyFp128 magnitudeMinusOne = DecimalIeee754MagnitudeMinusOne(decoded, ref argument); + if (DiyFp128MagnitudeExceedsOne(argument)) { return TDecimal.NaNMask; } - return DiyFp128ToDecimal(DiyFp128Acos(argument)); + + return DiyFp128ToDecimal(DiyFp128Acos(argument, magnitudeMinusOne)); } /// Computes atan2(y, x), the angle of the vector (x, y) in radians. @@ -1337,8 +1465,8 @@ internal static TValue SinPiDecimalIeee754(TValue x) return ConvertFloatToDecimalIeee754(double.SinPi(value)); } - DiyFp128 argument = DecimalToDiyFp128(decoded.Signed, decoded.UnbiasedExponent, decoded.Significand); - return DiyFp128ToDecimal(DiyFp128SinPi(argument)); + DiyFp128 argument = ReduceDecimalIeee754Pi(decoded, out int octant); + return DiyFp128ToDecimal(DiyFp128SinPi(argument, octant)); } /// Computes cos(pi * x). @@ -1371,8 +1499,8 @@ internal static TValue CosPiDecimalIeee754(TValue x) return ConvertFloatToDecimalIeee754(double.CosPi(value)); } - DiyFp128 argument = DecimalToDiyFp128(decoded.Signed, decoded.UnbiasedExponent, decoded.Significand); - return DiyFp128ToDecimal(DiyFp128CosPi(argument)); + DiyFp128 argument = ReduceDecimalIeee754Pi(decoded, out int octant); + return DiyFp128ToDecimal(DiyFp128CosPi(argument, octant)); } /// Computes tan(pi * x). @@ -1405,8 +1533,8 @@ internal static TValue TanPiDecimalIeee754(TValue x) return ConvertFloatToDecimalIeee754(double.TanPi(value)); } - DiyFp128 argument = DecimalToDiyFp128(decoded.Signed, decoded.UnbiasedExponent, decoded.Significand); - DiyFp128SinCosPi(argument, out DiyFp128 sin, out DiyFp128 cos); + DiyFp128 argument = ReduceDecimalIeee754Pi(decoded, out int octant); + DiyFp128SinCosPi(argument, octant, out DiyFp128 sin, out DiyFp128 cos); if (DiyFp128IsZero(cos)) { @@ -1414,6 +1542,11 @@ internal static TValue TanPiDecimalIeee754(TValue x) return (sin._sign != 0) ? TDecimal.NegativeInfinity : TDecimal.PositiveInfinity; } + if (DiyFp128IsZero(sin)) + { + return DecimalIeee754FiniteNumberBinaryEncoding(sin.IsNegative ^ cos.IsNegative, TValue.Zero, 0); + } + DiyFp128Divide(sin, cos, DiyFp128FullPrecision, out DiyFp128 tangent); return DiyFp128ToDecimal(tangent); } @@ -1453,8 +1586,8 @@ internal static (TValue SinPi, TValue CosPi) SinCosPiDecimalIeee754(cosValue)); } - DiyFp128 argument = DecimalToDiyFp128(decoded.Signed, decoded.UnbiasedExponent, decoded.Significand); - DiyFp128SinCosPi(argument, out DiyFp128 sin, out DiyFp128 cos); + DiyFp128 argument = ReduceDecimalIeee754Pi(decoded, out int octant); + DiyFp128SinCosPi(argument, octant, out DiyFp128 sin, out DiyFp128 cos); return (DiyFp128ToDecimal(sin), DiyFp128ToDecimal(cos)); } @@ -1536,12 +1669,14 @@ internal static TValue AsinPiDecimalIeee754(TValue x) DiyFp128 argument = DecimalToDiyFp128(decoded.Signed, decoded.UnbiasedExponent, decoded.Significand); + DiyFp128 magnitudeMinusOne = DecimalIeee754MagnitudeMinusOne(decoded, ref argument); + if (DiyFp128MagnitudeExceedsOne(argument)) { return TDecimal.NaNMask; } - DiyFp128Divide(DiyFp128Asin(argument), GetInvTrigConstant(4), DiyFp128FullPrecision, out DiyFp128 quotient); + DiyFp128Divide(DiyFp128Asin(argument, magnitudeMinusOne), GetInvTrigConstant(4), DiyFp128FullPrecision, out DiyFp128 quotient); return DiyFp128ToDecimal(quotient); } @@ -1583,12 +1718,14 @@ internal static TValue AcosPiDecimalIeee754(TValue x) DiyFp128 argument = DecimalToDiyFp128(decoded.Signed, decoded.UnbiasedExponent, decoded.Significand); + DiyFp128 magnitudeMinusOne = DecimalIeee754MagnitudeMinusOne(decoded, ref argument); + if (DiyFp128MagnitudeExceedsOne(argument)) { return TDecimal.NaNMask; } - DiyFp128Divide(DiyFp128Acos(argument), GetInvTrigConstant(4), DiyFp128FullPrecision, out DiyFp128 quotient); + DiyFp128Divide(DiyFp128Acos(argument, magnitudeMinusOne), GetInvTrigConstant(4), DiyFp128FullPrecision, out DiyFp128 quotient); return DiyFp128ToDecimal(quotient); } @@ -1829,11 +1966,14 @@ internal static TValue AcoshDecimalIeee754(TValue x) DiyFp128 argument = DecimalToDiyFp128(decoded.Signed, decoded.UnbiasedExponent, decoded.Significand); + DiyFp128 magnitudeMinusOne = DecimalIeee754MagnitudeMinusOne(decoded, ref argument); + if (!DiyFp128MagnitudeExceedsOne(argument) && !DiyFp128MagnitudeIsOne(argument)) { return TDecimal.NaNMask; } - return DiyFp128ToDecimal(DiyFp128Acosh(argument)); + + return DiyFp128ToDecimal(DiyFp128Acosh(argument, magnitudeMinusOne)); } /// Computes atanh(x). @@ -1870,6 +2010,8 @@ internal static TValue AtanhDecimalIeee754(TValue x) DiyFp128 argument = DecimalToDiyFp128(decoded.Signed, decoded.UnbiasedExponent, decoded.Significand); + DiyFp128 magnitudeMinusOne = DecimalIeee754MagnitudeMinusOne(decoded, ref argument); + if (DiyFp128MagnitudeIsOne(argument)) { // atanh(+/-1) = +/-inf (pole). @@ -1882,6 +2024,6 @@ internal static TValue AtanhDecimalIeee754(TValue x) return TDecimal.NaNMask; } - return DiyFp128ToDecimal(DiyFp128Atanh(argument)); + return DiyFp128ToDecimal(DiyFp128Atanh(argument, magnitudeMinusOne)); } } diff --git a/src/libraries/System.Runtime/tests/System.Runtime.Tests/System/Decimal128Tests.cs b/src/libraries/System.Runtime/tests/System.Runtime.Tests/System/Decimal128Tests.cs index a6fdb9b823733a..8231bc3f9c90e5 100644 --- a/src/libraries/System.Runtime/tests/System.Runtime.Tests/System/Decimal128Tests.cs +++ b/src/libraries/System.Runtime/tests/System.Runtime.Tests/System/Decimal128Tests.cs @@ -12,6 +12,289 @@ namespace System.Tests { public class Decimal128Tests { + [Theory] + [InlineData("0", "0", "1", "0")] + [InlineData("0.5", "1", "0", "Infinity")] + [InlineData("1", "0", "-1", "-0")] + [InlineData("1.5", "-1", "0", "-Infinity")] + [InlineData("2", "0", "1", "0")] + [InlineData("2.5", "1", "0", "Infinity")] + [InlineData("1000", "0", "1", "0")] + [InlineData("1000.5", "1", "0", "Infinity")] + [InlineData("1001", "0", "-1", "-0")] + [InlineData("1001.5", "-1", "0", "-Infinity")] + public static void PiExactCohortTest(string input, string expectedSin, string expectedCos, string expectedTan) + { + Test(7, input, expectedSin, expectedCos, expectedTan); + Test(16, input, expectedSin, expectedCos, expectedTan); + Test(34, input, expectedSin, expectedCos, expectedTan); + + static void Test(int precision, string input, string expectedSin, string expectedCos, string expectedTan) + where T : IFloatingPointIeee754 + { + T sin = T.Parse(expectedSin, CultureInfo.InvariantCulture); + T cos = T.Parse(expectedCos, CultureInfo.InvariantCulture); + T tan = T.Parse(expectedTan, CultureInfo.InvariantCulture); + int digits = input.Replace(".", "").TrimStart('0').Length; + + for (int padding = 0; padding <= precision - digits; padding++) + { + string padded = input + ((padding > 0 && !input.Contains('.')) ? "." : "") + new string('0', padding); + T x = T.Parse(padded, CultureInfo.InvariantCulture); + Check(x, sin, cos, tan); + Check(-x, -sin, cos, -tan); + } + + static void Check(T x, T sin, T cos, T tan) + { + AssertResult(sin, T.SinPi(x)); + AssertResult(cos, T.CosPi(x)); + AssertResult(tan, T.TanPi(x)); + (T actualSin, T actualCos) = T.SinCosPi(x); + AssertResult(sin, actualSin); + AssertResult(cos, actualCos); + } + + static void AssertResult(T expected, T actual) + { + Assert.Equal(expected, actual); + Assert.Equal(T.IsNegative(expected), T.IsNegative(actual)); + } + } + } + + [Theory] + [InlineData(32, 7)] + [InlineData(64, 16)] + [InlineData(128, 34)] + public static void TranscendentalUnitCohortTest(int width, int precision) + { + switch (width) + { + case 32: + Test(precision); + break; + case 64: + Test(precision); + break; + case 128: + Test(precision); + break; + default: + throw new InvalidOperationException($"Unexpected width '{width}'."); + } + + static void Test(int precision) where T : IFloatingPointIeee754 + { + T half = T.CreateChecked(0.5); + for (int padding = 0; padding < precision; padding++) + { + T one = T.Parse("1." + new string('0', padding), CultureInfo.InvariantCulture); + Assert.Equal(T.Zero, T.Acosh(one)); + Assert.Equal(T.Zero, T.Acos(one)); + Assert.Equal(T.Zero, T.AcosPi(one)); + Assert.Equal(T.One, T.AcosPi(-one)); + Assert.Equal(half, T.AsinPi(one)); + Assert.Equal(-half, T.AsinPi(-one)); + Assert.Equal(T.PositiveInfinity, T.Atanh(one)); + Assert.Equal(T.NegativeInfinity, T.Atanh(-one)); + Assert.Equal(T.Zero, T.Log(one)); + Assert.Equal(T.Zero, T.Log2(one)); + Assert.Equal(T.Zero, T.Log10(one)); + Assert.Equal(T.NegativeInfinity, T.LogP1(-one)); + Assert.Equal(T.NegativeInfinity, T.Log2P1(-one)); + Assert.Equal(T.NegativeInfinity, T.Log10P1(-one)); + Assert.Equal(T.One, T.Pow(one, T.CreateChecked(1000))); + Assert.Equal(T.One, T.RootN(one, 3)); + } + } + } + + [Theory] + [InlineData(nameof(Decimal128.Acos), "0.5", "1.047197551196597746154214461093167628065723133125")] + [InlineData(nameof(Decimal128.Acos), "0.75", "0.72273424781341561117837735264133336202521848642444")] + [InlineData(nameof(Decimal128.Asin), "0.75", "0.84806207898148100805294433899841808007336621326311")] + [InlineData(nameof(Decimal128.Atanh), "0.75", "0.97295507452765665255267637172158986481854236479093")] + [InlineData(nameof(Decimal128.Log), "0.75", "-0.28768207245178092743921900599382743150350971089776")] + [InlineData(nameof(Decimal128.Log2), "0.75", "-0.41503749927884381854626105605218349124018559230752")] + [InlineData(nameof(Decimal128.Log10), "0.75", "-0.12493873660829995313244988619387074433625089873352")] + [InlineData(nameof(Decimal128.LogP1), "0.75", "-1.3862943611198906188344642429163531361510002687205")] + [InlineData(nameof(Decimal128.Log2P1), "0.75", "-2")] + [InlineData(nameof(Decimal128.Log10P1), "0.75", "-0.60205999132796239042747778944898605353637976292422")] + [InlineData(nameof(Decimal128.Pow), "0.75", "0.69795364432657469920591406023742556581340316925824")] + [InlineData(nameof(Decimal128.Acosh), "1.25", "0.69314718055994530941723212145817656807550013436026")] + [InlineData(nameof(Decimal128.Log), "1.25", "0.22314355131420975576629509030983450337460108554801")] + [InlineData(nameof(Decimal128.Acos), "0.7", "0.79539883018414355549096833892476432854279596104639")] + [InlineData(nameof(Decimal128.Asin), "0.7", "0.77539749661075306374035335271498711355578873864116")] + [InlineData(nameof(Decimal128.Atanh), "0.7", "0.8673005276940531944271446904753004154703562273815")] + [InlineData(nameof(Decimal128.Log), "0.7", "-0.35667494393873237891263871124118447796401675904691")] + [InlineData(nameof(Decimal128.Log2), "0.7", "-0.51457317282975824042835011225755936722380476705844")] + [InlineData(nameof(Decimal128.Log10), "0.7", "-0.15490195998574316928778374140736380651642760367603")] + [InlineData(nameof(Decimal128.LogP1), "0.7", "-1.203972804325935992622746217761838502953610930806")] + [InlineData(nameof(Decimal128.Log2P1), "0.7", "-1.7369655941662061664165804855415736671050169853321")] + [InlineData(nameof(Decimal128.Log10P1), "0.7", "-0.5228787452803375627049720967448846907998711358093")] + [InlineData(nameof(Decimal128.Pow), "0.7", "0.64028385346008610598026830743665121678041972366247")] + [InlineData(nameof(Decimal128.Acosh), "1.3", "0.75643291085695958624207680696874177560093336091494")] + [InlineData(nameof(Decimal128.Log), "1.3", "0.26236426446749105203549598688095439720416645613143")] + public static void TranscendentalCohortAccuracyTest(string operation, string input, string oracle) + { + Test(7, operation, input, oracle); + Test(16, operation, input, oracle); + Test(34, operation, input, oracle); + + static void Test(int precision, string operation, string input, string oracle) + where T : IFloatingPointIeee754 + { + T expected = T.Parse(oracle, CultureInfo.InvariantCulture); + T ulp = T.Max(T.Abs(T.BitIncrement(expected) - expected), T.Abs(expected - T.BitDecrement(expected))); + int digits = input.Replace(".", "").TrimStart('0').Length; + + for (int padding = 0; padding <= precision - digits; padding++) + { + T x = T.Parse(input + new string('0', padding), CultureInfo.InvariantCulture); + T actual = operation switch + { + nameof(Decimal128.Acos) => T.Acos(x), + nameof(Decimal128.Asin) => T.Asin(x), + nameof(Decimal128.Acosh) => T.Acosh(x), + nameof(Decimal128.Atanh) => T.Atanh(x), + nameof(Decimal128.Log) => T.Log(x), + nameof(Decimal128.Log2) => T.Log2(x), + nameof(Decimal128.Log10) => T.Log10(x), + nameof(Decimal128.LogP1) => T.LogP1(-x), + nameof(Decimal128.Log2P1) => T.Log2P1(-x), + nameof(Decimal128.Log10P1) => T.Log10P1(-x), + nameof(Decimal128.Pow) => T.Pow(x, T.CreateChecked(1.25)), + _ => throw new InvalidOperationException($"Unexpected operation '{operation}'."), + }; + Assert.True(T.Abs(actual - expected) <= ulp, + $"{typeof(T).Name}.{operation}({x}): expected {expected}, actual {actual}, padding {padding}"); + } + } + } + + [Theory] + [InlineData(nameof(Decimal128.Asin), "0.9999999999999999999999999999999999", "1.5707963267948966050891860679088009540817")] + [InlineData(nameof(Decimal128.Acos), "0.9999999999999999999999999999999999", "1.41421356237309504880168872420969809035478e-17")] + [InlineData(nameof(Decimal128.Acos), "-0.9999999999999999999999999999999999", "3.14159265358979322432050775954855239618028")] + [InlineData(nameof(Decimal128.Acosh), "1.000000000000000000000000000000001", "4.47213595499957939281834733746255209820324e-17")] + [InlineData(nameof(Decimal128.Atanh), "0.9999999999999999999999999999999999", "39.4905201711787492830144707903632797882565")] + [InlineData(nameof(Decimal128.Log), "0.9999999999999999999999999999999999", "-1.00000000000000000000000000000000005e-34")] + [InlineData(nameof(Decimal128.Log), "1.000000000000000000000000000000001", "9.999999999999999999999999999999995e-34")] + [InlineData(nameof(Decimal128.Log2), "0.9999999999999999999999999999999999", "-1.4426950408889634073599246810018922095614e-34")] + [InlineData(nameof(Decimal128.Log10), "0.9999999999999999999999999999999999", "-4.34294481903251827651128918916605104009121e-35")] + [InlineData(nameof(Decimal128.LogP1), "-0.9999999999999999999999999999999999", "-78.2878931617975532566117094592683830584375")] + [InlineData(nameof(Decimal128.Log2P1), "-0.9999999999999999999999999999999999", "-112.945555226170319827590860602639265979404")] + [InlineData(nameof(Decimal128.Log10P1), "-0.9999999999999999999999999999999999", "-34.00000000000000000000000000000000")] + [InlineData(nameof(Decimal128.AsinPi), "0.9999999999999999999999999999999999", "0.499999999999999995498418419214469652224004")] + [InlineData(nameof(Decimal128.AcosPi), "0.9999999999999999999999999999999999", "4.50158158078553034777599595503370295083605e-18")] + [InlineData(nameof(Decimal128.SinPi), "0.9999999999999999999999999999999999", "3.14159265358979323846264338327950288419717e-34")] + [InlineData(nameof(Decimal128.SinPi), "1.000000000000000000000000000000001", "-3.14159265358979323846264338327950288419717e-33")] + [InlineData(nameof(Decimal128.SinPi), "10000000000000000000000000.00000001", "3.14159265358979272169136537828252546123232e-8")] + [InlineData(nameof(Decimal128.CosPi), "0.4999999999999999999999999999999999", "3.14159265358979323846264338327950288419717e-34")] + [InlineData(nameof(Decimal128.CosPi), "0.5000000000000000000000000000000001", "-3.14159265358979323846264338327950288419717e-34")] + [InlineData(nameof(Decimal128.TanPi), "0.4999999999999999999999999999999999", "3183098861837906715377675267450287.24068919")] + [InlineData(nameof(Decimal128.Exp2), "20413.2481430828416276631128942123", "9.99999999999999999999999999998176545067148e6144")] + [InlineData(nameof(Decimal128.Exp2M1), "20413.2481430828416276631128942123", "9.99999999999999999999999999998176545067148e6144")] + public static void TranscendentalBoundaryAccuracyTest(string operation, string input, string oracle) + { + Decimal128 x = Decimal128.Parse(input, CultureInfo.InvariantCulture); + Decimal128 expected = Decimal128.Parse(oracle, CultureInfo.InvariantCulture); + Decimal128 actual = operation switch + { + nameof(Decimal128.Asin) => Decimal128.Asin(x), + nameof(Decimal128.Acos) => Decimal128.Acos(x), + nameof(Decimal128.Acosh) => Decimal128.Acosh(x), + nameof(Decimal128.Atanh) => Decimal128.Atanh(x), + nameof(Decimal128.Log) => Decimal128.Log(x), + nameof(Decimal128.Log2) => Decimal128.Log2(x), + nameof(Decimal128.Log10) => Decimal128.Log10(x), + nameof(Decimal128.LogP1) => Decimal128.LogP1(x), + nameof(Decimal128.Log2P1) => Decimal128.Log2P1(x), + nameof(Decimal128.Log10P1) => Decimal128.Log10P1(x), + nameof(Decimal128.AsinPi) => Decimal128.AsinPi(x), + nameof(Decimal128.AcosPi) => Decimal128.AcosPi(x), + nameof(Decimal128.SinPi) => Decimal128.SinPi(x), + nameof(Decimal128.CosPi) => Decimal128.CosPi(x), + nameof(Decimal128.TanPi) => Decimal128.TanPi(x), + nameof(Decimal128.Exp2) => Decimal128.Exp2(x), + nameof(Decimal128.Exp2M1) => Decimal128.Exp2M1(x), + _ => throw new InvalidOperationException($"Unexpected operation '{operation}'."), + }; + + DecimalIeee754IntelTestData.AssertResultWithinUlp( + Unsafe.BitCast(actual), + Unsafe.BitCast(expected), recordedUlp: 0, limit: 1); + + if (operation == nameof(Decimal128.SinPi)) + { + Assert.Equal(actual, Decimal128.SinCosPi(x).SinPi); + Assert.Equal(-actual, Decimal128.SinPi(-x)); + } + else if (operation == nameof(Decimal128.CosPi)) + { + Assert.Equal(actual, Decimal128.SinCosPi(x).CosPi); + Assert.Equal(actual, Decimal128.CosPi(-x)); + } + } + + [Theory] + [InlineData("-1", "0.5", "-0.5")] + [InlineData("1", "2", "1")] + [InlineData("10", "1024", "1023")] + public static void Exp2IntegerAccuracyTest(string input, string expected, string expectedM1) + { + Decimal128 x = Decimal128.Parse(input, CultureInfo.InvariantCulture); + Assert.Equal(Decimal128.Parse(expected, CultureInfo.InvariantCulture), Decimal128.Exp2(x)); + Assert.Equal(Decimal128.Parse(expectedM1, CultureInfo.InvariantCulture), Decimal128.Exp2M1(x)); + } + + [Theory] + [InlineData("-14145.68261823660227840883486429809", "4.05598317787840255697517496726153e-6144")] + public static void ExpSubnormalRoundingTest(string input, string expected) + { + Assert.Equal(Decimal128.Parse(expected, CultureInfo.InvariantCulture), + Decimal128.Exp(Decimal128.Parse(input, CultureInfo.InvariantCulture))); + } + + [Theory] + [InlineData(-6177, 0)] + [InlineData(-6176, 1)] + [InlineData(-6175, 10)] + public static void Exp10SubnormalAccuracyTest(int input, int expectedUnits) + { + Assert.Equal(Decimal128.Epsilon * expectedUnits, Decimal128.Exp10(input)); + } + + [Theory] + [InlineData("0.999999999999999999999999999999999", "1e33", "0.367879441171442321595523770161460683506091")] + [InlineData("1.000000000000000000000000000000001", "1e33", "2.71828182845904523536028747135266113861633")] + public static void PowNearOneAccuracyTest(string input, string exponent, string oracle) + { + Decimal128 actual = Decimal128.Pow(Decimal128.Parse(input, CultureInfo.InvariantCulture), + Decimal128.Parse(exponent, CultureInfo.InvariantCulture)); + Decimal128 expected = Decimal128.Parse(oracle, CultureInfo.InvariantCulture); + DecimalIeee754IntelTestData.AssertResultWithinUlp( + Unsafe.BitCast(actual), + Unsafe.BitCast(expected), recordedUlp: 0, limit: 1); + } + + [Theory] + [InlineData("0.9999999999999999999999999999999999", 7, "0.999999999999999999999999999999999985714285714285714")] + [InlineData("1.000000000000000000000000000000001", -7, "0.999999999999999999999999999999999857142857142857143")] + [InlineData("0.7500000000000000000000000000000001", 7, "0.959735609788702600120982840125601770360737133179383")] + [InlineData("1.300000000000000000000000000000001", -7, "0.963213095002933238402631258970973215179971809330741")] + public static void RootNDecimalAccuracyTest(string input, int n, string oracle) + { + Decimal128 x = Decimal128.Parse(input, CultureInfo.InvariantCulture); + Decimal128 expected = Decimal128.Parse(oracle, CultureInfo.InvariantCulture); + Decimal128 actual = Decimal128.RootN(x, n); + DecimalIeee754IntelTestData.AssertResultWithinUlp( + Unsafe.BitCast(actual), + Unsafe.BitCast(expected), recordedUlp: 0, limit: 1); + Assert.Equal(-actual, Decimal128.RootN(-x, n)); + } + public static IEnumerable Parse_Valid_TestData() { NumberStyles defaultStyle = NumberStyles.Number; @@ -1999,6 +2282,7 @@ public static void LogAccuracyTest(double input) [InlineData(0xFC00000000000000UL, 0x0000000000000000UL, 0x3040000000000000UL, 0x0000000000000002UL, 0xFC00000000000000UL, 0x0000000000000000UL)] // log(NaN, 2) = NaN [InlineData(0x3040000000000000UL, 0x0000000000000002UL, 0xFC00000000000000UL, 0x0000000000000000UL, 0xFC00000000000000UL, 0x0000000000000000UL)] // log(2, NaN) = NaN [InlineData(0x3040000000000000UL, 0x0000000000000002UL, 0x3040000000000000UL, 0x0000000000000001UL, 0x7C00000000000000UL, 0x0000000000000000UL)] // log(2, 1) = NaN (base 1) + [InlineData(0x303A000000000000UL, 0x00000000000003E8UL, 0x303E000000000000UL, 0x0000000000000005UL, 0xB040000000000000UL, 0x0000000000000000UL)] // log(1.000, 0.5) = -0 public static void LogNewBaseTest(ulong valueUpper, ulong valueLower, ulong baseUpper, ulong baseLower, ulong expectedUpper, ulong expectedLower) { Decimal128 result = Decimal128.Log(Unsafe.BitCast(new UInt128(valueUpper, valueLower)), Unsafe.BitCast(new UInt128(baseUpper, baseLower))); @@ -2006,14 +2290,18 @@ public static void LogNewBaseTest(ulong valueUpper, ulong valueLower, ulong base } [Theory] - [InlineData(8.0, 2.0)] - [InlineData(100.0, 10.0)] - [InlineData(2.5, 3.0)] - public static void LogNewBaseAccuracyTest(double input, double newBase) + [InlineData("8", "2", "3.000000000000000000000000000000000")] + [InlineData("100", "10", "2.000000000000000000000000000000000")] + [InlineData("2.5", "3", "0.8340437671464697300975132933358795420083467265341692611822354509568071")] + [InlineData("1.000000000000000000000000000000001", "0.9999999999999999999999999999999999", "-9.999999999999999999999999999999994500000000000000000000000000000003575")] + public static void LogNewBaseAccuracyTest(string input, string newBase, string oracle) { - double expected = double.Log(input, newBase); - double actual = (double)Decimal128.Log((Decimal128)input, (Decimal128)newBase); - Assert.True(double.Abs(actual - expected) <= 1e-13 * double.Abs(double.MaxMagnitude(expected, 1.0)), $"log({input}, {newBase}): expected {expected}, got {actual}"); + Decimal128 actual = Decimal128.Log(Decimal128.Parse(input, CultureInfo.InvariantCulture), + Decimal128.Parse(newBase, CultureInfo.InvariantCulture)); + Decimal128 expected = Decimal128.Parse(oracle, CultureInfo.InvariantCulture); + DecimalIeee754IntelTestData.AssertResultWithinUlp( + Unsafe.BitCast(actual), + Unsafe.BitCast(expected), recordedUlp: 0, limit: 1); } [Theory] @@ -2310,6 +2598,9 @@ public static void HypotAccuracyTest(double x, double y) [InlineData(0x3040000000000000UL, 0x0000000000000000UL, -5, 0x7800000000000000UL, 0x0000000000000000UL)] // rootn(+0, n < 0) = +Infinity [InlineData(0xB040000000000000UL, 0x0000000000000000UL, -5, 0xF800000000000000UL, 0x0000000000000000UL)] // rootn(-0, odd < 0) = -Infinity [InlineData(0xB040000000000000UL, 0x0000000000000004UL, 2, 0x7C00000000000000UL, 0x0000000000000000UL)] // rootn(-4, even) = NaN + [InlineData(0x3038000000000000UL, 0x0000000000001B58UL, 1, 0x2FFD59206BDFDF06UL, 0x8D497D4600000000UL)] // rootn(0.7000, 1) uses the full-precision cohort + [InlineData(0xB038000000000000UL, 0x0000000000001B58UL, 1, 0xAFFD59206BDFDF06UL, 0x8D497D4600000000UL)] // rootn(-0.7000, 1) uses the full-precision cohort + [InlineData(0x0040000000000000UL, 0x0000000000000001UL, 1, 0x000004EE2D6D415BUL, 0x85ACEF8100000000UL)] // subnormal padding stops at the minimum quantum public static void RootNTest(ulong valueUpper, ulong valueLower, int n, ulong expectedUpper, ulong expectedLower) { Decimal128 result = Decimal128.RootN(Unsafe.BitCast(new UInt128(valueUpper, valueLower)), n); @@ -2715,7 +3006,7 @@ public static void SinPiAccuracyTest(string input, string oracle, double ulpLimi { // The engine evaluates in software binary128 (as Intel does), so the result is compared to a // high-precision oracle -- the true value rounded to Decimal128 by the independently tested parser -- - // in decimal ULPs. Exact identities use a 0 ULP limit; near-singular arguments a documented wider one. + // in decimal ULPs. Exact identities use a 0 ULP limit. Decimal128 actual = Decimal128.SinPi(Decimal128.Parse(input, CultureInfo.InvariantCulture)); Decimal128 expected = Decimal128.Parse(oracle, CultureInfo.InvariantCulture); DecimalIeee754IntelTestData.AssertResultWithinUlp( @@ -2761,7 +3052,7 @@ public static void CosPiHalfIntegerReturnsPositiveZero(double input) [InlineData("2.25", "0.707106781186547524400844362104849039284835938", 2.0)] [InlineData("0.1", "0.951056516295153572116439333379382143405698634", 2.0)] [InlineData("1234.567", "-0.208935890402411702274907259384464393664923236", 2.0)] - [InlineData("0.4999999", "0.000000314159265358974156133484288383422682765979151", 32.0)] // near a zero -> cancellation + [InlineData("0.4999999", "0.000000314159265358974156133484288383422682765979151", 1.0)] [InlineData("1", "-1.00000000000000000000000000000000000000000000", 0.0)] // cosPi(odd integer) = -1 exactly [InlineData("2", "1.00000000000000000000000000000000000000000000", 0.0)] // cosPi(even integer) = 1 exactly [InlineData("0.5", "0.0", 0.0)] // cosPi(half-integer) is an exact zero @@ -2944,7 +3235,7 @@ public static void AcosPiZeroTest() [InlineData("0.25", "0.419569376744833756229049806671515744415935569", 2.0)] [InlineData("-0.5", "0.666666666666666666666666666666666666666666667", 2.0)] [InlineData("0.999", "0.0142364374062396550708063524160105321570871313", 2.0)] - [InlineData("0.9999999", "0.000142352509869706344081743805264037099381195810", 32.0)] // near 1 -> cancellation + [InlineData("0.9999999", "0.000142352509869706344081743805264037099381195810", 1.0)] [InlineData("0.5", "0.333333333333333333333333333333333333333333333", 2.0)] [InlineData("0", "0.500000000000000000000000000000000000000000000", 0.0)] // acosPi(0) = 1/2 [InlineData("1", "0.0", 0.0)] // acosPi(1) = 0 diff --git a/src/libraries/System.Runtime/tests/System.Runtime.Tests/System/Decimal32Tests.cs b/src/libraries/System.Runtime/tests/System.Runtime.Tests/System/Decimal32Tests.cs index c6a19ae9d67f73..60e038c5e058d3 100644 --- a/src/libraries/System.Runtime/tests/System.Runtime.Tests/System/Decimal32Tests.cs +++ b/src/libraries/System.Runtime/tests/System.Runtime.Tests/System/Decimal32Tests.cs @@ -13,6 +13,35 @@ namespace System.Tests { public class Decimal32Tests { + [Theory] + [InlineData("-219.1358", "6.76911e-96")] + [InlineData("-219.6328", "4.11801e-96")] + public static void ExpSubnormalRoundingTest(string input, string expected) + { + Assert.Equal(Decimal32.Parse(expected, CultureInfo.InvariantCulture), + Decimal32.Exp(Decimal32.Parse(input, CultureInfo.InvariantCulture))); + } + + [Theory] + [InlineData(-102, 0)] + [InlineData(-101, 1)] + [InlineData(-100, 10)] + public static void Exp10SubnormalAccuracyTest(int input, int expectedUnits) + { + Assert.Equal(Decimal32.Epsilon * (Decimal32)expectedUnits, Decimal32.Exp10((Decimal32)input)); + } + + [Theory] + [InlineData("0.9999999", "0.000447213599226738", "-1.00000005e-7", "3.14159265358974e-7")] + public static void TranscendentalBoundaryAccuracyTest(string input, string acos, string log, string sinPi) + { + Decimal32 x = Decimal32.Parse(input, CultureInfo.InvariantCulture); + Assert.Equal(Decimal32.Parse(acos, CultureInfo.InvariantCulture), Decimal32.Acos(x)); + Assert.Equal(Decimal32.Parse(log, CultureInfo.InvariantCulture), Decimal32.Log(x)); + Assert.Equal(Decimal32.Parse(sinPi, CultureInfo.InvariantCulture), Decimal32.SinPi(x)); + Assert.Equal(Decimal32.SinPi(x), Decimal32.SinCosPi(x).SinPi); + } + public static IEnumerable Parse_Valid_TestData() { NumberStyles defaultStyle = NumberStyles.Number; @@ -2008,20 +2037,25 @@ public static void LogAccuracyTest(double input) [InlineData(0x7C000000U, 0x32800002U, 0x7C000000U)] // log(NaN, 2) = NaN [InlineData(0x32800002U, 0x7C000000U, 0x7C000000U)] // log(2, NaN) = NaN [InlineData(0x32800002U, 0x32800001U, 0x7C000000U)] // log(2, 1) = NaN (base 1) + [InlineData(0x310003E8U, 0x32000005U, 0xB2800000U)] // log(1.000, 0.5) = -0 public static void LogNewBaseTest(uint value, uint newBase, uint expected) { Assert.Equal(expected, Unsafe.BitCast(Decimal32.Log(Unsafe.BitCast(value), Unsafe.BitCast(newBase)))); } [Theory] - [InlineData(8.0, 2.0)] - [InlineData(100.0, 10.0)] - [InlineData(2.5, 3.0)] - public static void LogNewBaseAccuracyTest(double input, double newBase) + [InlineData("8", "2", "3.000000")] + [InlineData("100", "10", "2.000000")] + [InlineData("2.5", "3", "0.8340437671464697300975132933358795420083467265341692611822354509568071")] + [InlineData("1.000001", "0.9999999", "-9.99999450000357499733708545573573529")] + public static void LogNewBaseAccuracyTest(string input, string newBase, string oracle) { - double expected = double.Log(input, newBase); - double actual = (double)Decimal32.Log((Decimal32)input, (Decimal32)newBase); - Assert.True(double.Abs(actual - expected) <= 5e-7 * double.Abs(double.MaxMagnitude(expected, 1.0)), $"log({input}, {newBase}): expected {expected}, got {actual}"); + Decimal32 actual = Decimal32.Log(Decimal32.Parse(input, CultureInfo.InvariantCulture), + Decimal32.Parse(newBase, CultureInfo.InvariantCulture)); + Decimal32 expected = Decimal32.Parse(oracle, CultureInfo.InvariantCulture); + DecimalIeee754IntelTestData.AssertResultWithinUlp( + Unsafe.BitCast(actual), + Unsafe.BitCast(expected), recordedUlp: 0, limit: 1); } [Theory] @@ -2309,6 +2343,9 @@ public static void HypotAccuracyTest(double x, double y) [InlineData(0x32800000U, -5, 0x78000000U)] // rootn(+0, n < 0) = +Infinity [InlineData(0xB2800000U, -5, 0xF8000000U)] // rootn(-0, odd < 0) = -Infinity [InlineData(0xB2800004U, 2, 0x7C000000U)] // rootn(-4, even) = NaN + [InlineData(0x30801B58U, 1, 0x2F6ACFC0U)] // rootn(0.7000, 1) uses the full-precision cohort + [InlineData(0xB0801B58U, 1, 0xAF6ACFC0U)] // rootn(-0.7000, 1) uses the full-precision cohort + [InlineData(0x02800001U, 1, 0x000186A0U)] // subnormal padding stops at the minimum quantum public static void RootNTest(uint value, int n, uint expected) { Assert.Equal(expected, Unsafe.BitCast(Decimal32.RootN(Unsafe.BitCast(value), n))); diff --git a/src/libraries/System.Runtime/tests/System.Runtime.Tests/System/Decimal64Tests.cs b/src/libraries/System.Runtime/tests/System.Runtime.Tests/System/Decimal64Tests.cs index 8bd8bb011da8e8..fce5fc6eb0fec3 100644 --- a/src/libraries/System.Runtime/tests/System.Runtime.Tests/System/Decimal64Tests.cs +++ b/src/libraries/System.Runtime/tests/System.Runtime.Tests/System/Decimal64Tests.cs @@ -12,6 +12,34 @@ namespace System.Tests { public class Decimal64Tests { + [Theory] + [InlineData("-883.2925549776085", "2.45990008214473e-384")] + public static void ExpSubnormalRoundingTest(string input, string expected) + { + Assert.Equal(Decimal64.Parse(expected, CultureInfo.InvariantCulture), + Decimal64.Exp(Decimal64.Parse(input, CultureInfo.InvariantCulture))); + } + + [Theory] + [InlineData(-399, 0)] + [InlineData(-398, 1)] + [InlineData(-397, 10)] + public static void Exp10SubnormalAccuracyTest(int input, int expectedUnits) + { + Assert.Equal(Decimal64.Epsilon * expectedUnits, Decimal64.Exp10(input)); + } + + [Theory] + [InlineData("0.9999999999999999", "1.4142135623730950605868e-8", "-1.00000000000000005e-16", "3.14159265358979323846264e-16")] + public static void TranscendentalBoundaryAccuracyTest(string input, string acos, string log, string sinPi) + { + Decimal64 x = Decimal64.Parse(input, CultureInfo.InvariantCulture); + Assert.Equal(Decimal64.Parse(acos, CultureInfo.InvariantCulture), Decimal64.Acos(x)); + Assert.Equal(Decimal64.Parse(log, CultureInfo.InvariantCulture), Decimal64.Log(x)); + Assert.Equal(Decimal64.Parse(sinPi, CultureInfo.InvariantCulture), Decimal64.SinPi(x)); + Assert.Equal(Decimal64.SinPi(x), Decimal64.SinCosPi(x).SinPi); + } + public static IEnumerable Parse_Valid_TestData() { NumberStyles defaultStyle = NumberStyles.Number; @@ -1996,20 +2024,25 @@ public static void LogAccuracyTest(double input) [InlineData(0x7C00000000000000UL, 0x31C0000000000002UL, 0x7C00000000000000UL)] // log(NaN, 2) = NaN [InlineData(0x31C0000000000002UL, 0x7C00000000000000UL, 0x7C00000000000000UL)] // log(2, NaN) = NaN [InlineData(0x31C0000000000002UL, 0x31C0000000000001UL, 0x7C00000000000000UL)] // log(2, 1) = NaN (base 1) + [InlineData(0x31600000000003E8UL, 0x31A0000000000005UL, 0xB1C0000000000000UL)] // log(1.000, 0.5) = -0 public static void LogNewBaseTest(ulong value, ulong newBase, ulong expected) { Assert.Equal(expected, Unsafe.BitCast(Decimal64.Log(Unsafe.BitCast(value), Unsafe.BitCast(newBase)))); } [Theory] - [InlineData(8.0, 2.0)] - [InlineData(100.0, 10.0)] - [InlineData(2.5, 3.0)] - public static void LogNewBaseAccuracyTest(double input, double newBase) + [InlineData("8", "2", "3.000000000000000")] + [InlineData("100", "10", "2.000000000000000")] + [InlineData("2.5", "3", "0.8340437671464697300975132933358795420083467265341692611822354509568071")] + [InlineData("1.000000000000001", "0.9999999999999999", "-9.999999999999994500000000000003575")] + public static void LogNewBaseAccuracyTest(string input, string newBase, string oracle) { - double expected = double.Log(input, newBase); - double actual = (double)Decimal64.Log((Decimal64)input, (Decimal64)newBase); - Assert.True(double.Abs(actual - expected) <= 1e-13 * double.Abs(double.MaxMagnitude(expected, 1.0)), $"log({input}, {newBase}): expected {expected}, got {actual}"); + Decimal64 actual = Decimal64.Log(Decimal64.Parse(input, CultureInfo.InvariantCulture), + Decimal64.Parse(newBase, CultureInfo.InvariantCulture)); + Decimal64 expected = Decimal64.Parse(oracle, CultureInfo.InvariantCulture); + DecimalIeee754IntelTestData.AssertResultWithinUlp( + Unsafe.BitCast(actual), + Unsafe.BitCast(expected), recordedUlp: 0, limit: 1); } [Theory] @@ -2297,6 +2330,9 @@ public static void HypotAccuracyTest(double x, double y) [InlineData(0x31C0000000000000UL, -5, 0x7800000000000000UL)] // rootn(+0, n < 0) = +Infinity [InlineData(0xB1C0000000000000UL, -5, 0xF800000000000000UL)] // rootn(-0, odd < 0) = -Infinity [InlineData(0xB1C0000000000004UL, 2, 0x7C00000000000000UL)] // rootn(-4, even) = NaN + [InlineData(0x3140000000001B58UL, 1, 0x2FD8DE76816D8000UL)] // rootn(0.7000, 1) uses the full-precision cohort + [InlineData(0xB140000000001B58UL, 1, 0xAFD8DE76816D8000UL)] // rootn(-0.7000, 1) uses the full-precision cohort + [InlineData(0x01C0000000000001UL, 1, 0x00005AF3107A4000UL)] // subnormal padding stops at the minimum quantum public static void RootNTest(ulong value, int n, ulong expected) { Assert.Equal(expected, Unsafe.BitCast(Decimal64.RootN(Unsafe.BitCast(value), n)));