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,