From 2026fcde19c8f431ef865a2437bb8941fc387a25 Mon Sep 17 00:00:00 2001 From: Tanner Gooding Date: Tue, 15 Sep 2026 06:33:17 -0700 Subject: [PATCH] Improve IEEE decimal transcendental accuracy (#133921) Fixes https://github.com/dotnet/runtime/issues/133432. Preserve decimal precision through cancellation-sensitive argument reduction for `Decimal32`, `Decimal64`, and `Decimal128`, without importing Intel DFP's large binary/decimal conversion tables. Remove redundant work in the affected transcendental paths. - Form small endpoint residuals in decimal before converting them to the existing binary engine. This covers inverse trig, inverse hyperbolic functions, logarithms, and the compensated `Pow` path used by `RootN`. - Verify and restore exact dyadic operands against their original decimal coefficients, including alternate representations of `1`. Skip residual correction where it is unnecessary. - Keep both logarithms and their quotient in working precision for `Log(x, base)`, rounding to decimal only once. Preserve the existing special-value handling. - Reduce Pi-scaled trig arguments in decimal, preserving distances from integers and half-integers. Evaluate the already-reduced angle directly; share polynomial evaluation for `SinCosPi` and `TanPi`, and bypass unnecessary computation for exact results. - Reuse the existing accurate `ln(2)` constant in `Exp2` and `Exp2M1`; the previous reconstruction subtracted a low part that was already negative. - Select the final decimal quantum before rounding, avoiding a second rounding for subnormal results. - Reuse repeated power-of-ten chunk multipliers, classify `Pow` exponent parity only for negative bases, and convert `RootN` degrees directly to working precision. Handle `RootN(x, 1)` without log/exp evaluation while retaining the existing result cohort. The existing approximation algorithms and compact tables are retained. No public API is added. Rounded numerical results intentionally change; source and binary compatibility are unchanged. ## Accuracy: before / after Observed maximum absolute error against a high-precision mathematical reference, measured in decimal ULPs. **Bold after-values highlight improvements visible at two decimal places.** `0.50 / 0.50` is deliberately not highlighted; rounding to two decimals can hide small differences on either side of the halfway point. The baseline is `22b24296b6a2580cab3fcc89489bbda3eea79a23`. Both builds ran the same inputs in two seeded sweeps, including full-significand random values, endpoint probes, and exponent windows. The sweeps evaluated 211,483 finite mathematical results. References used 220 decimal digits, with each of the 360 post-change window maxima per sweep rechecked at 440 digits. The subsequent wide-log change was evaluated over the same corpora; only log-base results changed. Later performance-only changes were checked for bit-identical results separately. These are sample maxima, not proven bounds or a universal correct-rounding guarantee. The tables exclude exact-zero references and cases where the baseline returned a nonfinite result for a finite reference; those corrections are described separately below. ULP spacing is determined from the mathematical result and format precision, clamped to the minimum subnormal quantum. Original operands are parsed to the target format before computing the reference. | API | Decimal32 before / after | Decimal64 before / after | Decimal128 before / after | |---|---:|---:|---:| | `Asin` | 0.50 / 0.50 | 0.50 / 0.50 | 5.18e+11 / **0.50** | | `Acos` | 0.50 / 0.50 | 0.50 / 0.50 | 5.18e+28 / **0.50** | | `AsinPi` | 0.50 / 0.50 | 0.55 / **0.50** | 1.65e+12 / **0.50** | | `AcosPi` | 0.50 / 0.50 | 0.50 / 0.50 | 1.65e+29 / **0.50** | | `Acosh` | 0.50 / 0.50 | 0.50 / 0.50 | 5.39e+27 / **0.50** | | `Atanh` | 0.50 / 0.50 | 0.50 / 0.50 | 2.59e+27 / **0.50** | | `Log` | 0.50 / 0.50 | 0.50 / 0.50 | 1.14e+29 / **0.50** | | `Log2` | 0.50 / 0.50 | 0.50 / 0.50 | 1.50e+29 / **0.50** | | `Log10` | 0.50 / 0.50 | 0.50 / 0.50 | 4.50e+29 / **0.50** | | `LogP1` | 0.50 / 0.50 | 0.50 / 0.50 | 5.18e+27 / **0.50** | | `Log2P1` | 0.50 / 0.50 | 0.50 / 0.50 | 7.48e+26 / **0.50** | | `Log10P1` | 0.50 / 0.50 | 0.50 / 0.50 | 2.25e+27 / **0.50** | | `Log(x, base)` | 2.86 / **0.50** | 3.83 / **0.50** | 5.29e+29 / **0.50** | | `SinPi` | 0.50 / 0.50 | 0.50 / 0.50 | 1.63e+29 / **0.50** | | `CosPi` | 0.50 / 0.50 | 0.50 / 0.50 | 1.63e+29 / **0.50** | | `TanPi` | 0.50 / 0.50 | 0.50 / 0.50 | 2.20e+29 / **0.50** | | `SinCosPi` sine component | 0.50 / 0.50 | 0.50 / 0.50 | 1.63e+29 / **0.50** | | `SinCosPi` cosine component | 0.50 / 0.50 | 0.50 / 0.50 | 1.63e+29 / **0.50** | | `Exp` | 0.53 / **0.50** | 0.50 / 0.50 | 1.19 / 1.19 | | `ExpM1` | 0.50 / 0.50 | 0.50 / 0.50 | 1.31 / 1.31 | | `Exp2` | 0.53 / **0.50** | 0.74 / **0.50** | 4.68e+18 / **2.84** | | `Exp2M1` | 0.50 / 0.50 | 0.76 / **0.50** | 4.68e+18 / **2.84** | | `Exp10` | 0.54 / **0.50** | 0.52 / **0.50** | 1.54 / 1.54 | | `Exp10M1` | 0.50 / 0.50 | 0.50 / 0.50 | 1.54 / 1.54 | | `Pow` | 0.54 / **0.50** | 0.52 / **0.50** | 8.43e+32 / **2.37** | | `RootN` | 0.50 / 0.50 | 0.50 / 0.50 | 0.56 / 0.56 | The extra accuracy headroom of the shared 128-bit working significand explains why most Decimal32/64 results were already good. Their improvements include subnormal rounding, the base-two exponential constant at large exponents, avoiding rounded log-base intermediates, and exact-result handling. ### Cancellation-sensitive ranges: Decimal128 For endpoint probes, `delta = m * 10^-k`, with `m` in `{1, 2, 5, 9}` and `k` spanning the representable decimal precision. Interior/tiny/large-input windows not listed here remain approximately `0.50 / 0.50` for the inverse functions and single-argument logarithms. | API | Sampled range | Before | After | |---|---|---:|---:| | `Asin` | `x = +/- (1 - delta)` | 5.18e+11 | **0.49** | | `Acos` | `x = +/- (1 - delta)` | 5.18e+28 | **0.49** | | `AsinPi` | `x = +/- (1 - delta)` | 1.65e+12 | **0.49** | | `AcosPi` | `x = +/- (1 - delta)` | 1.65e+29 | **0.50** | | `Atanh` | `x = +/- (1 - delta)` | 2.59e+27 | **0.50** | | `Acosh` | `x = 1 + delta` | 5.39e+27 | **0.50** | | `Log` | `x = 1 +/- delta` | 1.14e+29 | **0.50** | | `Log2` | `x = 1 +/- delta` | 1.50e+29 | **0.50** | | `Log10` | `x = 1 +/- delta` | 4.50e+29 | **0.50** | | `LogP1` | `x = -1 + delta` | 5.18e+27 | **0.50** | | `Log2P1` | `x = -1 + delta` | 7.48e+26 | **0.49** | | `Log10P1` | `x = -1 + delta` | 2.25e+27 | **0.38** | | `Pow` | Near-one bases; exponents targeting ordinary, subnormal, and large finite results | 8.43e+32 | **2.37** | | `Pow` | Positive bases in `0.1..10`, same output-window sampling | 8.45 | **1.47** | `Log(x, base)` samples used bases `0.5`, `2`, `10`, and the adjacent representable values below/above `1`, over tiny, near-one, ordinary, and large positive arguments. Keeping both logarithms and the division wide reduces the sampled maxima to approximately `0.49999643`, `0.50000000`, and `0.50002430` ULP for Decimal32/64/128 respectively. The Decimal128 result is slightly above half an ULP; the table's `0.50` is not a correct-rounding claim. ### Pi reduction ranges: Decimal128 | API | Sampled range | Before | After | |---|---|---:|---:| | `SinPi` | Ordinary values in `-1000..1000` | 1.65 | **0.50** | | `CosPi` | Ordinary values in `-1000..1000` | 1.26 | **0.50** | | `TanPi` | Ordinary values in `-1000..1000` | 7.51 | **0.50** | | `SinPi`, `CosPi` | Beside integers/half-integers, down to the representable spacing | 1.63e+29 | **0.48** | | `TanPi` | Beside integers/half-integers, excluding exact poles | 1.65e+29 | **0.49** | | `SinPi`, `CosPi` | Full-precision fractional values with decimal exponents 3, 17, and 32 | 1.60e+29 | **0.50** | | `TanPi` | Same large fractional inputs | 2.20e+29 | **0.50** | The corresponding `SinCosPi` components have the same errors. Tiny-input maxima remain approximately `0.50 / 0.50`; the actual largest post-change sine/tangent error is slightly above half an ULP, so this is not a correct-rounding guarantee. ### Base-two exponential ranges: Decimal128 Here `u = x * log10(2)`, the decimal exponent of the mathematical `2^x` before subtracting one. Each row takes the maximum across `Exp2` and `Exp2M1`. | Sampled range | Before | After | |---|---:|---:| | Tiny inputs and `abs(x) < 1` | 4.55e+14 | **0.50** | | Subnormal exponential outputs, `u` in `-6177..-6143` | 3.62e+17 | **0.50** | | Negative normal exponential outputs, `u` in `-6143..-1` | 4.28e+18 | **2.48** | | `u` in `0..100` | 6.95e+16 | **0.51** | | `u` in `100..1000` | 7.18e+17 | **0.69** | | `u` in `1000..6144` | 4.47e+18 | **2.57** | | Last finite decade, `u` in `6144..6145` | 4.68e+18 | **2.56** | | Additional probes immediately beside the overflow boundary | 4.68e+18 | **2.84** | The last row excludes cases where the baseline overflowed despite a finite mathematical result; the new build returns finite results for those probes. ### Smaller-format and subnormal improvements | Format / API | Sampled range | Before | After | |---|---|---:|---:| | Decimal32 `Exp`, `Exp2` | Subnormal exponential outputs | 0.53 | **0.50** | | Decimal32 `Exp10` | Subnormal exponential outputs | 0.54 | **0.50** | | Decimal64 `AsinPi` | Tiny inputs spanning the negative exponent range | 0.55 | **0.50** | | Decimal64 `Exp2` | Large negative inputs producing normal results | 0.70 | **0.50** | | Decimal64 `Exp2`, `Exp2M1` | Last finite output decade | 0.76 | **0.50** | | Decimal64 `Exp10` | Subnormal exponential outputs | 0.52 | **0.50** | | Decimal128 `Exp10` | Subnormal exponential outputs | 0.52 | **0.50** | The remaining exponential families retain approximately the same overall maxima: Decimal128 `Exp`/`ExpM1` reach `1.19`/`1.31`, and `Exp10`/`Exp10M1` reach `1.54` at extreme exponents. Fixing final subnormal rounding does not remove all conversion/approximation error. ### Exact results and remaining limitations - Exact Pi zeros are now preserved for alternate decimal representations in all three formats. For example, Decimal32 `SinPi(1000.000)` previously returned approximately `9.45e-36`; it now returns zero. These are reported as exact-result corrections, not enormous ULP figures computed using the minimum subnormal spacing at zero. - Alternate representations of `1` retain the correct inverse-function endpoints and poles. For example, Decimal64/128 `Acosh(1.0000000)` no longer returns NaN. - The original issue's Decimal128 `Acos(1 - 1e-30)` improves from `2.24e+24` to `0.07` ULP, and `Asin(1 - 1e-30)` from `2.24e+09` to `0.36` ULP. - This does not promise universal correct rounding or a universal two-ULP bound. Extreme-exponent errors remain, as shown above. - A separate pre-existing cohort-dependent overflow issue remains: Decimal128 `Exp10` and `Exp10M1` of `6145.000000000000000000000000000000` return `MaxValue` instead of infinity, whereas the representation `6145` returns infinity. This was reproduced on both the baseline and this change and is outside this PR. ## Performance and size Local BenchmarkDotNet measurements used Windows x64 / Ryzen 9 7950X. The initial accuracy fixes, compared with the original Release baseline, generally had roughly 0-5% overhead for exact inputs and 7-22% for inexact inverse/hyperbolic and `LogP1` cases requiring decimal residual preservation. These are not claimed to be performance-neutral; the later optimizations recover work in specific paths rather than establishing a uniform end-to-end speedup. Incremental measurements of the later optimizations showed approximately 8-15% faster ordinary `SinCosPi` and 10-23% faster `TanPi` from paired evaluation. Wide log-base evaluation improved sampled Decimal128 cases by roughly 11-23%. These comparisons use the preceding implementation, not the original pre-PR baseline. The exact-result shortcuts showed the following incremental timings (ns/op, before / after), with two launches and fixed processor affinity: | Case | Decimal32 | Decimal64 | Decimal128 | |---|---:|---:|---:| | `RootN(0.7, 1)` | 384.36 / 6.01 | 389.04 / 6.41 | 388.54 / 8.74 | | `TanPi(1000.000)` | 38.23 / 7.84 | 38.03 / 8.21 | 39.36 / 9.61 | The measurements are mixed outside the targeted paths: ordinary `RootN` changes were small, and half-integer `TanPi` was about 1 ns slower for Decimal32/64. Earlier exact-result optimizations also showed small ordinary sine/cosine regressions. No universal performance improvement is claimed. No managed allocations were observed. The final measured Release CoreLib file is `18,305,024` bytes, versus `18,300,928` bytes for the original baseline (+4,096 bytes). No large lookup tables were added. ## Validation - Checked and Release CoreLib builds, including refreshed library testhosts. - Full `System.Runtime.Tests` passed after the product changes and test consolidation. After removing three redundant test rows, the final Decimal32/64/128 test run passed 7,407 tests, with 69 skipped and zero failures. - Strengthened the existing log-base accuracy tests with decimal oracles rather than adding parallel test methods. The existing `Log(2.5, 3)` case now detects the old two-ULP error in every format; all three fail on the pre-wide-log baseline and pass with this change. - Focused raw-bit tests cover signed-zero log-base results and the normal/subnormal result cohorts of the `RootN(x, 1)` shortcut. Existing endpoint, subnormal, Pi-cohort, and degree-boundary coverage is retained. - 211,483 finite-reference sampling evaluations across the two range sweeps described above; subsequent performance-only simplifications preserved all bits in a 224,518-evaluation comparison against the preceding optimized build. - 4,924 Intel DFP Update 3 vectors checked during the initial accuracy work; all 1,386 sampled exact-operation vectors remained bit-identical. Execution and performance measurements were Windows x64 only. The changes are in shared CoreLib source. > [!NOTE] > This PR description was drafted with GitHub Copilot. --------- Co-authored-by: Copilot App <223556219+Copilot@users.noreply.github.com> --- ...mber.DecimalIeee754.DiyFp128Conversions.cs | 64 ++-- .../Number.DecimalIeee754.DiyFp128Exp.cs | 15 +- .../Number.DecimalIeee754.DiyFp128InvHyper.cs | 22 +- .../Number.DecimalIeee754.DiyFp128InvTrig.cs | 13 +- .../Number.DecimalIeee754.DiyFp128PiTrig.cs | 228 ++++++------- .../Number.DecimalIeee754.DiyFp128Pow.cs | 10 +- .../Number.DecimalIeee754.Transcendental.cs | 226 ++++++++++--- .../System/Decimal128Tests.cs | 311 +++++++++++++++++- .../System/Decimal32Tests.cs | 51 ++- .../System/Decimal64Tests.cs | 50 ++- 10 files changed, 722 insertions(+), 268 deletions(-) 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 64dd72c5cbd4b3..bd193b4276ffd1 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 @@ -12,6 +12,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; @@ -1993,20 +2022,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] @@ -2294,6 +2328,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)));