fixed multiplication asymmetry

0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs

@@ -91,7 +91,10 @@ public readonly partial struct DoubleDouble :
     }
 
     /// <summary>Multiplies expansions using an FMA product residual and cross terms.</summary>
-    /// <remarks>Results are approximate double-double values, not universally correctly rounded.</remarks>
+    /// <remarks>
+    /// Swapping operands preserves both result components bit-for-bit.
+    /// Results are approximate double-double values, not universally correctly rounded.
+    /// </remarks>
     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);

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<DoubleDouble> 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);
+                    // Operand ordering intentionally changes the old unordered
-                    error = Math.FusedMultiplyAdd(left.High, right.Low, error);
+                    // result bits, not the FMA kernel or its normalization.
-                    error = Math.FusedMultiplyAdd(left.Low, right.High, error);
+                    DoubleDouble first = left;
-                    error = Math.FusedMultiplyAdd(left.Low, right.Low, error);
+                    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;

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<DoubleDouble> 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<double> scalars = [0.0, -0.0, 1.0, -1.0, 2.0, -2.0, 3.0, -3.0, 100.0, -100.0,