diff --git a/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs b/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs index ea5d214..0c7ac16 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs @@ -147,7 +147,14 @@ public readonly partial struct DoubleDouble : // the first product and subtraction gives |remainder.High| <= 2^457 // at the first correction, hence |correction| <= 2^908 and a finite sum // below 2^909 in magnitude. Retain normalization before the final add. - return PreciseMathHelper.NormalizeFinite(quotient, correction) + finalCorrection; + // The correlated correction product (including its exact fallback) has + // |High| <= 2^460, giving |finalRemainder| <= 2^462 and |finalCorrection| + // <= 2^912. Both final-add highs therefore pass the addition range guard. + // More tightly, the first product differs from left.High by < 4u*|left.High|, + // so |remainder.High| < 16u*|left.High| and |correction/quotient| < 32u + // (u=2^-53, including underflow). The normalized quotient cannot be zero. + DoubleDouble partial = PreciseMathHelper.NormalizeFinite(quotient, correction); + return PreciseMathHelper.AddScalarFinite(partial._high, partial._low, finalCorrection); } /// Applies the expansion operation without discarding the low component. @@ -231,8 +238,10 @@ public readonly partial struct DoubleDouble : double finalRemainder = PreciseMathHelper.SubtractDivisionCorrectionHigh(remainder, correctedProduct); double finalCorrection = finalRemainder / right; // The same finite bounds as DD/DD apply. Normalize the first correction - // before adding the last, preserving the DD/DD evaluation order. - return PreciseMathHelper.NormalizeFinite(quotient, correction) + finalCorrection; + // before adding the last, preserving the DD/DD evaluation order. Its high + // is nonzero and < 2^909; |finalCorrection| <= 2^912 fits scalar addition. + DoubleDouble partial = PreciseMathHelper.NormalizeFinite(quotient, correction); + return PreciseMathHelper.AddScalarFinite(partial._high, partial._low, finalCorrection); } /// Divides a scalar using a quotient estimate and two residual corrections. diff --git a/0-source/Just.PreciseMath/DoubleDouble.Comparison.cs b/0-source/Just.PreciseMath/DoubleDouble.Comparison.cs index 2d0a8d9..f93d97b 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.Comparison.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.Comparison.cs @@ -1,10 +1,18 @@ namespace Just.PreciseMath; -public readonly partial struct DoubleDouble : IComparable, IComparable, +public readonly partial struct DoubleDouble : IEquatable, IEqualityOperators, + IComparable, IComparable, IComparisonOperators { + /// + [Pure] + public bool Equals(DoubleDouble other) + { + return _high.Equals(other._high) && _low.Equals(other._low); + } + /// /// Tests numerical equality; NaN operands are never equal. /// diff --git a/0-source/Just.PreciseMath/DoubleDouble.RootFunctions.cs b/0-source/Just.PreciseMath/DoubleDouble.RootFunctions.cs index 067436f..606594c 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.RootFunctions.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.RootFunctions.cs @@ -285,13 +285,18 @@ public readonly partial struct DoubleDouble : IRootFunctions // where u=2^-53; both components of the input contribute to the residual. double residual = Math.FusedMultiplyAdd(-estimate, estimate, high) + low; double correction = residual / (2.0 * estimate); - DoubleDouble root = DoubleDouble.FromComponents(estimate, correction); + // Scaled |low| <= 2^-51 and the rounded sqrt seed's squared error + // is < 2^-49, so |correction| < 2^-48. Both inputs and their sum + // are finite (sum < 4); retain magnitude-aware normalization. + DoubleDouble root = PreciseMathHelper.NormalizeFinite(estimate, correction); // The root high stays normal, even if rounding reaches the next binade. // Normalize again to canonicalize a zero low, including any underflow // of an exceptionally sparse correction during rescaling. int rootExponent = exponent / 2; - root = DoubleDouble.FromComponents(Math.ScaleB(root.High, rootExponent), Math.ScaleB(root.Low, rootExponent)); + // rootExponent is in [-537, 511]. With root.High < 4 and normalized + // |root.Low| <= 2^-52, scaled components and their sum stay below 2^514. + root = PreciseMathHelper.NormalizeFinite(Math.ScaleB(root.High, rootExponent), Math.ScaleB(root.Low, rootExponent)); // In this branch |originalLow/originalHigh| < 2^-1020. The omitted // quadratic term is negligible even relative to the low correction. return separateLow ? root + RootLowCorrection(value, root, 2) : root; diff --git a/0-source/Just.PreciseMath/DoubleDouble.cs b/0-source/Just.PreciseMath/DoubleDouble.cs index 39ac13e..313b92f 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.cs @@ -81,12 +81,6 @@ public readonly partial struct DoubleDouble : public readonly double Low => _low; - /// - [Pure] - public bool Equals(DoubleDouble other) - { - return _high.Equals(other._high) && _low.Equals(other._low); - } /// [Pure] public override bool Equals(object? obj) diff --git a/0-source/Just.PreciseMath/PreciseMathHelper.cs b/0-source/Just.PreciseMath/PreciseMathHelper.cs index b31b79e..06e4aed 100644 --- a/0-source/Just.PreciseMath/PreciseMathHelper.cs +++ b/0-source/Just.PreciseMath/PreciseMathHelper.cs @@ -118,8 +118,11 @@ internal static class PreciseMathHelper double correction = remainder._high / right._high; double finalRemainder = SubtractDivisionCorrectionHigh(remainder, right * correction); double finalCorrection = finalRemainder / right._high; - // DD/DD's guard-derived finite bounds and normalization order apply. - return NormalizeFinite(quotient, correction) + finalCorrection; + // DD/DD's guard-derived bounds apply, including exact correction products: + // the normalized high is nonzero and < 2^909; |finalCorrection| <= 2^912. + // Keep normalization before scalar addition, without repeating its dispatch. + DoubleDouble partial = NormalizeFinite(quotient, correction); + return AddScalarFinite(partial._high, partial._low, finalCorrection); } // Only for division's second remainder: normalized finite inputs with highs @@ -149,7 +152,7 @@ internal static class PreciseMathHelper // Both components and their rounded sum must be finite. Unlike QuickTwoSum, // this entry point permits either magnitude order, including cancellation. - // AddFinite and DD division establish these bounds; this helper does not validate + // Arithmetic and Sqrt callers establish these bounds; this helper does not validate // them or canonicalize NaN/infinity. Use DoubleDouble.FromComponents for arbitrary pairs. // Retain the high zero's sign when low is zero, as FromComponents does. [MethodImpl(MethodImplOptions.AggressiveInlining)] @@ -181,6 +184,16 @@ internal static class PreciseMathHelper return AddScalarBoundary(high, low, value); } + return AddScalarFinite(high, low, value); + } + + // Normalized finite (high, low), with a negated zero low also allowed. + // high and value must pass the addition range guard and not both be zero. + // Public scalar addition checks these conditions; division establishes them + // before its final correction. Preserve all three transforms and grouping. + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static DoubleDouble AddScalarFinite(double high, double low, double value) + { (double sum, double error) = TwoAdd(high, value); // Specialize AddFinite for a zero-low scalar: adding that zero to low // has no residual. Retain the rounding error when combining the high diff --git a/1-tests/Just.PreciseMath.Tests/ArithmeticRangeTests.cs b/1-tests/Just.PreciseMath.Tests/ArithmeticRangeTests.cs index 8cf9653..7869245 100644 --- a/1-tests/Just.PreciseMath.Tests/ArithmeticRangeTests.cs +++ b/1-tests/Just.PreciseMath.Tests/ArithmeticRangeTests.cs @@ -141,7 +141,21 @@ public class ArithmeticRangeTests DoubleDouble finiteNormalized = PreciseMathHelper.NormalizeFinite(quotient, correction); BitConverter.DoubleToInt64Bits(finiteNormalized.High).ShouldBe(BitConverter.DoubleToInt64Bits(normalized.High)); BitConverter.DoubleToInt64Bits(finiteNormalized.Low).ShouldBe(BitConverter.DoubleToInt64Bits(normalized.Low)); - expected = normalized + finalCorrection; + // Pin the final-add caller invariants before removing public + // dispatch. Sparse correction products above can use the exact + // fallback; these bounds must cover their results too. + normalized.High.ShouldNotBe(0.0); + (Math.Abs(correction) < Math.Abs(quotient) * 0.5).ShouldBeTrue(); + (Math.Abs(normalized.High) < Math.ScaleB(1.0, 909)).ShouldBeTrue(); + (Math.Abs(remainder.High) < Math.ScaleB(1.0, 470)).ShouldBeTrue(); + (Math.Abs(finalCorrection) < Math.ScaleB(1.0, 921)).ShouldBeTrue(); + PreciseMathHelper.IsAdditionWithinFastRange(normalized.High).ShouldBeTrue(); + PreciseMathHelper.IsAdditionWithinFastRange(finalCorrection).ShouldBeTrue(); + (Math.ILogB(normalized.High) <= 1020).ShouldBeTrue(); + (Math.ILogB(finalCorrection) <= 1020).ShouldBeTrue(); + expected = PreviousFiniteScalarAdd(normalized.High, normalized.Low, finalCorrection); + AssertComponentBits(normalized + finalCorrection, expected); + AssertComponentBits(PreciseMathHelper.AddScalarFinite(normalized.High, normalized.Low, finalCorrection), expected); actual = left / right; } else if (operation == "*") @@ -231,6 +245,58 @@ public class ArithmeticRangeTests } } + [Fact] + public void ScalarAdditionPreservesThePreviousFiniteSequence() + { + // Characterization, not an accuracy oracle. Include cancellation, sparse + // lows, signed high zeros, and the negated zero low used by subtraction. + double[] highs = [0.0, -0.0, double.Epsilon, -double.Epsilon, + Math.ScaleB(1.0, -1022), 1.0, -1.0, Math.BitDecrement(2.0), + Math.BitDecrement(Math.ScaleB(1.0, 909)), + -Math.BitDecrement(Math.ScaleB(1.0, 909)), + Math.BitDecrement(Math.ScaleB(1.0, 921)), + -Math.BitDecrement(Math.ScaleB(1.0, 921)), + Math.BitDecrement(Math.ScaleB(1.0, 1021))]; + foreach (double high in highs) + { + foreach (double low in new[] { 0.0, -0.0, double.Epsilon, -double.Epsilon, + Math.ScaleB(high, -54), -Math.ScaleB(high, -54) }) + { + DoubleDouble value = DoubleDouble.FromComponents(high, low); + foreach (double scalar in highs.Concat([-high, Math.BitIncrement(high), Math.BitDecrement(high)])) + { + if ((value.High == 0.0 && scalar == 0.0) + || Math.ILogB(value.High) > 1020 || Math.ILogB(scalar) > 1020) + { + continue; // Public signed-zero and boundary dispatch remain. + } + foreach (double signedLow in value.Low == 0.0 ? new[] { 0.0, -0.0 } : new[] { value.Low }) + { + DoubleDouble expected = PreviousFiniteScalarAdd(value.High, signedLow, scalar); + AssertComponentBits(PreciseMathHelper.AddScalar(value.High, signedLow, scalar), expected); + AssertComponentBits(PreciseMathHelper.AddScalarFinite(value.High, signedLow, scalar), expected); + } + } + } + } + } + + private static DoubleDouble PreviousFiniteScalarAdd(double high, double low, double value) + { + // Freeze the pre-extraction scalar sequence and public normalization. + // Do not redirect this reference to the extracted production helper. + (double sum, double error) = PreciseMathHelper.TwoAdd(high, value); + (double middle, double middleError) = PreciseMathHelper.TwoAdd(error, low); + (double result, double residual) = PreciseMathHelper.TwoAdd(sum, middle); + return DoubleDouble.FromComponents(result, residual + middleError); + } + + private static void AssertComponentBits(DoubleDouble actual, DoubleDouble expected) + { + BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(BitConverter.DoubleToInt64Bits(expected.High)); + BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(BitConverter.DoubleToInt64Bits(expected.Low)); + } + private static void AssertDivisionResidualCancellation(DoubleDouble remainder, DoubleDouble product) { DoubleDouble expected = remainder - product; diff --git a/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs b/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs index 327bc39..6282fa3 100644 --- a/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs +++ b/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs @@ -261,6 +261,91 @@ public class PreciseMathSqrtTests } } + [Fact] + public void FiniteNormalizationMatchesBothPreviousSqrtSitesAcrossTheRange() + { + // Replay the unchanged scaled Newton sequence with public factories. + // This pins normalization compatibility, not root accuracy: the exact + // squared-inequality tests above remain the independent accuracy oracle. + bool sawSeparateLow = false; + bool sawZeroCorrection = false; + bool sawNonzeroCorrection = false; + foreach (DoubleDouble input in NormalizationInputs()) + { + int exponent = Math.ILogB(input.High) & ~1; + double high = Math.ScaleB(input.High, -exponent); + double low = Math.ScaleB(input.Low, -exponent); + bool separateLow = exponent > 0 && input.Low != 0.0 && Math.Abs(low) < Math.ScaleB(1.0, -1020); + sawSeparateLow |= separateLow; + if (separateLow) + { + low = 0.0; + } + double estimate = Math.Sqrt(high); + double residual = Math.FusedMultiplyAdd(-estimate, estimate, high) + low; + double correction = residual / (2.0 * estimate); + sawZeroCorrection |= correction == 0.0; + sawNonzeroCorrection |= correction != 0.0; + (estimate >= 1.0 && estimate <= 2.0).ShouldBeTrue(); + (Math.Abs(correction) < Math.ScaleB(1.0, -48)).ShouldBeTrue(); + (estimate + correction < 4.0).ShouldBeTrue(); + DoubleDouble root = AssertFiniteNormalization(estimate, correction); + + double scaledHigh = Math.ScaleB(root.High, exponent / 2); + double scaledLow = Math.ScaleB(root.Low, exponent / 2); + (Math.Abs(scaledHigh) < Math.ScaleB(1.0, 514)).ShouldBeTrue(); + (Math.Abs(scaledLow) < Math.ScaleB(1.0, 461)).ShouldBeTrue(); + (scaledHigh + scaledLow < Math.ScaleB(1.0, 514)).ShouldBeTrue(); + DoubleDouble expected = AssertFiniteNormalization(scaledHigh, scaledLow); + if (!separateLow) + { + DoubleDouble actual = DoubleDouble.Sqrt(input); + BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(BitConverter.DoubleToInt64Bits(expected.High)); + BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(BitConverter.DoubleToInt64Bits(expected.Low)); + } + // The sparse branch also passed both site checks. Its additional + // correction is pinned independently by AssertSparseRoot above. + } + sawSeparateLow.ShouldBeTrue(); + sawZeroCorrection.ShouldBeTrue(); + sawNonzeroCorrection.ShouldBeTrue(); + } + + private static IEnumerable NormalizationInputs() + { + for (int exponent = -1074; exponent <= 1023; ++exponent) + { + foreach (double significand in new[] { 1.0, Math.BitIncrement(1.0), 1.5, Math.BitDecrement(2.0) }) + { + double high = Math.ScaleB(significand, exponent); + double dense = Math.ScaleB(1.0, exponent - 54); + double threshold = Math.ScaleB(1.0, (exponent & ~1) - 1020); + foreach (double low in new[] { 0.0, dense, -dense, double.Epsilon, -double.Epsilon, + Math.BitDecrement(threshold), threshold, Math.BitIncrement(threshold), -threshold }) + { + DoubleDouble input = DoubleDouble.FromComponents(high, low); + if (input.High > 0.0 && DoubleDouble.IsFinite(input)) + { + yield return input; + } + } + } + } + yield return DoubleDouble.FromComponents(double.MaxValue, Math.BitDecrement(Math.ScaleB(1.0, 970))); + } + + private static DoubleDouble AssertFiniteNormalization(double high, double low) + { + double.IsFinite(high).ShouldBeTrue(); + double.IsFinite(low).ShouldBeTrue(); + double.IsFinite(high + low).ShouldBeTrue(); + DoubleDouble expected = DoubleDouble.FromComponents(high, low); + DoubleDouble actual = PreciseMathHelper.NormalizeFinite(high, low); + BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(BitConverter.DoubleToInt64Bits(expected.High)); + BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(BitConverter.DoubleToInt64Bits(expected.Low)); + return expected; + } + private static void AssertSparseRoot(DoubleDouble input, double high, double low) { BigInteger x = (Units(input.High) + Units(input.Low)) << 1074; diff --git a/2-benchmarks/Directory.Build.props b/2-benchmarks/Directory.Build.props index ef46d91..0fb9729 100644 --- a/2-benchmarks/Directory.Build.props +++ b/2-benchmarks/Directory.Build.props @@ -5,4 +5,9 @@ false false + + + + + diff --git a/2-benchmarks/Just.PreciseMath.Benchmarks/Just.PreciseMath.Benchmarks.csproj b/2-benchmarks/Just.PreciseMath.Benchmarks/Just.PreciseMath.Benchmarks.csproj index 3503c0a..0cd995a 100644 --- a/2-benchmarks/Just.PreciseMath.Benchmarks/Just.PreciseMath.Benchmarks.csproj +++ b/2-benchmarks/Just.PreciseMath.Benchmarks/Just.PreciseMath.Benchmarks.csproj @@ -5,9 +5,4 @@ BenchmarkDotNet comparisons of DoubleDouble, decimal, and double arithmetic. Exe - - - - - diff --git a/2-benchmarks/Just.PreciseMath.PerformanceComparison/ArithmeticComparisonBenchmarks.cs b/2-benchmarks/Just.PreciseMath.PerformanceComparison/ArithmeticComparisonBenchmarks.cs new file mode 100644 index 0000000..dbad57f --- /dev/null +++ b/2-benchmarks/Just.PreciseMath.PerformanceComparison/ArithmeticComparisonBenchmarks.cs @@ -0,0 +1,164 @@ +using BenchmarkDotNet.Attributes; +using BenchmarkDotNet.Configs; +using Just.PreciseMath.PerformanceComparison.Reference; + +namespace Just.PreciseMath.PerformanceComparison; + +/// +/// Compares the four arithmetic operators on the same ordinary finite operands across +/// binary64, the Just.PreciseMath double-double, and the blog reference +/// implementation. Every category measures one operator, and its binary64 method is +/// the baseline, so the ratio columns read as "times the cost of a double". +/// +/// +/// Operands are independent array elements, so the measurements include loads, stores, +/// and loop overhead on purpose; input construction, conversion, and fixture validation +/// stay outside the timed methods. +/// +[MemoryDiagnoser] +[GroupBenchmarksBy(BenchmarkLogicalGroupRule.ByCategory)] +[CategoriesColumn] +public class ArithmeticComparisonBenchmarks +{ + private double[] _doubleLeft = []; + private double[] _doubleRight = []; + private DoubleDouble[] _left = []; + private DoubleDouble[] _right = []; + private DoubleDoubleSample[] _sampleLeft = []; + private DoubleDoubleSample[] _sampleRight = []; + private readonly double[] _doubleResults = new double[PerformanceInputs.Count]; + private readonly DoubleDouble[] _results = new DoubleDouble[PerformanceInputs.Count]; + private readonly DoubleDoubleSample[] _sampleResults = new DoubleDoubleSample[PerformanceInputs.Count]; + + [GlobalSetup] + public void Setup() + { + PerformanceInputs inputs = PerformanceInputs.Create(); + _doubleLeft = inputs.DoubleLeft; + _doubleRight = inputs.DoubleRight; + _left = inputs.Left; + _right = inputs.Right; + _sampleLeft = inputs.SampleLeft; + _sampleRight = inputs.SampleRight; + } + + [Benchmark(Baseline = true, OperationsPerInvoke = PerformanceInputs.Count), BenchmarkCategory("Addition")] + public double[] DoubleAdd() + { + for (int i = 0; i < PerformanceInputs.Count; i++) + { + _doubleResults[i] = _doubleLeft[i] + _doubleRight[i]; + } + return _doubleResults; + } + + [Benchmark(OperationsPerInvoke = PerformanceInputs.Count), BenchmarkCategory("Addition")] + public DoubleDouble[] DoubleDoubleAdd() + { + for (int i = 0; i < PerformanceInputs.Count; i++) + { + _results[i] = _left[i] + _right[i]; + } + return _results; + } + + [Benchmark(OperationsPerInvoke = PerformanceInputs.Count), BenchmarkCategory("Addition")] + public DoubleDoubleSample[] DoubleDoubleSampleAdd() + { + for (int i = 0; i < PerformanceInputs.Count; i++) + { + _sampleResults[i] = _sampleLeft[i] + _sampleRight[i]; + } + return _sampleResults; + } + + [Benchmark(Baseline = true, OperationsPerInvoke = PerformanceInputs.Count), BenchmarkCategory("Subtraction")] + public double[] DoubleSubtract() + { + for (int i = 0; i < PerformanceInputs.Count; i++) + { + _doubleResults[i] = _doubleLeft[i] - _doubleRight[i]; + } + return _doubleResults; + } + + [Benchmark(OperationsPerInvoke = PerformanceInputs.Count), BenchmarkCategory("Subtraction")] + public DoubleDouble[] DoubleDoubleSubtract() + { + for (int i = 0; i < PerformanceInputs.Count; i++) + { + _results[i] = _left[i] - _right[i]; + } + return _results; + } + + [Benchmark(OperationsPerInvoke = PerformanceInputs.Count), BenchmarkCategory("Subtraction")] + public DoubleDoubleSample[] DoubleDoubleSampleSubtract() + { + for (int i = 0; i < PerformanceInputs.Count; i++) + { + _sampleResults[i] = _sampleLeft[i] - _sampleRight[i]; + } + return _sampleResults; + } + + [Benchmark(Baseline = true, OperationsPerInvoke = PerformanceInputs.Count), BenchmarkCategory("Multiplication")] + public double[] DoubleMultiply() + { + for (int i = 0; i < PerformanceInputs.Count; i++) + { + _doubleResults[i] = _doubleLeft[i] * _doubleRight[i]; + } + return _doubleResults; + } + + [Benchmark(OperationsPerInvoke = PerformanceInputs.Count), BenchmarkCategory("Multiplication")] + public DoubleDouble[] DoubleDoubleMultiply() + { + for (int i = 0; i < PerformanceInputs.Count; i++) + { + _results[i] = _left[i] * _right[i]; + } + return _results; + } + + [Benchmark(OperationsPerInvoke = PerformanceInputs.Count), BenchmarkCategory("Multiplication")] + public DoubleDoubleSample[] DoubleDoubleSampleMultiply() + { + for (int i = 0; i < PerformanceInputs.Count; i++) + { + _sampleResults[i] = _sampleLeft[i] * _sampleRight[i]; + } + return _sampleResults; + } + + [Benchmark(Baseline = true, OperationsPerInvoke = PerformanceInputs.Count), BenchmarkCategory("Division")] + public double[] DoubleDivide() + { + for (int i = 0; i < PerformanceInputs.Count; i++) + { + _doubleResults[i] = _doubleLeft[i] / _doubleRight[i]; + } + return _doubleResults; + } + + [Benchmark(OperationsPerInvoke = PerformanceInputs.Count), BenchmarkCategory("Division")] + public DoubleDouble[] DoubleDoubleDivide() + { + for (int i = 0; i < PerformanceInputs.Count; i++) + { + _results[i] = _left[i] / _right[i]; + } + return _results; + } + + [Benchmark(OperationsPerInvoke = PerformanceInputs.Count), BenchmarkCategory("Division")] + public DoubleDoubleSample[] DoubleDoubleSampleDivide() + { + for (int i = 0; i < PerformanceInputs.Count; i++) + { + _sampleResults[i] = _sampleLeft[i] / _sampleRight[i]; + } + return _sampleResults; + } +} diff --git a/2-benchmarks/Just.PreciseMath.PerformanceComparison/Just.PreciseMath.PerformanceComparison.csproj b/2-benchmarks/Just.PreciseMath.PerformanceComparison/Just.PreciseMath.PerformanceComparison.csproj new file mode 100644 index 0000000..ed9781c --- /dev/null +++ b/2-benchmarks/Just.PreciseMath.PerformanceComparison/Just.PreciseMath.PerformanceComparison.csproj @@ -0,0 +1,10 @@ + + + + Exe + net10.0 + enable + enable + + + diff --git a/2-benchmarks/Just.PreciseMath.PerformanceComparison/MathComparisonBenchmarks.cs b/2-benchmarks/Just.PreciseMath.PerformanceComparison/MathComparisonBenchmarks.cs new file mode 100644 index 0000000..abbb08c --- /dev/null +++ b/2-benchmarks/Just.PreciseMath.PerformanceComparison/MathComparisonBenchmarks.cs @@ -0,0 +1,135 @@ +using BenchmarkDotNet.Attributes; +using BenchmarkDotNet.Configs; +using Just.PreciseMath.PerformanceComparison.Reference; + +namespace Just.PreciseMath.PerformanceComparison; + +/// +/// Compares Square, Sqrt, and Abs on the same ordinary finite operands across binary64, +/// the Just.PreciseMath double-double, and the blog reference implementation. +/// Every category measures one function, and its binary64 method is the baseline, so +/// the ratio columns read as "times the cost of a double". +/// +/// +/// Operands are independent array elements, so the measurements include loads, stores, +/// and loop overhead on purpose; input construction and fixture validation stay outside +/// the timed methods. +/// +[MemoryDiagnoser] +[GroupBenchmarksBy(BenchmarkLogicalGroupRule.ByCategory)] +[CategoriesColumn] +public class MathComparisonBenchmarks +{ + private double[] _doublePositive = []; + private DoubleDouble[] _positive = []; + private DoubleDoubleSample[] _samplePositive = []; + private double[] _doubleMixed = []; + private DoubleDouble[] _mixed = []; + private DoubleDoubleSample[] _sampleMixed = []; + private readonly double[] _doubleResults = new double[PerformanceInputs.Count]; + private readonly DoubleDouble[] _results = new DoubleDouble[PerformanceInputs.Count]; + private readonly DoubleDoubleSample[] _sampleResults = new DoubleDoubleSample[PerformanceInputs.Count]; + + [GlobalSetup] + public void Setup() + { + PerformanceInputs inputs = PerformanceInputs.Create(); + _doublePositive = inputs.DoublePositive; + _positive = inputs.Positive; + _samplePositive = inputs.SamplePositive; + // Square and Sqrt need positive operands; Abs also covers negative ones. + _doubleMixed = inputs.DoubleLeft; + _mixed = inputs.Left; + _sampleMixed = inputs.SampleLeft; + } + + [Benchmark(Baseline = true, OperationsPerInvoke = PerformanceInputs.Count), BenchmarkCategory("Square")] + public double[] DoubleSquare() + { + for (int i = 0; i < PerformanceInputs.Count; i++) + { + _doubleResults[i] = _doublePositive[i] * _doublePositive[i]; + } + return _doubleResults; + } + + [Benchmark(OperationsPerInvoke = PerformanceInputs.Count), BenchmarkCategory("Square")] + public DoubleDouble[] DoubleDoubleSquare() + { + for (int i = 0; i < PerformanceInputs.Count; i++) + { + _results[i] = DoubleDouble.Square(_positive[i]); + } + return _results; + } + + [Benchmark(OperationsPerInvoke = PerformanceInputs.Count), BenchmarkCategory("Square")] + public DoubleDoubleSample[] DoubleDoubleSampleSquare() + { + for (int i = 0; i < PerformanceInputs.Count; i++) + { + _sampleResults[i] = _samplePositive[i].Square(); + } + return _sampleResults; + } + + [Benchmark(Baseline = true, OperationsPerInvoke = PerformanceInputs.Count), BenchmarkCategory("Sqrt")] + public double[] DoubleSqrt() + { + for (int i = 0; i < PerformanceInputs.Count; i++) + { + _doubleResults[i] = Math.Sqrt(_doublePositive[i]); + } + return _doubleResults; + } + + [Benchmark(OperationsPerInvoke = PerformanceInputs.Count), BenchmarkCategory("Sqrt")] + public DoubleDouble[] DoubleDoubleSqrt() + { + for (int i = 0; i < PerformanceInputs.Count; i++) + { + _results[i] = DoubleDouble.Sqrt(_positive[i]); + } + return _results; + } + + [Benchmark(OperationsPerInvoke = PerformanceInputs.Count), BenchmarkCategory("Sqrt")] + public DoubleDoubleSample[] DoubleDoubleSampleSqrt() + { + for (int i = 0; i < PerformanceInputs.Count; i++) + { + _sampleResults[i] = _samplePositive[i].Sqrt(); + } + return _sampleResults; + } + + [Benchmark(Baseline = true, OperationsPerInvoke = PerformanceInputs.Count), BenchmarkCategory("Abs")] + public double[] DoubleAbs() + { + for (int i = 0; i < PerformanceInputs.Count; i++) + { + _doubleResults[i] = Math.Abs(_doubleMixed[i]); + } + return _doubleResults; + } + + [Benchmark(OperationsPerInvoke = PerformanceInputs.Count), BenchmarkCategory("Abs")] + public DoubleDouble[] DoubleDoubleAbs() + { + for (int i = 0; i < PerformanceInputs.Count; i++) + { + _results[i] = DoubleDouble.Abs(_mixed[i]); + } + return _results; + } + + [Benchmark(OperationsPerInvoke = PerformanceInputs.Count), BenchmarkCategory("Abs")] + public DoubleDoubleSample[] DoubleDoubleSampleAbs() + { + for (int i = 0; i < PerformanceInputs.Count; i++) + { + _sampleResults[i] = _sampleMixed[i].Abs(); + } + return _sampleResults; + } +} diff --git a/2-benchmarks/Just.PreciseMath.PerformanceComparison/PerformanceInputs.cs b/2-benchmarks/Just.PreciseMath.PerformanceComparison/PerformanceInputs.cs new file mode 100644 index 0000000..b7838c5 --- /dev/null +++ b/2-benchmarks/Just.PreciseMath.PerformanceComparison/PerformanceInputs.cs @@ -0,0 +1,167 @@ +using Just.PreciseMath.PerformanceComparison.Reference; + +namespace Just.PreciseMath.PerformanceComparison; + +/// +/// Deterministic operand sets for the implementation comparison: +/// values per operation, each with a nonzero high and low component, built from exact +/// dyadic multiples and held inside a moderate exponent range so that no benchmark +/// path takes an overflow, underflow, or boundary branch. +/// +/// +/// Both double-double implementations receive bit-identical operands, and the plain +/// doubles equal the high components, so the binary64 reference measures the leading +/// component of the same values. The mixed-sign set feeds +, -, *, / and Abs; the +/// nonnegative set feeds Square and Sqrt. +/// +internal sealed class PerformanceInputs +{ + /// Operands per benchmark invocation. One double-double operation takes only + /// a few nanoseconds, so an invocation covers many values and the measurement is + /// normalized back to a single operation. + internal const int Count = 200; + + /// Mixed-sign operands for +, -, *, / and Abs. + internal DoubleDouble[] Left { get; } = new DoubleDouble[Count]; + + internal DoubleDouble[] Right { get; } = new DoubleDouble[Count]; + + /// Nonnegative operands for Square and Sqrt, with both components positive. + internal DoubleDouble[] Positive { get; } = new DoubleDouble[Count]; + + /// The same operands for the blog reference implementation. + internal DoubleDoubleSample[] SampleLeft { get; } = new DoubleDoubleSample[Count]; + + internal DoubleDoubleSample[] SampleRight { get; } = new DoubleDoubleSample[Count]; + + internal DoubleDoubleSample[] SamplePositive { get; } = new DoubleDoubleSample[Count]; + + /// The high components as plain doubles, for the binary64 reference. + internal double[] DoubleLeft { get; } = new double[Count]; + + internal double[] DoubleRight { get; } = new double[Count]; + + internal double[] DoublePositive { get; } = new double[Count]; + + internal static PerformanceInputs Create() + { + PerformanceInputs inputs = new PerformanceInputs(); + for (int i = 0; i < Count; i++) + { + // Highs are small dyadic multiples of a power of two: exact in binary64, + // with the significand needing at most nine bits. Lows are multiples of + // 2^(exponent-57) or 2^(exponent-60), at least two bits below half an ulp + // of their high, so every pair is normalized and both components survive + // every operation. The exponents keep highs within 2^+-48, which leaves + // products, quotients, and squares far from the binary64 exponent limits. + int leftExponent = (((i / 4) % 5) - 2) * 24; + int rightExponent = (((i / 4) % 3) - 1) * 16; + int positiveExponent = ((i % 5) - 2) * 24; + double leftSign = (i % 2 == 0) ? 1.0 : -1.0; + double rightSign = (i % 3 == 0) ? 1.0 : -1.0; + double leftHigh = leftSign * Math.ScaleB(1.25 + ((i % 7) / 32.0), leftExponent); + double leftLow = leftSign * Math.ScaleB(1.0 + (i % 4), leftExponent - 57); + double rightHigh = rightSign * Math.ScaleB(1.0 + (((i % 5) - 2) / 128.0), rightExponent); + // Opposite sign, and a three-bit finer scale than the left low: no + // addition or multiplication of the two lows can cancel to a value that + // binary64 would represent without a residual. + double rightLow = -rightSign * Math.ScaleB(1.0 + ((i / 2) % 3), rightExponent - 60); + double positiveHigh = Math.ScaleB(1.25 + ((i % 7) / 32.0), positiveExponent); + double positiveLow = Math.ScaleB(1.0 + (i % 4), positiveExponent - 57); + + inputs.Left[i] = DoubleDouble.FromComponents(leftHigh, leftLow); + inputs.Right[i] = DoubleDouble.FromComponents(rightHigh, rightLow); + inputs.Positive[i] = DoubleDouble.FromComponents(positiveHigh, positiveLow); + inputs.SampleLeft[i] = DoubleDoubleSample.FromTwoDoubles(leftHigh, leftLow); + inputs.SampleRight[i] = DoubleDoubleSample.FromTwoDoubles(rightHigh, rightLow); + inputs.SamplePositive[i] = DoubleDoubleSample.FromTwoDoubles(positiveHigh, positiveLow); + inputs.DoubleLeft[i] = leftHigh; + inputs.DoubleRight[i] = rightHigh; + inputs.DoublePositive[i] = positiveHigh; + } + + inputs.Validate(); + return inputs; + } + + /// + /// Setup-only fixture guard. It checks the operand shape and that every benchmarked + /// operation consumes both components; it is not an accuracy oracle and does not + /// compare the two implementations against each other. + /// + private void Validate() + { + for (int i = 0; i < Count; i++) + { + ValidateOperand(Left[i], nameof(Left), i); + ValidateOperand(Right[i], nameof(Right), i); + ValidateOperand(Positive[i], nameof(Positive), i); + ValidateMatch(Left[i], SampleLeft[i], nameof(Left), i); + ValidateMatch(Right[i], SampleRight[i], nameof(Right), i); + ValidateMatch(Positive[i], SamplePositive[i], nameof(Positive), i); + if (DoubleLeft[i] != Left[i].High || DoubleRight[i] != Right[i].High || DoublePositive[i] != Positive[i].High) + { + throw new InvalidOperationException($"[{i}] the binary64 operands do not match the leading components."); + } + if (!(Positive[i].High > 0.0) || !(Positive[i].Low > 0.0)) + { + throw new InvalidOperationException($"[{i}] Square and Sqrt operands must be positive in both components."); + } + } + + for (int i = 0; i < Count; i++) + { + CheckResidual("Left + Right", Left[i] + Right[i], i); + CheckResidual("Left - Right", Left[i] - Right[i], i); + CheckResidual("Left * Right", Left[i] * Right[i], i); + CheckResidual("Left / Right", Left[i] / Right[i], i); + CheckResidual("Positive * Positive", DoubleDouble.Square(Positive[i]), i); + CheckResidual("Sqrt(Positive)", DoubleDouble.Sqrt(Positive[i]), i); + CheckResidual("Abs(Left)", DoubleDouble.Abs(Left[i]), i); + CheckResidual("SampleLeft + SampleRight", SampleLeft[i] + SampleRight[i], i); + CheckResidual("SampleLeft - SampleRight", SampleLeft[i] - SampleRight[i], i); + CheckResidual("SampleLeft * SampleRight", SampleLeft[i] * SampleRight[i], i); + CheckResidual("SampleLeft / SampleRight", SampleLeft[i] / SampleRight[i], i); + CheckResidual("SamplePositive.Square()", SamplePositive[i].Square(), i); + CheckResidual("SamplePositive.Sqrt()", SamplePositive[i].Sqrt(), i); + CheckResidual("SampleLeft.Abs()", SampleLeft[i].Abs(), i); + } + } + + private static void ValidateOperand(DoubleDouble value, string name, int index) + { + // FromComponents renormalizes an arbitrary pair, so reproducing the value + // bit-for-bit is the normalization check for a purely dyadic construction. + bool normalized = DoubleDouble.FromComponents(value.High, value.Low).Equals(value); + if (!DoubleDouble.IsFinite(value) || value.High == 0.0 || value.Low == 0.0 || !normalized + || Math.Abs(Math.ILogB(value.High)) > 100) + { + throw new InvalidOperationException($"{name}[{index}] is not an ordinary normalized value with both components used."); + } + } + + private static void ValidateMatch(DoubleDouble value, DoubleDoubleSample sample, string name, int index) + { + bool normalized = DoubleDoubleSample.FromTwoDoubles(sample.Hi, sample.Lo).Equals(sample); + if (sample.Hi != value.High || sample.Lo != value.Low || !normalized) + { + throw new InvalidOperationException($"{name}[{index}] differs between the two double-double implementations."); + } + } + + private static void CheckResidual(string operation, DoubleDouble result, int index) + { + if (!DoubleDouble.IsFinite(result) || result.Low == 0.0) + { + throw new InvalidOperationException($"{operation}[{index}] left the ordinary range or lost its low component."); + } + } + + private static void CheckResidual(string operation, DoubleDoubleSample result, int index) + { + if (!double.IsFinite(result.Hi) || result.Lo == 0.0) + { + throw new InvalidOperationException($"{operation}[{index}] left the ordinary range or lost its low component."); + } + } +} diff --git a/2-benchmarks/Just.PreciseMath.PerformanceComparison/Program.cs b/2-benchmarks/Just.PreciseMath.PerformanceComparison/Program.cs new file mode 100644 index 0000000..c981e50 --- /dev/null +++ b/2-benchmarks/Just.PreciseMath.PerformanceComparison/Program.cs @@ -0,0 +1,23 @@ +using BenchmarkDotNet.Reports; +using BenchmarkDotNet.Running; + +BenchmarkSwitcher switcher = BenchmarkSwitcher.FromAssembly(typeof(Program).Assembly); +// Only these standalone discovery/help forms intentionally produce no reports. +if (args is ["--list", "flat" or "tree"] or ["--help"] or ["-h"]) +{ + switcher.Run(args); + return 0; +} + +Summary[] summaries = [.. switcher.Run(args.Length == 0 ? ["--filter", "*"] : args)]; +// BenchmarkDotNet can return normally after validation, build, or execution +// failures. Do not let CI mistake an empty or failed run for a passing smoke test. +if (summaries.Length == 0 || summaries.Any(summary => + summary.HasCriticalValidationErrors || summary.Reports.Length == 0 || + summary.Reports.Any(report => !report.Success || report.ResultStatistics is null))) +{ + Console.Error.WriteLine("Benchmark run failed or produced no measurements."); + return 1; +} + +return 0; diff --git a/2-benchmarks/Just.PreciseMath.PerformanceComparison/Reference/DoubleDoubleSample.Format.cs b/2-benchmarks/Just.PreciseMath.PerformanceComparison/Reference/DoubleDoubleSample.Format.cs new file mode 100644 index 0000000..9fe8145 --- /dev/null +++ b/2-benchmarks/Just.PreciseMath.PerformanceComparison/Reference/DoubleDoubleSample.Format.cs @@ -0,0 +1,278 @@ +// Printing and parsing of DoubleDoubleSample, by Marek Fiser. +// Source: https://marekfiser.com/blog/double-double-arithmetic +// Licensed under CC BY 4.0: https://creativecommons.org/licenses/by/4.0/ +// Changes for this benchmark: copied into this project, kept in the +// Just.PreciseMath.PerformanceComparison.Reference namespace, and the file and type +// renamed to DoubleDoubleSample.Format.cs and DoubleDoubleSample. +// Printing 31 correct decimal digits requires computing them exactly, which is done with +// BigInteger (part of the BCL, so still no NuGet dependencies). +using System.Globalization; +using System.Numerics; +using System.Text; + +namespace Just.PreciseMath.PerformanceComparison.Reference; + +public readonly partial struct DoubleDoubleSample { + + /// Significant digits prints: enough to round-trip any + /// value whose two components are contiguous (see ). + private const int DEFAULT_DIGITS = 36; + + + /// Formats with significant digits. + public override string ToString() { + return ToString(DEFAULT_DIGITS); + } + + /// Formats with the given number of significant decimal digits, correctly rounded. + public string ToString(int significantDigits) { + ArgumentOutOfRangeException.ThrowIfLessThan(significantDigits, 1); + + if (!double.IsFinite(m_hi)) { + return m_hi.ToString(CultureInfo.InvariantCulture); + } + + (BigInteger numerator, BigInteger denominator) = toRational(this); + return formatRational(numerator, denominator, significantDigits); + } + + /// + /// Formats the exact value. A double-double is a sum of two binary fractions, so it + /// always has a finite decimal expansion, but it can be hundreds of digits long when the + /// components are far apart. + /// + public string ToStringExact() { + if (!double.IsFinite(m_hi)) { + return m_hi.ToString(CultureInfo.InvariantCulture); + } + + (BigInteger numerator, BigInteger denominator) = toRational(this); + if (numerator.IsZero) { + return "0"; + } + + // The denominator is a power of two, so multiplying by the matching power of five + // turns it into a power of ten, which gives the decimal digits directly. + int twos = (int)(denominator.GetBitLength() - 1); + BigInteger scaled = BigInteger.Abs(numerator) * BigInteger.Pow(5, twos); + string digits = scaled.ToString(CultureInfo.InvariantCulture).PadLeft(twos + 1, '0'); + string sign = numerator.Sign < 0 ? "-" : ""; + if (twos == 0) { + return sign + digits; + } + + string result = sign + digits[..^twos] + "." + digits[^twos..]; + return result.TrimEnd('0').TrimEnd('.'); + } + + /// + /// Parses a decimal string to the nearest double-double. The value is built exactly as a + /// rational, then rounded once into each component. Round-trips with + /// for every value whose components are contiguous; + /// a hand-built pair such as (1.0, 1e-300) holds more information than 36 digits can + /// represent. + /// + public static DoubleDoubleSample Parse(string text) { + (BigInteger numerator, BigInteger denominator) = parseRational(text); + double hi = toNearestDouble(numerator, denominator); + if (!double.IsFinite(hi)) { + return FromDouble(hi); + } + + // Subtract the leading component exactly and round the remainder into the trailing one. + (BigInteger hiNumerator, BigInteger hiDenominator) = toRational(FromDouble(hi)); + BigInteger restNumerator = numerator * hiDenominator - hiNumerator * denominator; + BigInteger restDenominator = denominator * hiDenominator; + double lo = toNearestDouble(restNumerator, restDenominator); + return FromTwoDoubles(hi, lo); + } + + /// The exact value of a double-double as a rational (the denominator is a power of two). + private static (BigInteger Numerator, BigInteger Denominator) toRational(DoubleDoubleSample value) { + (BigInteger hiNumerator, BigInteger hiDenominator) = toRational(value.m_hi); + (BigInteger loNumerator, BigInteger loDenominator) = toRational(value.m_lo); + return (hiNumerator * loDenominator + loNumerator * hiDenominator, hiDenominator * loDenominator); + } + + /// The exact value of a finite double as a rational, from its bits. + private static (BigInteger Numerator, BigInteger Denominator) toRational(double value) { + long bits = BitConverter.DoubleToInt64Bits(value); + int exponent = (int)((bits >> 52) & 0x7FF); + long mantissa = bits & 0xF_FFFF_FFFF_FFFFL; + if (exponent == 0) { + exponent = -1074; // subnormal: no implicit leading bit + } else { + mantissa |= 1L << 52; + exponent -= 1075; + } + + BigInteger numerator = bits < 0 ? -mantissa : mantissa; + return exponent >= 0 + ? (numerator << exponent, BigInteger.One) + : (numerator, BigInteger.One << -exponent); + } + + /// Rounds an exact rational to the nearest double, ties to even. + private static double toNearestDouble(BigInteger numerator, BigInteger denominator) { + if (numerator.IsZero) { + return 0.0; + } + + int sign = numerator.Sign; + numerator = BigInteger.Abs(numerator); + + // Scale so that the quotient is in [2^52, 2^54). One integer division then gives all + // bits of the significand plus the remainder needed for the rounding. + long exponent = numerator.GetBitLength() - denominator.GetBitLength() - 53; + if (exponent > 0) { + denominator <<= (int)exponent; + } else { + numerator <<= (int)-exponent; + } + + BigInteger quotient = BigInteger.DivRem(numerator, denominator, out BigInteger remainder); + if (quotient.GetBitLength() > 53) { + // One bit too many. Move it to the exponent and add the dropped bit to the + // remainder so that the rounding decision below stays exact. + if (!quotient.IsEven) { + remainder += denominator; + } + quotient >>= 1; + denominator <<= 1; + exponent++; + } + + int comparison = (remainder << 1).CompareTo(denominator); + if (comparison > 0 || (comparison == 0 && !quotient.IsEven)) { + quotient++; + if (quotient.GetBitLength() > 53) { + quotient >>= 1; + exponent++; + } + } + + return sign * Math.ScaleB((double)quotient, (int)exponent); + } + + /// Parses a decimal string such as "-1.25e-7" into an exact rational. + private static (BigInteger Numerator, BigInteger Denominator) parseRational(string text) { + text = text.Trim(); + if (text.Length == 0) { + throw new FormatException("Empty input."); + } + + int index = 0; + bool negative = text[index] == '-'; + if (negative || text[index] == '+') { + index++; + } + + BigInteger mantissa = BigInteger.Zero; + int fractionDigits = 0; + bool seenDot = false; + bool seenDigit = false; + + for (; index < text.Length; index++) { + char c = text[index]; + if (c == '.') { + if (seenDot) { + throw new FormatException($"Two decimal points in '{text}'."); + } + seenDot = true; + } else if (c is >= '0' and <= '9') { + mantissa = mantissa * 10 + (c - '0'); + seenDigit = true; + if (seenDot) { + fractionDigits++; + } + } else if (c is 'e' or 'E') { + break; + } else { + throw new FormatException($"Unexpected character '{c}' in '{text}'."); + } + } + + if (!seenDigit) { + throw new FormatException($"No digits in '{text}'."); + } + + int exponent = -fractionDigits; + if (index < text.Length) { + exponent += int.Parse(text[(index + 1)..], CultureInfo.InvariantCulture); + } + + if (negative) { + mantissa = -mantissa; + } + + return exponent >= 0 + ? (mantissa * BigInteger.Pow(10, exponent), BigInteger.One) + : (mantissa, BigInteger.Pow(10, -exponent)); + } + + /// Rounds an exact rational to N significant digits and formats them. + private static string formatRational(BigInteger numerator, BigInteger denominator, int digits) { + if (numerator.IsZero) { + return "0"; + } + + bool negative = numerator.Sign < 0; + numerator = BigInteger.Abs(numerator); + + // Estimate of the decimal exponent of the leading digit from the bit lengths + // (log10(2) ~ 0.30103). The loop below corrects it when the estimate is off by one. + int exponent10 = (int)Math.Floor((numerator.GetBitLength() - denominator.GetBitLength()) * 0.30103); + string text; + + while (true) { + int scale = digits - 1 - exponent10; + BigInteger scaledNumerator = numerator; + BigInteger scaledDenominator = denominator; + if (scale >= 0) { + scaledNumerator *= BigInteger.Pow(10, scale); + } else { + scaledDenominator *= BigInteger.Pow(10, -scale); + } + BigInteger rounded = roundToNearest(scaledNumerator, scaledDenominator); + text = rounded.ToString(CultureInfo.InvariantCulture); + if (text.Length == digits) { + break; + } + // Off by one either way (estimate too low, or rounding carried into a new digit). + exponent10 += text.Length - digits; + } + + text = text.TrimEnd('0'); + if (text.Length == 0) { + text = "0"; + } + + StringBuilder result = new(negative ? "-" : ""); + if (exponent10 >= -5 && exponent10 < digits) { + // Positional notation for the range where it is readable. + if (exponent10 >= 0) { + text = text.PadRight(exponent10 + 1, '0'); + result.Append(text[..(exponent10 + 1)]); + if (text.Length > exponent10 + 1) { + result.Append('.').Append(text[(exponent10 + 1)..]); + } + } else { + result.Append("0.").Append('0', -exponent10 - 1).Append(text); + } + } else { + result.Append(text[0]); + if (text.Length > 1) { + result.Append('.').Append(text[1..]); + } + result.Append('e').Append(exponent10.ToString(CultureInfo.InvariantCulture)); + } + + return result.ToString(); + } + + /// Integer nearest to numerator/denominator, ties away from zero (both positive). + private static BigInteger roundToNearest(BigInteger numerator, BigInteger denominator) { + BigInteger quotient = BigInteger.DivRem(numerator, denominator, out BigInteger remainder); + return (remainder << 1) >= denominator ? quotient + 1 : quotient; + } +} diff --git a/2-benchmarks/Just.PreciseMath.PerformanceComparison/Reference/DoubleDoubleSample.cs b/2-benchmarks/Just.PreciseMath.PerformanceComparison/Reference/DoubleDoubleSample.cs new file mode 100644 index 0000000..5292c60 --- /dev/null +++ b/2-benchmarks/Just.PreciseMath.PerformanceComparison/Reference/DoubleDoubleSample.cs @@ -0,0 +1,213 @@ +// Double-double arithmetic in about 200 lines, by Marek Fiser. +// Source: https://marekfiser.com/blog/double-double-arithmetic +// Licensed under CC BY 4.0: https://creativecommons.org/licenses/by/4.0/ +// Changes for this benchmark: copied into this project, kept in the +// Just.PreciseMath.PerformanceComparison.Reference namespace, and the file and type +// renamed to DoubleDoubleSample. +// This file depends only on System so that it can be copied into any project. Printing and +// parsing use BigInteger and are in DoubleDoubleSample.Format.cs. +namespace Just.PreciseMath.PerformanceComparison.Reference; + +/// +/// A value stored as the unevaluated sum of two doubles: ~106 significand bits +/// (~31 decimal digits) at a few times the cost of a double, with no heap allocation. +/// +/// Invariant: Hi == RoundToNearest(Hi + Lo), i.e. |Lo| is at most half an ulp of Hi. +/// Every factory and operation below maintains it; comparisons and ToDouble rely on it. +/// +/// Operations are not correctly rounded but stay within a few units of 2^-106 of the +/// exact result. The exponent range is the same as a plain double's, and infinities are +/// not special-cased: the first operation on one returns NaN (the error terms compute +/// inf - inf). +/// +/// By Marek Fiser, from marekfiser.com/blog/double-double-arithmetic/ (CC BY 4.0; the file +/// header carries the full attribution and license notice). +/// +public readonly partial struct DoubleDoubleSample : IEquatable, IComparable { + + private readonly double m_hi; + private readonly double m_lo; + + /// The leading component; equals the value rounded to nearest double. + public double Hi => m_hi; + + /// The trailing component (the rounding error of ). + public double Lo => m_lo; + + public bool IsZero => m_hi == 0.0; + + /// -1, 0 or +1 (0 also for NaN). + public int Sign => m_hi > 0.0 ? 1 : m_hi < 0.0 ? -1 : 0; + + public static DoubleDoubleSample Zero { get; } = new DoubleDoubleSample(0.0, 0.0); + + public static DoubleDoubleSample One { get; } = new DoubleDoubleSample(1.0, 0.0); + + + private DoubleDoubleSample(double hi, double lo) { + m_hi = hi; + m_lo = lo; + } + + + /// Exact conversion from a double. + public static DoubleDoubleSample FromDouble(double value) { + return new DoubleDoubleSample(value, 0.0); + } + + /// The exact sum + of two arbitrary doubles. + public static DoubleDoubleSample FromTwoDoubles(double hi, double lo) { + // Renormalizes so that the invariant holds. + (double sum, double error) = twoSum(hi, lo); + return new DoubleDoubleSample(sum, error); + } + + public DoubleDoubleSample Add(DoubleDoubleSample right) { + // Add the two components pairwise, then incorporate lowSum and lowError one by one. + // Adding lowError can round, so a second quickTwoSum renormalizes the result. + (double sum, double error) = twoSum(m_hi, right.m_hi); + (double lowSum, double lowError) = twoSum(m_lo, right.m_lo); + error += lowSum; + (sum, error) = quickTwoSum(sum, error); + error += lowError; + (sum, error) = quickTwoSum(sum, error); + return new DoubleDoubleSample(sum, error); + } + + public DoubleDoubleSample Subtract(DoubleDoubleSample right) { + return Add(right.Negate()); + } + + public DoubleDoubleSample Multiply(DoubleDoubleSample right) { + // hi*hi exactly, then the three cross terms. Each cross term is at least 2^-53 smaller + // than the product, so plain double arithmetic is accurate enough for them. + (double product, double error) = twoProduct(m_hi, right.m_hi); + error += m_hi * right.m_lo + m_lo * right.m_hi + m_lo * right.m_lo; + (product, error) = quickTwoSum(product, error); + return new DoubleDoubleSample(product, error); + } + + /// Multiplies by a plain double. Faster than the full double-double product. + public DoubleDoubleSample Multiply(double right) { + (double product, double error) = twoProduct(m_hi, right); + error += m_lo * right; + (product, error) = quickTwoSum(product, error); + return new DoubleDoubleSample(product, error); + } + + /// Slower and slightly less accurate than the other operations. + public DoubleDoubleSample Divide(DoubleDoubleSample right) { + // Three steps of long division. + double quotient1 = m_hi / right.m_hi; + DoubleDoubleSample remainder = Subtract(right.Multiply(quotient1)); + double quotient2 = remainder.m_hi / right.m_hi; + remainder = remainder.Subtract(right.Multiply(quotient2)); + double quotient3 = remainder.m_hi / right.m_hi; + + (double sum, double error) = quickTwoSum(quotient1, quotient2); + return new DoubleDoubleSample(sum, error).Add(FromDouble(quotient3)); + } + + /// Slightly faster than the general product. + public DoubleDoubleSample Square() { + (double product, double error) = twoProduct(m_hi, m_hi); + error += 2.0 * m_hi * m_lo + m_lo * m_lo; + (product, error) = quickTwoSum(product, error); + return new DoubleDoubleSample(product, error); + } + + public DoubleDoubleSample Negate() { + // Negation never rounds, so this is exact. + return new DoubleDoubleSample(-m_hi, -m_lo); + } + + public DoubleDoubleSample Abs() { + return Sign < 0 ? Negate() : this; + } + + /// Slightly less accurate than the other operations. Negative and NaN inputs + /// return NaN, as Math.Sqrt does. + public DoubleDoubleSample Sqrt() { + if (!(m_hi > 0.0)) { // Not "m_hi <= 0.0", because this form is also true for NaN. + return FromDouble(Math.Sqrt(m_hi)); + } + // One Newton step from the double square root. + double approx = Math.Sqrt(m_hi); + DoubleDoubleSample residual = Subtract(FromDouble(approx).Square()); + double correction = residual.m_hi / (approx + approx); + (double hi, double lo) = quickTwoSum(approx, correction); + return new DoubleDoubleSample(hi, lo); + } + + public double ToDouble() { + return m_hi; + } + + public int CompareTo(DoubleDoubleSample other) { + // Correct only for normalized pairs. + int hiComparison = m_hi.CompareTo(other.m_hi); + return hiComparison != 0 ? hiComparison : m_lo.CompareTo(other.m_lo); + } + + public bool Equals(DoubleDoubleSample other) { + return m_hi == other.m_hi && m_lo == other.m_lo; + } + + public override bool Equals(object? obj) { + return obj is DoubleDoubleSample other && Equals(other); + } + + public override int GetHashCode() { + return HashCode.Combine(m_hi, m_lo); + } + + public static DoubleDoubleSample operator +(DoubleDoubleSample left, DoubleDoubleSample right) => left.Add(right); + public static DoubleDoubleSample operator -(DoubleDoubleSample left, DoubleDoubleSample right) => left.Subtract(right); + public static DoubleDoubleSample operator *(DoubleDoubleSample left, DoubleDoubleSample right) => left.Multiply(right); + public static DoubleDoubleSample operator *(DoubleDoubleSample left, double right) => left.Multiply(right); + public static DoubleDoubleSample operator /(DoubleDoubleSample left, DoubleDoubleSample right) => left.Divide(right); + public static DoubleDoubleSample operator -(DoubleDoubleSample value) => value.Negate(); + + public static bool operator ==(DoubleDoubleSample left, DoubleDoubleSample right) => left.Equals(right); + public static bool operator !=(DoubleDoubleSample left, DoubleDoubleSample right) => !left.Equals(right); + public static bool operator <(DoubleDoubleSample left, DoubleDoubleSample right) => left.CompareTo(right) < 0; + public static bool operator >(DoubleDoubleSample left, DoubleDoubleSample right) => left.CompareTo(right) > 0; + public static bool operator <=(DoubleDoubleSample left, DoubleDoubleSample right) => left.CompareTo(right) <= 0; + public static bool operator >=(DoubleDoubleSample left, DoubleDoubleSample right) => left.CompareTo(right) >= 0; + + // Conversion from a double is exact, so it is implicit. Conversion to a double drops the + // low component, so it is explicit. + public static implicit operator DoubleDoubleSample(double value) => FromDouble(value); + public static explicit operator double(DoubleDoubleSample value) => value.ToDouble(); + + /// Knuth's two-sum: sum + error == a + b exactly, for any two doubles. + private static (double Sum, double Error) twoSum(double a, double b) { + // The rounding error of a floating-point addition is itself representable as a double. + // The intermediate results below may round, but their rounding errors cancel out. + double sum = a + b; + double bKept = sum - a; // The part of b that made it into the sum. + double aKept = sum - bKept; // The part of a that made it into the sum. + double aLost = a - aKept; + double bLost = b - bKept; + return (sum, aLost + bLost); + } + + /// Dekker's fast two-sum: the same guarantee as twoSum, valid only when + /// |a| >= |b| (or a == 0). + private static (double Sum, double Error) quickTwoSum(double a, double b) { + double sum = a + b; + double bKept = sum - a; + return (sum, b - bKept); // a lost nothing, so only b's loss remains. + } + + /// Exact product: product + error == a * b exactly, as long as the product is a + /// normal double. Below about 2e-292 the error no longer fits a double. + private static (double Product, double Error) twoProduct(double a, double b) { + // The fused multiply-add computes a*b - product with a single rounding, which is the + // low half of the 106-bit product. On a CPU without an FMA instruction, + // Math.FusedMultiplyAdd computes the same result in software, much slower. + double product = a * b; + double error = Math.FusedMultiplyAdd(a, b, -product); + return (product, error); + } +} \ No newline at end of file diff --git a/2-benchmarks/Just.PreciseMath.PerformanceComparison/packages.lock.json b/2-benchmarks/Just.PreciseMath.PerformanceComparison/packages.lock.json new file mode 100644 index 0000000..dc77c2b --- /dev/null +++ b/2-benchmarks/Just.PreciseMath.PerformanceComparison/packages.lock.json @@ -0,0 +1,174 @@ +{ + "version": 1, + "dependencies": { + "net10.0": { + "BenchmarkDotNet": { + "type": "Direct", + "requested": "[0.15.8, )", + "resolved": "0.15.8", + "contentHash": "paCfrWxSeHqn3rUZc0spYXVFnHCF0nzRhG0nOLnyTjZYs8spsimBaaNmb3vwqvALKIplbYq/TF393vYiYSnh/Q==", + "dependencies": { + "BenchmarkDotNet.Annotations": "0.15.8", + "CommandLineParser": "2.9.1", + "Gee.External.Capstone": "2.3.0", + "Iced": "1.21.0", + "Microsoft.CodeAnalysis.CSharp": "4.14.0", + "Microsoft.Diagnostics.Runtime": "3.1.512801", + "Microsoft.Diagnostics.Tracing.TraceEvent": "3.1.21", + "Microsoft.DotNet.PlatformAbstractions": "3.1.6", + "Perfolizer": "[0.6.1]", + "System.Management": "9.0.5" + } + }, + "BenchmarkDotNet.Annotations": { + "type": "Transitive", + "resolved": "0.15.8", + "contentHash": "hfucY0ycAsB0SsoaZcaAp9oq5wlWBJcylvEJb9pmvdYUx6PD6S4mDiYnZWjdjAlLhIpe/xtGCwzORfzAzPqvzA==" + }, + "CommandLineParser": { + "type": "Transitive", + "resolved": "2.9.1", + "contentHash": "OE0sl1/sQ37bjVsPKKtwQlWDgqaxWgtme3xZz7JssWUzg5JpMIyHgCTY9MVMxOg48fJ1AgGT3tgdH5m/kQ5xhA==" + }, + "Gee.External.Capstone": { + "type": "Transitive", + "resolved": "2.3.0", + "contentHash": "2ap/rYmjtzCOT8hxrnEW/QeiOt+paD8iRrIcdKX0cxVwWLFa1e+JDBNeECakmccXrSFeBQuu5AV8SNkipFMMMw==" + }, + "Iced": { + "type": "Transitive", + "resolved": "1.21.0", + "contentHash": "dv5+81Q1TBQvVMSOOOmRcjJmvWcX3BZPZsIq31+RLc5cNft0IHAyNlkdb7ZarOWG913PyBoFDsDXoCIlKmLclg==" + }, + "Microsoft.CodeAnalysis.Analyzers": { + "type": "Transitive", + "resolved": "3.11.0", + "contentHash": "v/EW3UE8/lbEYHoC2Qq7AR/DnmvpgdtAMndfQNmpuIMx/Mto8L5JnuCfdBYtgvalQOtfNCnxFejxuRrryvUTsg==" + }, + "Microsoft.CodeAnalysis.Common": { + "type": "Transitive", + "resolved": "4.14.0", + "contentHash": "PC3tuwZYnC+idaPuoC/AZpEdwrtX7qFpmnrfQkgobGIWiYmGi5MCRtl5mx6QrfMGQpK78X2lfIEoZDLg/qnuHg==", + "dependencies": { + "Microsoft.CodeAnalysis.Analyzers": "3.11.0" + } + }, + "Microsoft.CodeAnalysis.CSharp": { + "type": "Transitive", + "resolved": "4.14.0", + "contentHash": "568a6wcTivauIhbeWcCwfWwIn7UV7MeHEBvFB2uzGIpM2OhJ4eM/FZ8KS0yhPoNxnSpjGzz7x7CIjTxhslojQA==", + "dependencies": { + "Microsoft.CodeAnalysis.Analyzers": "3.11.0", + "Microsoft.CodeAnalysis.Common": "[4.14.0]" + } + }, + "Microsoft.Diagnostics.NETCore.Client": { + "type": "Transitive", + "resolved": "0.2.510501", + "contentHash": "juoqJYMDs+lRrrZyOkXXMImJHneCF23cuvO4waFRd2Ds7j+ZuGIPbJm0Y/zz34BdeaGiiwGWraMUlln05W1PCQ==", + "dependencies": { + "Microsoft.Extensions.Logging": "6.0.0" + } + }, + "Microsoft.Diagnostics.Runtime": { + "type": "Transitive", + "resolved": "3.1.512801", + "contentHash": "0lMUDr2oxNZa28D6NH5BuSQEe5T9tZziIkvkD44YkkCGQXPJqvFjLq5ZQq1hYLl3RjQJrY+hR0jFgap+EWPDTw==", + "dependencies": { + "Microsoft.Diagnostics.NETCore.Client": "0.2.410101" + } + }, + "Microsoft.Diagnostics.Tracing.TraceEvent": { + "type": "Transitive", + "resolved": "3.1.21", + "contentHash": "/OrJFKaojSR6TkUKtwh8/qA9XWNtxLrXMqvEb89dBSKCWjaGVTbKMYodIUgF5deCEtmd6GXuRerciXGl5bhZ7Q==", + "dependencies": { + "Microsoft.Diagnostics.NETCore.Client": "0.2.510501", + "System.Reflection.TypeExtensions": "4.7.0" + } + }, + "Microsoft.DotNet.PlatformAbstractions": { + "type": "Transitive", + "resolved": "3.1.6", + "contentHash": "jek4XYaQ/PGUwDKKhwR8K47Uh1189PFzMeLqO83mXrXQVIpARZCcfuDedH50YDTepBkfijCZN5U/vZi++erxtg==" + }, + "Microsoft.Extensions.DependencyInjection": { + "type": "Transitive", + "resolved": "6.0.0", + "contentHash": "k6PWQMuoBDGGHOQTtyois2u4AwyVcIwL2LaSLlTZQm2CYcJ1pxbt6jfAnpWmzENA/wfrYRI/X9DTLoUkE4AsLw==", + "dependencies": { + "Microsoft.Extensions.DependencyInjection.Abstractions": "6.0.0" + } + }, + "Microsoft.Extensions.DependencyInjection.Abstractions": { + "type": "Transitive", + "resolved": "6.0.0", + "contentHash": "xlzi2IYREJH3/m6+lUrQlujzX8wDitm4QGnUu6kUXTQAWPuZY8i+ticFJbzfqaetLA6KR/rO6Ew/HuYD+bxifg==" + }, + "Microsoft.Extensions.Logging": { + "type": "Transitive", + "resolved": "6.0.0", + "contentHash": "eIbyj40QDg1NDz0HBW0S5f3wrLVnKWnDJ/JtZ+yJDFnDj90VoPuoPmFkeaXrtu+0cKm5GRAwoDf+dBWXK0TUdg==", + "dependencies": { + "Microsoft.Extensions.DependencyInjection": "6.0.0", + "Microsoft.Extensions.DependencyInjection.Abstractions": "6.0.0", + "Microsoft.Extensions.Logging.Abstractions": "6.0.0", + "Microsoft.Extensions.Options": "6.0.0" + } + }, + "Microsoft.Extensions.Logging.Abstractions": { + "type": "Transitive", + "resolved": "6.0.0", + "contentHash": "/HggWBbTwy8TgebGSX5DBZ24ndhzi93sHUBDvP1IxbZD7FDokYzdAr6+vbWGjw2XAfR2EJ1sfKUotpjHnFWPxA==" + }, + "Microsoft.Extensions.Options": { + "type": "Transitive", + "resolved": "6.0.0", + "contentHash": "dzXN0+V1AyjOe2xcJ86Qbo233KHuLEY0njf/P2Kw8SfJU+d45HNS2ctJdnEnrWbM9Ye2eFgaC5Mj9otRMU6IsQ==", + "dependencies": { + "Microsoft.Extensions.DependencyInjection.Abstractions": "6.0.0", + "Microsoft.Extensions.Primitives": "6.0.0" + } + }, + "Microsoft.Extensions.Primitives": { + "type": "Transitive", + "resolved": "6.0.0", + "contentHash": "9+PnzmQFfEFNR9J2aDTfJGGupShHjOuGw4VUv+JB044biSHrnmCIMD+mJHmb2H7YryrfBEXDurxQ47gJZdCKNQ==" + }, + "Perfolizer": { + "type": "Transitive", + "resolved": "0.6.1", + "contentHash": "CR1QmWg4XYBd1Pb7WseP+sDmV8nGPwvmowKynExTqr3OuckIGVMhvmN4LC5PGzfXqDlR295+hz/T7syA1CxEqA==", + "dependencies": { + "Pragmastat": "3.2.4" + } + }, + "Pragmastat": { + "type": "Transitive", + "resolved": "3.2.4", + "contentHash": "I5qFifWw/gaTQT52MhzjZpkm/JPlfjSeO/DTZJjO7+hTKI+0aGRgOgZ3NN6D96dDuuqbIAZSeA5RimtHjqrA2A==" + }, + "System.CodeDom": { + "type": "Transitive", + "resolved": "9.0.5", + "contentHash": "cuzLM2MWutf9ZBEMPYYfd0DXwYdvntp7VCT6a/wvbKCa2ZuvGmW74xi+YBa2mrfEieAXqM4TNKlMmSnfAfpUoQ==" + }, + "System.Management": { + "type": "Transitive", + "resolved": "9.0.5", + "contentHash": "n6o9PZm9p25+zAzC3/48K0oHnaPKTInRrxqFq1fi/5TPbMLjuoCm/h//mS3cUmSy+9AO1Z+qsC/Ilt/ZFatv5Q==", + "dependencies": { + "System.CodeDom": "9.0.5" + } + }, + "System.Reflection.TypeExtensions": { + "type": "Transitive", + "resolved": "4.7.0", + "contentHash": "VybpaOQQhqE6siHppMktjfGBw1GCwvCqiufqmP8F1nj7fTUNtW35LOEt3UZTEsECfo+ELAl/9o9nJx3U91i7vA==" + }, + "just.precisemath": { + "type": "Project" + } + } + } +} \ No newline at end of file diff --git a/Just.PreciseMath.slnx b/Just.PreciseMath.slnx index 6c3926c..b4ce7ac 100644 --- a/Just.PreciseMath.slnx +++ b/Just.PreciseMath.slnx @@ -7,6 +7,7 @@ +