From 34076c246605c09baa45869c30bff156f338e726 Mon Sep 17 00:00:00 2001 From: just Date: Sat, 19 Sep 2026 20:02:33 +0400 Subject: [PATCH] fixed multiplication asymmetry --- .../DoubleDouble.Arithmetic.cs | 13 ++- .../ArithmeticRangeTests.cs | 19 +++-- .../DoubleDoubleArithmeticTests.cs | 79 +++++++++++++++++++ 3 files changed, 105 insertions(+), 6 deletions(-) diff --git a/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs b/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs index e66ae4e..ea5d214 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs @@ -91,7 +91,10 @@ public readonly partial struct DoubleDouble : } /// Multiplies expansions using an FMA product residual and cross terms. - /// Results are approximate double-double values, not universally correctly rounded. + /// + /// Swapping operands preserves both result components bit-for-bit. + /// Results are approximate double-double values, not universally correctly rounded. + /// public static DoubleDouble operator *(DoubleDouble left, DoubleDouble right) { if (!IsFinite(left) || !IsFinite(right) || left._high == 0.0 || right._high == 0.0) @@ -103,6 +106,14 @@ public readonly partial struct DoubleDouble : return PreciseMathHelper.MultiplyBoundary(left, right); } + // Sequential cross-term FMAs round in order. Use the same operand order + // in both directions, including when highs tie. Inputs here are finite + // and nonzero with canonical lows; the exact boundary path is symmetric. + if (left._high > right._high || (left._high == right._high && left._low > right._low)) + { + (left, right) = (right, left); + } + (double product, double error) = PreciseMathHelper.TwoMultiply(left._high, right._high); error = Math.FusedMultiplyAdd(left._high, right._low, error); error = Math.FusedMultiplyAdd(left._low, right._high, error); diff --git a/1-tests/Just.PreciseMath.Tests/ArithmeticRangeTests.cs b/1-tests/Just.PreciseMath.Tests/ArithmeticRangeTests.cs index c500019..8cf9653 100644 --- a/1-tests/Just.PreciseMath.Tests/ArithmeticRangeTests.cs +++ b/1-tests/Just.PreciseMath.Tests/ArithmeticRangeTests.cs @@ -81,7 +81,8 @@ public class ArithmeticRangeTests public void FiniteKernelsPreservePreviousComponentBits(string operation) { // Differential characterization, not an independent accuracy oracle. - // Freeze the previous four-TwoSum/FMA expressions and public normalization; + // Freeze the previous four-TwoSum/FMA expressions and public normalization, + // with multiplication evaluated in its deterministic operand order; // independent exact/rational accuracy cases remain in the arithmetic suites. List values = [new(0.0), new(-0.0)]; int[] exponents = [-1074, -1022, -901, -900, -899, -451, -450, -1, 0, 1, @@ -154,10 +155,18 @@ public class ArithmeticRangeTests { continue; } - (double product, double error) = PreciseMathHelper.TwoMultiply(left.High, right.High); - error = Math.FusedMultiplyAdd(left.High, right.Low, error); - error = Math.FusedMultiplyAdd(left.Low, right.High, error); - error = Math.FusedMultiplyAdd(left.Low, right.Low, error); + // Operand ordering intentionally changes the old unordered + // result bits, not the FMA kernel or its normalization. + DoubleDouble first = left; + DoubleDouble second = right; + if (first.High > second.High || (first.High == second.High && first.Low > second.Low)) + { + (first, second) = (second, first); + } + (double product, double error) = PreciseMathHelper.TwoMultiply(first.High, second.High); + error = Math.FusedMultiplyAdd(first.High, second.Low, error); + error = Math.FusedMultiplyAdd(first.Low, second.High, error); + error = Math.FusedMultiplyAdd(first.Low, second.Low, error); Math.Abs(error).ShouldBeLessThan(Math.Abs(product)); expected = DoubleDouble.FromComponents(product, error); actual = left * right; diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs index 4208975..0a26ad7 100644 --- a/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs @@ -33,6 +33,56 @@ public class DoubleDoubleArithmeticTests CheckBoundary(left / b, quotient, 0.0); } + [Theory] + [InlineData(1.0, 5.551115123125784e-17, 1.25, -5.551115123125783e-17)] + [InlineData(1.25, -5.551115123125783e-17, 1.25, -5.551115123125784e-17)] + public void MultiplicationIsBitwiseCommutativeWhenCrossTermRoundingDependsOnOrder( + double leftHigh, double leftLow, double rightHigh, double rightLow) + { + // The old sequential FMA accumulation rounded differently when reversed. + // Equal highs in the second row require a low-component ordering tie-break. + foreach (double leftSign in new[] { -1.0, 1.0 }) + { + foreach (double rightSign in new[] { -1.0, 1.0 }) + { + DoubleDouble left = DoubleDouble.FromComponents(leftSign * leftHigh, leftSign * leftLow); + DoubleDouble right = DoubleDouble.FromComponents(rightSign * rightHigh, rightSign * rightLow); + DoubleDouble forward = left * right; + DoubleDouble reverse = right * left; + CheckBoundary(reverse, forward.High, forward.Low); + // Symmetry is a compatibility contract, not an accuracy oracle. + BigInteger numerator = Units(left) * Units(right); + BigInteger denominator = BigInteger.One << 1074; + AssertRelative(forward, numerator, denominator); + AssertRelative(reverse, numerator, denominator); + } + } + } + + [Fact] + public void MultiplicationIsBitwiseCommutativeAcrossSpecialValuesAndRangeBoundaries() + { + DoubleDouble[] values = MultiplicationCommutativityCases(); + foreach (DoubleDouble left in values) + { + foreach (DoubleDouble right in values) + { + DoubleDouble forward = left * right; + DoubleDouble reverse = right * left; + CheckBoundary(reverse, forward.High, forward.Low); + DoubleDouble.IsCanonical(forward).ShouldBeTrue(); + if (!DoubleDouble.IsFinite(left) || !DoubleDouble.IsFinite(right) + || left.High == 0.0 || right.High == 0.0) + { + // Independent binary64 special-value rule, including zero signs + // and the library's canonical NaN / positive-zero low contract. + double high = left.High * right.High; + CheckBoundary(forward, double.IsNaN(high) ? double.NaN : high, 0.0); + } + } + } + } + [Theory] [InlineData(30.0, 0.5235987755982989, -5.360408832255455e-17)] [InlineData(45.0, 0.7853981633974483, 3.061616997868383e-17)] @@ -254,6 +304,12 @@ public class DoubleDoubleArithmeticTests ("/", "double/DD") => left.High / right, _ => throw new ArgumentOutOfRangeException(nameof(operation)) }; + if (operation == "*" && overload == "DD/DD") + { + // Pin symmetry for the seeded dense/sparse, cancellation and range + // matrix separately from its independent exact-rational accuracy bound. + CheckBoundary(right * left, actual.High, actual.Low); + } // Units decodes both IEEE-754 components independently of library arithmetic. // Expected numerator/denominator is measured in units of epsilon, not doubles. BigInteger x = overload == "double/DD" ? Units(left.High) : Units(left); @@ -1086,6 +1142,29 @@ public class DoubleDoubleArithmeticTests CheckBoundary(DDMath.Square(value), actual.High, actual.Low); } + private static DoubleDouble[] MultiplicationCommutativityCases() + { + List values = [new(0.0), new(-0.0), new(double.PositiveInfinity), + new(double.NegativeInfinity), new(double.NaN), new(double.MaxValue), new(double.MinValue)]; + // Pair exponent edges with +/-1 and one another to cross the fast-path + // exponent-sum bounds, underflow and overflow. Repeated highs exercise + // low tie-breaks; a value paired with itself exercises exact operand ties. + int[] exponents = [-1074, -1022, -901, -900, -899, -1, 0, 1, 899, 900, 901, 1023]; + foreach (int exponent in exponents) + { + foreach (double sign in new[] { -1.0, 1.0 }) + { + double high = Math.ScaleB(sign, exponent); + double denseLow = Math.ScaleB(1.0, exponent - 54); + foreach (double low in new[] { 0.0, denseLow, -denseLow, double.Epsilon, -double.Epsilon }) + { + values.Add(DoubleDouble.FromComponents(high, low)); + } + } + } + return values.ToArray(); + } + private static (double[] Scalars, DoubleDouble[] Values) ScalarArithmeticCompatibilityCases() { List scalars = [0.0, -0.0, 1.0, -1.0, 2.0, -2.0, 3.0, -3.0, 100.0, -100.0,