diff --git a/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs b/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs
index 221f2f0..1afd017 100644
--- a/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs
+++ b/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs
@@ -35,22 +35,30 @@ public readonly partial struct DoubleDouble :
{
return new DoubleDouble(left._high + right._high);
}
- if (Math.Max(Math.ILogB(left._high), Math.ILogB(right._high)) > 1020)
+ if (!PreciseMathHelper.IsAdditionWithinFastRange(left._high) || !PreciseMathHelper.IsAdditionWithinFastRange(right._high))
{
- return PreciseMathHelper.ArithmeticFromRatio(PreciseMathHelper.ArithmeticUnits(left) + PreciseMathHelper.ArithmeticUnits(right), BigInteger.One << 1074);
+ return PreciseMathHelper.AddBoundary(left, right);
}
- (double high, double highError) = PreciseMathHelper.TwoAdd(left._high, right._high);
- (double low, double lowError) = PreciseMathHelper.TwoAdd(left._low, right._low);
- (double middle, double middleError) = PreciseMathHelper.TwoAdd(highError, low);
- (double sum, double sumError) = PreciseMathHelper.TwoAdd(high, middle);
- return FromComponents(sum, sumError + (middleError + lowError));
+ return PreciseMathHelper.AddFinite(left._high, left._low, right._high, right._low);
}
/// Subtracts normalized expansions.
public static DoubleDouble operator -(DoubleDouble left, DoubleDouble right)
{
- return left + (-right);
+ if (!IsFinite(left) || !IsFinite(right) || (left._high == 0.0 && right._high == 0.0))
+ {
+ return new DoubleDouble(left._high - right._high);
+ }
+ if (!PreciseMathHelper.IsAdditionWithinFastRange(left._high) || !PreciseMathHelper.IsAdditionWithinFastRange(right._high))
+ {
+ return PreciseMathHelper.AddBoundary(left, -right);
+ }
+
+ // With a canonical left low (never -0), the low TwoSum absorbs the
+ // negated right zero low without changing either output component's bits.
+ // Avoid unary negation's intermediate canonicalization on this finite path.
+ return PreciseMathHelper.AddFinite(left._high, left._low, -right._high, -right._low);
}
/// Multiplies expansions using an FMA product residual and cross terms.
@@ -61,17 +69,20 @@ public readonly partial struct DoubleDouble :
{
return new DoubleDouble(left._high * right._high);
}
- int exponent = Math.ILogB(left._high) + Math.ILogB(right._high);
- if (exponent < -900 || exponent > 900)
+ if (!PreciseMathHelper.IsMultiplicationWithinFastRange(left._high, right._high))
{
- return PreciseMathHelper.ArithmeticFromRatio(PreciseMathHelper.ArithmeticUnits(left) * PreciseMathHelper.ArithmeticUnits(right), BigInteger.One << 2148);
+ return PreciseMathHelper.MultiplyBoundary(left, right);
}
(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);
- return FromComponents(product, error);
+ // With u = 2^-53, normalized inputs and the exponent-sum guard give
+ // |error| < 4u*|product|, including rounding at the subnormal floor.
+ // The product is normal and nonzero; its corrected sum remains finite.
+ (double high, double low) = PreciseMathHelper.TwoQuickAdd(product, error);
+ return new DoubleDouble(high, low == 0.0 ? 0.0 : low);
}
/// Divides expansions using a quotient estimate and two residual corrections.
@@ -82,11 +93,9 @@ public readonly partial struct DoubleDouble :
{
return new DoubleDouble(left._high / right._high);
}
- int leftExponent = Math.ILogB(left._high);
- int rightExponent = Math.ILogB(right._high);
- if (Math.Abs(leftExponent) > 450 || Math.Abs(rightExponent) > 450)
+ if (!PreciseMathHelper.IsDivisionWithinFastRange(left._high) || !PreciseMathHelper.IsDivisionWithinFastRange(right._high))
{
- return PreciseMathHelper.ArithmeticFromRatio(PreciseMathHelper.ArithmeticUnits(left), PreciseMathHelper.ArithmeticUnits(right));
+ return PreciseMathHelper.DivideBoundary(left, right);
}
double quotient = left._high / right._high;
@@ -129,10 +138,9 @@ public readonly partial struct DoubleDouble :
{
return new DoubleDouble(left._high * right);
}
- int exponent = Math.ILogB(left._high) + Math.ILogB(right);
- if (exponent < -900 || exponent > 900)
+ if (!PreciseMathHelper.IsMultiplicationWithinFastRange(left._high, right))
{
- return PreciseMathHelper.ArithmeticFromRatio(PreciseMathHelper.ArithmeticUnits(left) * PreciseMathHelper.ArithmeticUnits(right), BigInteger.One << 2148);
+ return PreciseMathHelper.MultiplyBoundary(left, right);
}
(double product, double error) = PreciseMathHelper.TwoMultiply(left._high, right);
@@ -156,11 +164,9 @@ public readonly partial struct DoubleDouble :
{
return new DoubleDouble(left._high / right);
}
- int leftExponent = Math.ILogB(left._high);
- int rightExponent = Math.ILogB(right);
- if (Math.Abs(leftExponent) > 450 || Math.Abs(rightExponent) > 450)
+ if (!PreciseMathHelper.IsDivisionWithinFastRange(left._high) || !PreciseMathHelper.IsDivisionWithinFastRange(right))
{
- return PreciseMathHelper.ArithmeticFromRatio(PreciseMathHelper.ArithmeticUnits(left), PreciseMathHelper.ArithmeticUnits(right));
+ return PreciseMathHelper.DivideBoundary(left, right);
}
double quotient = left._high / right;
@@ -179,11 +185,9 @@ public readonly partial struct DoubleDouble :
{
return new DoubleDouble(left / right._high);
}
- int leftExponent = Math.ILogB(left);
- int rightExponent = Math.ILogB(right._high);
- if (Math.Abs(leftExponent) > 450 || Math.Abs(rightExponent) > 450)
+ if (!PreciseMathHelper.IsDivisionWithinFastRange(left) || !PreciseMathHelper.IsDivisionWithinFastRange(right._high))
{
- return PreciseMathHelper.ArithmeticFromRatio(PreciseMathHelper.ArithmeticUnits(left), PreciseMathHelper.ArithmeticUnits(right));
+ return PreciseMathHelper.DivideBoundary(left, right);
}
double quotient = left / right._high;
diff --git a/0-source/Just.PreciseMath/PreciseMathHelper.cs b/0-source/Just.PreciseMath/PreciseMathHelper.cs
index 5c909d3..5532e03 100644
--- a/0-source/Just.PreciseMath/PreciseMathHelper.cs
+++ b/0-source/Just.PreciseMath/PreciseMathHelper.cs
@@ -2,6 +2,45 @@ namespace Just.PreciseMath;
internal static class PreciseMathHelper
{
+ // Finite operands only; callers retain their special-value/zero handling.
+ // ILogB(value) <= 1020 is exactly biasedExponent <= 2043. Zero and
+ // subnormal operands also qualify without computing their true exponents.
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ internal static bool IsAdditionWithinFastRange(double value)
+ {
+ int exponent = (int)((BitConverter.DoubleToUInt64Bits(value) >> 52) & 0x7ff);
+ return exponent <= 2043;
+ }
+
+ // Finite nonzero operands only. The original inclusive ILogB interval
+ // [-450, 450] becomes [573, 1473] with the binary64 bias of 1023.
+ // Unsigned subtraction rejects smaller exponents, including subnormals.
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ internal static bool IsDivisionWithinFastRange(double value)
+ {
+ int exponent = (int)((BitConverter.DoubleToUInt64Bits(value) >> 52) & 0x7ff);
+ return unchecked((uint)(exponent - 573)) <= 900;
+ }
+
+ // Finite nonzero operands only. For normal operands the exponent sum
+ // interval [-900, 900] becomes [1146, 2946] after adding both biases.
+ // A subnormal times a large normal can still be in range: preserve the
+ // original ILogB calculation for those operands, not an allocating detour.
+ // Both multiplication overloads reject nonfinite and zero operands first;
+ // this predicate does not validate them (ILogB(0) is an integer sentinel).
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ internal static bool IsMultiplicationWithinFastRange(double left, double right)
+ {
+ int leftExponent = (int)((BitConverter.DoubleToUInt64Bits(left) >> 52) & 0x7ff);
+ int rightExponent = (int)((BitConverter.DoubleToUInt64Bits(right) >> 52) & 0x7ff);
+ if (leftExponent == 0 || rightExponent == 0)
+ {
+ int exponent = Math.ILogB(left) + Math.ILogB(right);
+ return exponent >= -900 && exponent <= 900;
+ }
+ return unchecked((uint)((leftExponent + rightExponent) - 1146)) <= 1800;
+ }
+
// General TwoSum: no magnitude ordering required, but inputs, sum, and
// intermediate subtractions must stay finite. Arithmetic callers bound the
// exponents; arbitrary-component normalization uses magnitude ordering instead.
@@ -47,6 +86,44 @@ internal static class PreciseMathHelper
return (r, Math.FusedMultiplyAdd(a, a, -r));
}
+ // Normalized finite operands passing the addition range guard, excluding two
+ // zero highs. The left low must be canonical (nonzero or +0, never -0).
+ // The right low may be a negated zero: under that left-low precondition,
+ // TwoAdd(leftLow, +0) and TwoAdd(leftLow, -0) have bit-identical outputs.
+ // Keep all four transforms and grouping.
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ internal static DoubleDouble AddFinite(double leftHigh, double leftLow, double rightHigh, double rightLow)
+ {
+ (double high, double highError) = TwoAdd(leftHigh, rightHigh);
+ (double low, double lowError) = TwoAdd(leftLow, rightLow);
+ (double middle, double middleError) = TwoAdd(highError, low);
+ (double sum, double sumError) = TwoAdd(high, middle);
+ // Input highs are < 2^1021 and lows <= 2^967 in magnitude. These
+ // transforms keep |sum| <= 2^1022 and the correction <= 2^970,
+ // so finite-only normalization is safe even under cancellation.
+ return NormalizeFinite(sum, sumError + (middleError + lowError));
+ }
+
+ // Both components and their rounded sum must be finite. Unlike QuickTwoSum,
+ // this entry point permits either magnitude order, including cancellation.
+ // AddFinite establishes 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)]
+ internal static DoubleDouble NormalizeFinite(double high, double low)
+ {
+ if (low == 0.0)
+ {
+ return new DoubleDouble(high, 0.0);
+ }
+
+ double sum = high + low;
+ double error = Math.Abs(high) >= Math.Abs(low)
+ ? low - (sum - high)
+ : high - (sum - low);
+ return new DoubleDouble(sum, error == 0.0 ? 0.0 : error);
+ }
+
// The first two arguments are normalized components (a negated zero low is
// also allowed). Sharing this path preserves both subtraction orders without
// constructing a temporary expansion for the scalar or the negated operand.
@@ -56,9 +133,9 @@ internal static class PreciseMathHelper
{
return new DoubleDouble(high + value);
}
- if (Math.Max(Math.ILogB(high), Math.ILogB(value)) > 1020)
+ if (!IsAdditionWithinFastRange(high) || !IsAdditionWithinFastRange(value))
{
- return ArithmeticFromRatio(ArithmeticUnits(high) + ArithmeticUnits(low) + ArithmeticUnits(value), BigInteger.One << 1074);
+ return AddScalarBoundary(high, low, value);
}
(double sum, double error) = TwoAdd(high, value);
@@ -69,6 +146,54 @@ internal static class PreciseMathHelper
return new DoubleDouble(result, residual == 0.0 ? 0.0 : residual);
}
+ // Keep the complete BigInteger expressions out of ordinary arithmetic bodies,
+ // including operand conversion and denominator construction. NoInlining isolates
+ // this setup even when the public operators are inlined by their callers.
+ // All inputs must be finite and normalized; division denominators must be nonzero.
+ // Callers retain the special-value, signed-zero, and exponent-range dispatch.
+ [MethodImpl(MethodImplOptions.NoInlining)]
+ internal static DoubleDouble AddBoundary(DoubleDouble left, DoubleDouble right)
+ {
+ return ArithmeticFromRatio(ArithmeticUnits(left) + ArithmeticUnits(right), BigInteger.One << 1074);
+ }
+
+ [MethodImpl(MethodImplOptions.NoInlining)]
+ private static DoubleDouble AddScalarBoundary(double high, double low, double value)
+ {
+ // As in AddScalar, high/low are normalized but a negated zero low is allowed.
+ return ArithmeticFromRatio(ArithmeticUnits(high) + ArithmeticUnits(low) + ArithmeticUnits(value), BigInteger.One << 1074);
+ }
+
+ [MethodImpl(MethodImplOptions.NoInlining)]
+ internal static DoubleDouble MultiplyBoundary(DoubleDouble left, DoubleDouble right)
+ {
+ return ArithmeticFromRatio(ArithmeticUnits(left) * ArithmeticUnits(right), BigInteger.One << 2148);
+ }
+
+ [MethodImpl(MethodImplOptions.NoInlining)]
+ internal static DoubleDouble MultiplyBoundary(DoubleDouble left, double right)
+ {
+ return ArithmeticFromRatio(ArithmeticUnits(left) * ArithmeticUnits(right), BigInteger.One << 2148);
+ }
+
+ [MethodImpl(MethodImplOptions.NoInlining)]
+ internal static DoubleDouble DivideBoundary(DoubleDouble left, DoubleDouble right)
+ {
+ return ArithmeticFromRatio(ArithmeticUnits(left), ArithmeticUnits(right));
+ }
+
+ [MethodImpl(MethodImplOptions.NoInlining)]
+ internal static DoubleDouble DivideBoundary(DoubleDouble left, double right)
+ {
+ return ArithmeticFromRatio(ArithmeticUnits(left), ArithmeticUnits(right));
+ }
+
+ [MethodImpl(MethodImplOptions.NoInlining)]
+ internal static DoubleDouble DivideBoundary(double left, DoubleDouble right)
+ {
+ return ArithmeticFromRatio(ArithmeticUnits(left), ArithmeticUnits(right));
+ }
+
// The boundary path uses bounded binary integers (at most about 4200 bits), not
// arbitrary-precision storage. It avoids overflow and double rounding in EFTs
// at the binary64 exponent limits. The common path remains allocation-free.
diff --git a/1-tests/Just.PreciseMath.Tests/ArithmeticRangeTests.cs b/1-tests/Just.PreciseMath.Tests/ArithmeticRangeTests.cs
new file mode 100644
index 0000000..7fedfba
--- /dev/null
+++ b/1-tests/Just.PreciseMath.Tests/ArithmeticRangeTests.cs
@@ -0,0 +1,195 @@
+using Shouldly;
+using Xunit;
+
+namespace Just.PreciseMath.Tests;
+
+public class ArithmeticRangeTests
+{
+ [Fact]
+ public void AdditionRangeMatchesTheOriginalExponentGuard()
+ {
+ // Independent BCL predicate: include both signed zeros, subnormals, and
+ // the first/last significands of every finite exponent field.
+ foreach (double value in FiniteExponentSamples())
+ {
+ bool expected = Math.ILogB(value) <= 1020;
+ PreciseMathHelper.IsAdditionWithinFastRange(value).ShouldBe(expected, $"{value:R}");
+ }
+ }
+
+ [Fact]
+ public void DivisionRangeMatchesTheOriginalExponentGuard()
+ {
+ // Public operators handle zero before this guard; ILogB(0) is the
+ // int.MinValue sentinel, whose absolute value is not representable.
+ foreach (double value in FiniteExponentSamples())
+ {
+ if (value != 0.0)
+ {
+ bool expected = Math.Abs(Math.ILogB(value)) <= 450;
+ PreciseMathHelper.IsDivisionWithinFastRange(value).ShouldBe(expected, $"{value:R}");
+ }
+ }
+ }
+
+ [Fact]
+ public void MultiplicationRangeMatchesEveryNormalExponentPair()
+ {
+ // Exhaust all normal exponent combinations. The original guard depends
+ // only on these exponents, never on the sign or fractional significand.
+ double[] values = new double[2046];
+ for (int i = 0; i < values.Length; ++i)
+ {
+ values[i] = Math.ScaleB(1.0, i - 1022);
+ }
+ foreach (double left in values)
+ {
+ foreach (double right in values)
+ {
+ AssertMultiplicationRange(left, right);
+ }
+ }
+ }
+
+ [Fact]
+ public void MultiplicationRangePreservesSubnormalAndSignedOperandPaths()
+ {
+ // Include all subnormal binades and both significand edges, paired
+ // with large normals that can bring the exponent sum into range.
+ double[] partners = [double.Epsilon, Math.BitDecrement(Math.ScaleB(1.0, -1022)),
+ Math.ScaleB(1.0, -1022), Math.ScaleB(1.0, -901), Math.ScaleB(1.0, -900),
+ Math.ScaleB(1.0, -1), 1.0, 2.0, Math.ScaleB(1.0, 900), Math.ScaleB(1.0, 901),
+ Math.ScaleB(1.0, 1023), double.MaxValue];
+ foreach (double value in FiniteExponentSamples())
+ {
+ if (value == 0.0)
+ {
+ continue;
+ }
+ foreach (double partner in partners)
+ {
+ AssertMultiplicationRange(value, partner);
+ AssertMultiplicationRange(partner, value);
+ AssertMultiplicationRange(value, -partner);
+ AssertMultiplicationRange(-partner, value);
+ }
+ }
+ }
+
+ [Theory]
+ [InlineData("+")]
+ [InlineData("-")]
+ [InlineData("*")]
+ public void FiniteKernelsPreservePreviousComponentBits(string operation)
+ {
+ // Differential characterization, not an independent accuracy oracle.
+ // Freeze the previous four-TwoSum/FMA expressions and public normalization;
+ // 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,
+ 450, 451, 899, 900, 901, 1020, 1021, 1023];
+ foreach (int exponent in exponents)
+ {
+ foreach (double significand in new[] { 1.0, Math.BitIncrement(1.0), 1.5, Math.BitDecrement(2.0) })
+ {
+ foreach (double sign in new[] { -1.0, 1.0 })
+ {
+ double high = sign * Math.ScaleB(significand, exponent);
+ foreach (double low in new[] { 0.0, Math.ScaleB(high, -53), -Math.ScaleB(high, -53),
+ double.Epsilon, -double.Epsilon })
+ {
+ DoubleDouble value = DoubleDouble.FromComponents(high, low);
+ if (DoubleDouble.IsFinite(value))
+ {
+ values.Add(value);
+ }
+ }
+ }
+ }
+ }
+ foreach (DoubleDouble left in values)
+ {
+ foreach (DoubleDouble right in values)
+ {
+ DoubleDouble expected;
+ DoubleDouble actual;
+ if (operation == "*")
+ {
+ if (left.High == 0.0 || right.High == 0.0)
+ {
+ continue;
+ }
+ int exponent = Math.ILogB(left.High) + Math.ILogB(right.High);
+ if (exponent < -900 || exponent > 900)
+ {
+ 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);
+ Math.Abs(error).ShouldBeLessThan(Math.Abs(product));
+ expected = DoubleDouble.FromComponents(product, error);
+ actual = left * right;
+ }
+ else
+ {
+ if (Math.ILogB(left.High) > 1020 || Math.ILogB(right.High) > 1020)
+ {
+ continue;
+ }
+ DoubleDouble operand = operation == "-" ? -right : right;
+ if (left.High == 0.0 && operand.High == 0.0)
+ {
+ expected = new DoubleDouble(left.High + operand.High);
+ }
+ else
+ {
+ (double high, double highError) = PreciseMathHelper.TwoAdd(left.High, operand.High);
+ (double low, double lowError) = PreciseMathHelper.TwoAdd(left.Low, operand.Low);
+ (double middle, double middleError) = PreciseMathHelper.TwoAdd(highError, low);
+ (double sum, double sumError) = PreciseMathHelper.TwoAdd(high, middle);
+ double correction = sumError + (middleError + lowError);
+ double.IsFinite(sum + correction).ShouldBeTrue();
+ expected = DoubleDouble.FromComponents(sum, correction);
+ }
+ actual = operation == "-" ? left - right : left + right;
+ }
+ BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(BitConverter.DoubleToInt64Bits(expected.High));
+ BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(BitConverter.DoubleToInt64Bits(expected.Low));
+ }
+ }
+ }
+
+ private static void AssertMultiplicationRange(double left, double right)
+ {
+ int exponent = Math.ILogB(left) + Math.ILogB(right);
+ bool expected = exponent >= -900 && exponent <= 900;
+ bool actual = PreciseMathHelper.IsMultiplicationWithinFastRange(left, right);
+ // Only format diagnostics on failure in this exhaustive matrix.
+ if (actual != expected)
+ {
+ actual.ShouldBe(expected, $"({left:R}, {right:R})");
+ }
+ }
+
+ private static IEnumerable FiniteExponentSamples()
+ {
+ ulong[] fractions = [0, 1, 0x0008_0000_0000_0000, 0x000f_ffff_ffff_ffff];
+ foreach (ulong sign in new[] { 0UL, 0x8000_0000_0000_0000UL })
+ {
+ for (int bit = 0; bit < 52; ++bit)
+ {
+ yield return BitConverter.UInt64BitsToDouble(sign | (1UL << bit));
+ yield return BitConverter.UInt64BitsToDouble(sign | ((1UL << (bit + 1)) - 1));
+ }
+ for (ulong exponent = 0; exponent < 0x7ff; ++exponent)
+ {
+ foreach (ulong fraction in fractions)
+ {
+ yield return BitConverter.UInt64BitsToDouble(sign | (exponent << 52) | fraction);
+ }
+ }
+ }
+ }
+}
diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs
index b7bec28..104afab 100644
--- a/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs
+++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs
@@ -35,6 +35,53 @@ public class DoubleDoubleArithmeticTests
Check(6.0 / new DoubleDouble(2.0), 3.0, 0.0);
}
+ [Fact]
+ public void ExpansionCancellationRetainsExactComponentsAcrossTheAdditionGuard()
+ {
+ // (2^e + 2^(e-54)) - (2^e - 2^(e-108)) is exactly the
+ // normalized pair (2^(e-54), 2^(e-108)). The smallest residual
+ // is epsilon; the largest case exercises the boundary fallback.
+ foreach (int exponent in new[] { -966, -450, 0, 450, 1020, 1021 })
+ {
+ foreach (double sign in new[] { -1.0, 1.0 })
+ {
+ double high = sign * Math.ScaleB(1.0, exponent);
+ double small = sign * Math.ScaleB(1.0, exponent - 54);
+ double tiny = sign * Math.ScaleB(1.0, exponent - 108);
+ DoubleDouble left = DoubleDouble.FromComponents(high, small);
+ DoubleDouble right = DoubleDouble.FromComponents(high, -tiny);
+ CheckBoundary(left - right, small, tiny);
+ CheckBoundary(right - left, -small, -tiny);
+ CheckBoundary(left + (-right), small, tiny);
+ CheckBoundary((-right) + left, small, tiny);
+ CheckBoundary(left - left, 0.0, 0.0);
+ }
+ }
+ }
+
+ [Fact]
+ public void ExpansionSubtractionHandlesSpecialValuesAlongsideNonzeroResiduals()
+ {
+ foreach (double sign in new[] { -1.0, 1.0 })
+ {
+ DoubleDouble value = DoubleDouble.FromComponents(sign, sign * double.Epsilon);
+ foreach (double special in new[] { 0.0, -0.0, double.NegativeInfinity, double.PositiveInfinity, double.NaN })
+ {
+ DoubleDouble other = new(special);
+ if (special == 0.0)
+ {
+ CheckBoundary(value - other, value.High, value.Low);
+ CheckBoundary(other - value, -value.High, -value.Low);
+ }
+ else
+ {
+ CheckBits(value - other, sign - special);
+ CheckBits(other - value, special - sign);
+ }
+ }
+ }
+ }
+
[Fact]
public void ScalarLeftSubtractionAppliesTheRequestedOperandOrder()
{
@@ -135,6 +182,41 @@ public class DoubleDoubleArithmeticTests
Check(scalar * value, high, low);
Check(value * (-scalar), -high, -low);
Check((-scalar) * value, -high, -low);
+ CheckBoundary(value * new DoubleDouble(scalar), high, low);
+ CheckBoundary(new DoubleDouble(scalar) * value, high, low);
+ CheckBoundary(value * new DoubleDouble(-scalar), -high, -low);
+ CheckBoundary(new DoubleDouble(-scalar) * value, -high, -low);
+ }
+
+ [Fact]
+ public void ExpansionProductRetainsTheLowLowTermAtFastRangeEndpoints()
+ {
+ // (1 + 2^-53)(1 - 2^-54) = 1 + 2^-54 - 2^-107 exactly.
+ // Its residual is BitDecrement(2^-54); omitting low*low loses that bit.
+ // Power-of-two scaling keeps both expected components representable,
+ // including exponent sums at each inclusive fast-path endpoint.
+ foreach (int leftExponent in new[] { -450, 0, 450 })
+ {
+ foreach (int rightExponent in new[] { -450, 0, 450 })
+ {
+ foreach (double leftSign in new[] { -1.0, 1.0 })
+ {
+ foreach (double rightSign in new[] { -1.0, 1.0 })
+ {
+ DoubleDouble left = DoubleDouble.FromComponents(leftSign * Math.ScaleB(1.0, leftExponent),
+ leftSign * Math.ScaleB(1.0, leftExponent - 53));
+ DoubleDouble right = DoubleDouble.FromComponents(rightSign * Math.ScaleB(1.0, rightExponent),
+ -rightSign * Math.ScaleB(1.0, rightExponent - 54));
+ int exponent = leftExponent + rightExponent;
+ double sign = leftSign * rightSign;
+ double high = sign * Math.ScaleB(1.0, exponent);
+ double low = sign * Math.ScaleB(Math.BitDecrement(Math.ScaleB(1.0, -54)), exponent);
+ CheckBoundary(left * right, high, low);
+ CheckBoundary(right * left, high, low);
+ }
+ }
+ }
+ }
}
[Fact]
@@ -150,6 +232,47 @@ public class DoubleDoubleArithmeticTests
AssertRelative(-1.0 / value, -numerator, Units(value));
}
+ [Fact]
+ public void DivisionRetainsSubnormalCorrectionsWithOrdinaryHighComponents()
+ {
+ // Dividing (1 + epsilon) by +/-1 is exact. The outer division is
+ // ordinary, but its second residual product can use the boundary path.
+ foreach (double sign in new[] { -1.0, 1.0 })
+ {
+ DoubleDouble numerator = DoubleDouble.FromComponents(sign, sign * double.Epsilon);
+ foreach (double denominator in new[] { -1.0, 1.0 })
+ {
+ double high = sign / denominator;
+ double low = high * double.Epsilon;
+ CheckBoundary(numerator / new DoubleDouble(denominator), high, low);
+ CheckBoundary(numerator / denominator, high, low);
+ }
+ }
+ }
+
+ [Fact]
+ public void SubnormalProductsWithLargeNormalsRetainExactResultsInBothOrders()
+ {
+ // 2^-1074 * 2^1023 = 2^-51 exactly, despite the subnormal input.
+ foreach (double leftSign in new[] { -1.0, 1.0 })
+ {
+ foreach (double rightSign in new[] { -1.0, 1.0 })
+ {
+ double tiny = leftSign * double.Epsilon;
+ double large = rightSign * Math.ScaleB(1.0, 1023);
+ double expected = (leftSign * rightSign) * Math.ScaleB(1.0, -51);
+ DoubleDouble left = new(tiny);
+ DoubleDouble right = new(large);
+ CheckBoundary(left * right, expected, 0.0);
+ CheckBoundary(right * left, expected, 0.0);
+ CheckBoundary(left * large, expected, 0.0);
+ CheckBoundary(large * left, expected, 0.0);
+ CheckBoundary(tiny * right, expected, 0.0);
+ CheckBoundary(right * tiny, expected, 0.0);
+ }
+ }
+ }
+
[Fact]
public void ExtremeFiniteOperationsDoNotOverflowIntermediates()
{
@@ -170,6 +293,29 @@ public class DoubleDoubleArithmeticTests
Check(below + Math.ScaleB(1.0, 969), double.MaxValue, 0.0);
}
+ [Fact]
+ public void SubnormalHighProductsRetainANormalPartnersResidual()
+ {
+ // (2^-1074, 0) * (2^1023, 2^969) = (2^-51, 2^-105), exactly.
+ // The high product is ordinary despite the subnormal input high.
+ foreach (double leftSign in new[] { -1.0, 1.0 })
+ {
+ foreach (double rightSign in new[] { -1.0, 1.0 })
+ {
+ double scalar = leftSign * double.Epsilon;
+ DoubleDouble tiny = new(scalar);
+ DoubleDouble large = DoubleDouble.FromComponents(rightSign * Math.ScaleB(1.0, 1023),
+ rightSign * Math.ScaleB(1.0, 969));
+ double high = (leftSign * rightSign) * Math.ScaleB(1.0, -51);
+ double low = (leftSign * rightSign) * Math.ScaleB(1.0, -105);
+ CheckBoundary(tiny * large, high, low);
+ CheckBoundary(large * tiny, high, low);
+ CheckBoundary(scalar * large, high, low);
+ CheckBoundary(large * scalar, high, low);
+ }
+ }
+ }
+
[Theory]
[InlineData(1.0)]
[InlineData(-1.0)]
diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleBoundaryTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleBoundaryTests.cs
index 1d0965a..0f7e90d 100644
--- a/1-tests/Just.PreciseMath.Tests/DoubleDoubleBoundaryTests.cs
+++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleBoundaryTests.cs
@@ -6,6 +6,60 @@ namespace Just.PreciseMath.Tests;
public class DoubleDoubleBoundaryTests
{
+ [Fact]
+ public void BoundarySubtractionPreservesComponentsOfTheOriginalExpression()
+ {
+ DoubleDouble[] magnitudes =
+ [
+ new(0.0), new(double.Epsilon), new(1.0),
+ new(Math.BitDecrement(Math.ScaleB(1.0, 1021))),
+ new(Math.ScaleB(1.0, 1021)),
+ DoubleDouble.FromComponents(Math.ScaleB(1.0, 1021), double.Epsilon),
+ DoubleDouble.FromComponents(Math.ScaleB(1.0, 1021), Math.ScaleB(1.0, 967)),
+ new(double.MaxValue),
+ DoubleDouble.FromComponents(double.MaxValue, Math.ScaleB(1.0, 969)),
+ DoubleDouble.FromComponents(double.MaxValue, -Math.ScaleB(1.0, 969))
+ ];
+ foreach (DoubleDouble leftMagnitude in magnitudes)
+ {
+ foreach (DoubleDouble rightMagnitude in magnitudes)
+ {
+ // Use the original exponent-domain rule, not the production predicate.
+ if (Math.ILogB(leftMagnitude.High) <= 1020 && Math.ILogB(rightMagnitude.High) <= 1020)
+ {
+ continue;
+ }
+ foreach (double leftSign in new[] { -1.0, 1.0 })
+ {
+ foreach (double rightSign in new[] { -1.0, 1.0 })
+ {
+ DoubleDouble left = leftSign < 0.0 ? -leftMagnitude : leftMagnitude;
+ DoubleDouble right = rightSign < 0.0 ? -rightMagnitude : rightMagnitude;
+ DoubleDouble actual = left - right;
+ // Bitwise characterization of the expression being extracted,
+ // supplemented below by the independent exact-rational oracle.
+ DoubleDouble original = PreciseMathHelper.AddBoundary(left, -right);
+ string context = Describe(left, right, "-");
+ BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(BitConverter.DoubleToInt64Bits(original.High), context);
+ BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(BitConverter.DoubleToInt64Bits(original.Low), context);
+ Rational expected = Exact(left) - Exact(right);
+ if (BelowOverflowMidpoint(expected))
+ {
+ AssertAccurate(actual, expected, context);
+ AssertNormalized(actual);
+ }
+ else
+ {
+ actual.High.ShouldBe(expected.CompareTo(Exact(0.0)) < 0
+ ? double.NegativeInfinity : double.PositiveInfinity, context);
+ BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L, context);
+ }
+ }
+ }
+ }
+ }
+ }
+
[Theory]
[InlineData(1.0, "+")]
[InlineData(-1.0, "+")]
diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleRepresentationTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleRepresentationTests.cs
index 945fb62..695dc78 100644
--- a/1-tests/Just.PreciseMath.Tests/DoubleDoubleRepresentationTests.cs
+++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleRepresentationTests.cs
@@ -89,6 +89,38 @@ public class DoubleDoubleRepresentationTests
(value != DoubleDouble.NaN).ShouldBeTrue();
}
+ [Fact]
+ public void FiniteNormalizationPreservesOrderingResidualsAndZeroSigns()
+ {
+ // Exact dyadic sums, independent of the public factory. Include both
+ // magnitude orders, cancellation, subnormals, and signed zero lows.
+ double unit = Math.ScaleB(1.0, 970);
+ (double High, double Low, double ExpectedHigh, double ExpectedLow)[] cases =
+ [
+ (-0.0, 0.0, -0.0, 0.0), (-0.0, -0.0, -0.0, 0.0),
+ (0.0, -0.0, 0.0, 0.0), (1.0, -1.0, 0.0, 0.0),
+ (-1.0, 1.0, 0.0, 0.0), (double.Epsilon, -double.Epsilon, 0.0, 0.0),
+ (double.Epsilon, double.Epsilon, 2 * double.Epsilon, 0.0),
+ (0.0, double.Epsilon, double.Epsilon, 0.0),
+ (double.Epsilon, -1.0, -1.0, double.Epsilon),
+ (1.0, double.Epsilon, 1.0, double.Epsilon),
+ (double.Epsilon, 1.0, 1.0, double.Epsilon),
+ (-1.0, -double.Epsilon, -1.0, -double.Epsilon),
+ (-double.Epsilon, -1.0, -1.0, -double.Epsilon),
+ (double.MaxValue, 0.0, double.MaxValue, 0.0),
+ // MaxValue - 3*2^970 = (MaxValue - 2^971) - 2^970.
+ // An unordered TwoSum would overflow an intermediate in the reversed case.
+ (double.MaxValue, -3.0 * unit, Math.BitDecrement(double.MaxValue), -unit),
+ (-3.0 * unit, double.MaxValue, Math.BitDecrement(double.MaxValue), -unit)
+ ];
+ foreach ((double high, double low, double expectedHigh, double expectedLow) in cases)
+ {
+ DoubleDouble actual = PreciseMathHelper.NormalizeFinite(high, low);
+ BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(BitConverter.DoubleToInt64Bits(expectedHigh));
+ BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(BitConverter.DoubleToInt64Bits(expectedLow));
+ }
+ }
+
[Fact]
public void ZeroSignsArePreservedButEqual()
{
diff --git a/2-benchmarks/Just.PreciseMath.Benchmarks/AddSubtractMatchedBenchmarks.cs b/2-benchmarks/Just.PreciseMath.Benchmarks/AddSubtractMatchedBenchmarks.cs
new file mode 100644
index 0000000..79db244
--- /dev/null
+++ b/2-benchmarks/Just.PreciseMath.Benchmarks/AddSubtractMatchedBenchmarks.cs
@@ -0,0 +1,134 @@
+using BenchmarkDotNet.Attributes;
+
+namespace Just.PreciseMath.Benchmarks;
+
+/// Compares a+b with a-(-b) using matched results and normalization paths.
+/// Includes array and loop costs. Negation and fixture validation are outside timing.
+[MemoryDiagnoser]
+public class AddSubtractMatchedBenchmarks
+{
+ private const int Count = 64;
+ private readonly DoubleDouble[] _left = new DoubleDouble[Count];
+ private readonly DoubleDouble[] _right = new DoubleDouble[Count];
+ private readonly DoubleDouble[] _results = new DoubleDouble[Count];
+
+ [Params("ZeroCorrection", "NonzeroCorrection", "Cancellation", "MixedCorrection")]
+ public string Scenario { get; set; } = "NonzeroCorrection";
+
+ [GlobalSetup(Target = nameof(Add))]
+ public void SetupAdd()
+ {
+ Setup(false);
+ }
+
+ [GlobalSetup(Target = nameof(Subtract))]
+ public void SetupSubtract()
+ {
+ Setup(true);
+ }
+
+ [Benchmark(Baseline = true, OperationsPerInvoke = Count)]
+ public DoubleDouble[] Add()
+ {
+ for (int i = 0; i < Count; i++)
+ {
+ _results[i] = _left[i] + _right[i];
+ }
+ return _results;
+ }
+
+ [Benchmark(OperationsPerInvoke = Count)]
+ public DoubleDouble[] Subtract()
+ {
+ for (int i = 0; i < Count; i++)
+ {
+ _results[i] = _left[i] - _right[i];
+ }
+ return _results;
+ }
+
+ private void Setup(bool subtract)
+ {
+ if (Scenario is not ("ZeroCorrection" or "NonzeroCorrection" or "Cancellation" or "MixedCorrection"))
+ {
+ throw new InvalidOperationException($"Unknown matched add/subtract scenario: {Scenario}.");
+ }
+ for (int i = 0; i < Count; i++)
+ {
+ // Cycle both high signs and moderate scales. High sums/differences
+ // are exact: the significands are small integer multiples of 1/128.
+ int exponent = (((i / 4) % 3) - 1) * 80;
+ double leftSign = i % 2 == 0 ? 1.0 : -1.0;
+ double rightSign = i % 4 < 2 ? 1.0 : -1.0;
+ double leftHigh = leftSign * Math.ScaleB(1.25 + ((i % 7) / 32.0), exponent);
+ double rightHigh = rightSign * Math.ScaleB(1.0 + (((i % 5) - 2) / 128.0), exponent);
+ bool zeroCorrection = Scenario == "ZeroCorrection" || (Scenario == "MixedCorrection" && i % 3 == 0);
+ double leftLow;
+ double rightLow;
+ double expectedHigh;
+ double expectedLow;
+ if (Scenario == "Cancellation")
+ {
+ // Equal/opposite highs leave the exact normalized expansion
+ // (sign*2^(e-54), sign*2^(e-108)), not a binary64-only result.
+ rightHigh = -leftHigh;
+ leftLow = leftSign * Math.ScaleB(1.0, exponent - 54);
+ rightLow = leftSign * Math.ScaleB(1.0, exponent - 108);
+ expectedHigh = leftLow;
+ expectedLow = rightLow;
+ }
+ else
+ {
+ leftLow = leftSign * Math.ScaleB(1.0, exponent - 80);
+ rightLow = zeroCorrection ? -leftLow : leftLow / 2.0;
+ expectedHigh = leftHigh + rightHigh;
+ // This exact low sum is far below half an ulp of the high sum.
+ expectedLow = zeroCorrection ? 0.0 : leftSign * Math.ScaleB(3.0, exponent - 81);
+ }
+
+ DoubleDouble left = DoubleDouble.FromComponents(leftHigh, leftLow);
+ DoubleDouble right = DoubleDouble.FromComponents(rightHigh, rightLow);
+ DoubleDouble negativeRight = -right;
+ CheckComponents(left, leftHigh, leftLow, i, "left input");
+ CheckComponents(right, rightHigh, rightLow, i, "right input");
+ CheckComponents(negativeRight, -rightHigh, -rightLow, i, "negated right input");
+ if (!DoubleDouble.IsFinite(left) || !DoubleDouble.IsFinite(right)
+ || left.High == 0.0 || right.High == 0.0
+ || Math.Abs(Math.ILogB(left.High)) > 200 || Math.Abs(Math.ILogB(right.High)) > 200)
+ {
+ throw new InvalidOperationException($"{Scenario}[{i}] is outside the ordinary finite input range.");
+ }
+ // Independent dyadic expectations, not one DD operator as the other's oracle.
+ CheckComponents(left + right, expectedHigh, expectedLow, i, "addition");
+ CheckComponents(left - negativeRight, expectedHigh, expectedLow, i, "subtraction");
+ // Both operators feed these same effective components into AddFinite.
+ if ((FinalCorrection(left, right) == 0.0) != zeroCorrection)
+ {
+ throw new InvalidOperationException($"{Scenario}[{i}] does not select the expected normalization path.");
+ }
+ _left[i] = left;
+ // Target-specific setup keeps array fields and timed memory access identical.
+ _right[i] = subtract ? negativeRight : right;
+ }
+ }
+
+ private void CheckComponents(DoubleDouble actual, double high, double low, int index, string operation)
+ {
+ if (BitConverter.DoubleToInt64Bits(actual.High) != BitConverter.DoubleToInt64Bits(high)
+ || BitConverter.DoubleToInt64Bits(actual.Low) != BitConverter.DoubleToInt64Bits(low))
+ {
+ throw new InvalidOperationException($"{Scenario}[{index}] has unexpected {operation} components.");
+ }
+ }
+
+ // Setup-only characterization, not an accuracy oracle. Keep this path check
+ // synchronized with AddFinite if the production kernel changes in the future.
+ private static double FinalCorrection(DoubleDouble left, DoubleDouble right)
+ {
+ (double high, double highError) = PreciseMathHelper.TwoAdd(left.High, right.High);
+ (double low, double lowError) = PreciseMathHelper.TwoAdd(left.Low, right.Low);
+ (double middle, double middleError) = PreciseMathHelper.TwoAdd(highError, low);
+ (double _, double sumError) = PreciseMathHelper.TwoAdd(high, middle);
+ return sumError + (middleError + lowError);
+ }
+}
diff --git a/2-benchmarks/Just.PreciseMath.Benchmarks/ArithmeticBoundaryBenchmarks.cs b/2-benchmarks/Just.PreciseMath.Benchmarks/ArithmeticBoundaryBenchmarks.cs
new file mode 100644
index 0000000..17f3b3c
--- /dev/null
+++ b/2-benchmarks/Just.PreciseMath.Benchmarks/ArithmeticBoundaryBenchmarks.cs
@@ -0,0 +1,114 @@
+using BenchmarkDotNet.Attributes;
+
+namespace Just.PreciseMath.Benchmarks;
+
+/// Isolates ordinary/fallback exponent transitions, not special-value handling.
+[MemoryDiagnoser]
+[CategoriesColumn]
+public class ArithmeticBoundaryBenchmarks
+{
+ private DoubleDouble _addLeft;
+ private DoubleDouble _addRight;
+ private DoubleDouble _multiplyLeft;
+ private DoubleDouble _divideLeft;
+ private DoubleDouble _factor;
+
+ [Params(false, true)]
+ public bool Fallback { get; set; }
+
+ [GlobalSetup]
+ public void Setup()
+ {
+ int offset = Fallback ? 1 : 0;
+ _addLeft = Pair(1020 + offset);
+ _addRight = Pair(1019);
+ _multiplyLeft = Pair(900 + offset);
+ _divideLeft = Pair(450 + offset);
+ _factor = DoubleDouble.FromComponents(1.125, Math.ScaleB(1.0, -56));
+ // Pin dispatch assumptions to the current source guards. Both signs/orders
+ // remain finite; boundary timing is intentionally separate from ordinary data.
+ if ((Math.ILogB(_addLeft.High) > 1020) != Fallback
+ || (Math.ILogB(_multiplyLeft.High) + Math.ILogB(_factor.High) > 900) != Fallback
+ || (Math.Abs(Math.ILogB(_divideLeft.High)) > 450) != Fallback)
+ {
+ throw new InvalidOperationException("Boundary fixture does not select the requested guard branch.");
+ }
+ }
+
+ private static DoubleDouble Pair(int exponent)
+ {
+ return DoubleDouble.FromComponents(Math.ScaleB(1.0, exponent), Math.ScaleB(1.0, exponent - 56));
+ }
+
+ [Benchmark, BenchmarkCategory("Addition")]
+ public DoubleDouble DDAdd()
+ {
+ return _addLeft + _addRight;
+ }
+
+ [Benchmark, BenchmarkCategory("Addition")]
+ public DoubleDouble DDScalarAdd()
+ {
+ return _addLeft + _addRight.High;
+ }
+
+ [Benchmark, BenchmarkCategory("Addition")]
+ public DoubleDouble ScalarDDAdd()
+ {
+ return _addLeft.High + _addRight;
+ }
+
+ [Benchmark, BenchmarkCategory("Subtraction")]
+ public DoubleDouble DDSubtract()
+ {
+ return _addLeft - _addRight;
+ }
+
+ [Benchmark, BenchmarkCategory("Subtraction")]
+ public DoubleDouble DDScalarSubtract()
+ {
+ return _addLeft - _addRight.High;
+ }
+
+ [Benchmark, BenchmarkCategory("Subtraction")]
+ public DoubleDouble ScalarDDSubtract()
+ {
+ return _addLeft.High - _addRight;
+ }
+
+ [Benchmark, BenchmarkCategory("Multiplication")]
+ public DoubleDouble DDMultiply()
+ {
+ return _multiplyLeft * _factor;
+ }
+
+ [Benchmark, BenchmarkCategory("Multiplication")]
+ public DoubleDouble DDScalarMultiply()
+ {
+ return _multiplyLeft * _factor.High;
+ }
+
+ [Benchmark, BenchmarkCategory("Multiplication")]
+ public DoubleDouble ScalarDDMultiply()
+ {
+ return _multiplyLeft.High * _factor;
+ }
+
+ [Benchmark, BenchmarkCategory("Division")]
+ public DoubleDouble DDDivide()
+ {
+ return _divideLeft / _factor;
+ }
+
+ [Benchmark, BenchmarkCategory("Division")]
+ public DoubleDouble DDScalarDivide()
+ {
+ return _divideLeft / _factor.High;
+ }
+
+ [Benchmark, BenchmarkCategory("Division")]
+ public DoubleDouble ScalarDDDivide()
+ {
+ return _factor.High / _divideLeft;
+ }
+}
diff --git a/2-benchmarks/Just.PreciseMath.Benchmarks/ArithmeticInputs.cs b/2-benchmarks/Just.PreciseMath.Benchmarks/ArithmeticInputs.cs
new file mode 100644
index 0000000..3a86b83
--- /dev/null
+++ b/2-benchmarks/Just.PreciseMath.Benchmarks/ArithmeticInputs.cs
@@ -0,0 +1,129 @@
+namespace Just.PreciseMath.Benchmarks;
+
+// Input construction stays outside timed benchmark methods.
+internal sealed class ArithmeticInputs
+{
+ internal const int Count = 64;
+
+ internal DoubleDouble[] Left { get; } = new DoubleDouble[Count];
+ internal DoubleDouble[] Right { get; } = new DoubleDouble[Count];
+
+ internal static ArithmeticInputs Create(string scenario, bool chain)
+ {
+ if (chain && scenario is not ("BothResidual" or "Mixed"))
+ {
+ throw new ArgumentOutOfRangeException(nameof(scenario));
+ }
+ if (scenario is not ("BinaryExact" or "DecimalResidual" or "BothResidual" or "Mixed" or "SignsAndScales" or "Cancellation"))
+ {
+ throw new ArgumentOutOfRangeException(nameof(scenario));
+ }
+
+ ArithmeticInputs inputs = new();
+ for (int i = 0; i < Count; i++)
+ {
+ // Exact dyadic construction, not decimal-to-double residual estimation.
+ double high = 1.25 + ((i % 7) / 32.0);
+ double low = Math.ScaleB(1.0, -56);
+ DoubleDouble left = DoubleDouble.FromComponents(high, low);
+ DoubleDouble right = DoubleDouble.FromComponents(1.0 + (((i % 5) - 2) / 128.0), -low);
+ if (!chain)
+ {
+ switch (scenario)
+ {
+ case "BinaryExact":
+ left = new DoubleDouble(high);
+ right = new DoubleDouble(0.75);
+ break;
+ case "DecimalResidual":
+ left = (DoubleDouble)1.1m;
+ right = new DoubleDouble(0.75);
+ break;
+ case "SignsAndScales":
+ int exponent = ((i % 7) - 3) * 40;
+ double sign = i % 2 == 0 ? 1.0 : -1.0;
+ left = DoubleDouble.FromComponents(Math.ScaleB(sign * high, exponent), Math.ScaleB(low, exponent));
+ right = DoubleDouble.FromComponents(Math.ScaleB(1.125, -exponent), Math.ScaleB(-low, -exponent));
+ break;
+ case "Cancellation":
+ // Alternate near cancellation in addition and subtraction.
+ right = DoubleDouble.FromComponents(i % 2 == 0 ? -high : high, -low / 2.0);
+ break;
+ }
+ }
+ if (scenario == "Mixed")
+ {
+ // Alternate both residuals, right zero-low, left zero-low, both zero-low.
+ if (i % 4 is 2 or 3)
+ {
+ left = new DoubleDouble(left.High);
+ }
+ if (i % 4 is 1 or 3)
+ {
+ right = new DoubleDouble(right.High);
+ }
+ }
+ inputs.Left[i] = left;
+ inputs.Right[i] = right;
+ }
+ inputs.Validate(scenario, chain);
+ return inputs;
+ }
+
+ private void Validate(string scenario, bool chain)
+ {
+ for (int i = 0; i < Count; i++)
+ {
+ ValidateOrdinary(Left[i]);
+ ValidateOrdinary(Right[i]);
+ bool expectedLeftLow = scenario != "BinaryExact" && (scenario != "Mixed" || i % 4 < 2);
+ bool expectedRightLow = scenario is not ("BinaryExact" or "DecimalResidual") && (scenario != "Mixed" || i % 2 == 0);
+ if ((Left[i].Low != 0.0) != expectedLeftLow || (Right[i].Low != 0.0) != expectedRightLow)
+ {
+ throw new InvalidOperationException("Benchmark residual shape does not match its scenario.");
+ }
+ }
+ if (!chain)
+ {
+ return;
+ }
+
+ // Setup-only guard against a chain drifting into special or exponent-boundary
+ // paths. This is fixture validation, not an independent accuracy oracle.
+ for (int operation = 0; operation < 12; operation++)
+ {
+ DoubleDouble value = Left[0];
+ for (int i = 0; i < Count; i++)
+ {
+ DoubleDouble right = Right[i];
+ value = operation switch
+ {
+ 0 => value + right,
+ 1 => value - right,
+ 2 => value * right,
+ 3 => value / right,
+ 4 => value + right.High,
+ 5 => value - right.High,
+ 6 => value * right.High,
+ 7 => value / right.High,
+ 8 => right.High + value,
+ 9 => right.High - value,
+ 10 => right.High * value,
+ 11 => right.High / value,
+ _ => throw new InvalidOperationException()
+ };
+ ValidateOrdinary(value);
+ }
+ }
+ }
+
+ private static void ValidateOrdinary(DoubleDouble value)
+ {
+ if (!DoubleDouble.IsFinite(value) || value.High == 0.0 || Math.Abs(Math.ILogB(value.High)) > 200
+ || DoubleDouble.FromComponents(value.High, value.Low) != value
+ || (value.Low == 0.0 && BitConverter.DoubleToInt64Bits(value.Low) != 0))
+ {
+ throw new InvalidOperationException("Benchmark input/chain left the normalized ordinary finite range.");
+ }
+ }
+}
diff --git a/2-benchmarks/Just.PreciseMath.Benchmarks/ArithmeticLatencyBenchmarks.cs b/2-benchmarks/Just.PreciseMath.Benchmarks/ArithmeticLatencyBenchmarks.cs
new file mode 100644
index 0000000..beb74e9
--- /dev/null
+++ b/2-benchmarks/Just.PreciseMath.Benchmarks/ArithmeticLatencyBenchmarks.cs
@@ -0,0 +1,157 @@
+using BenchmarkDotNet.Attributes;
+
+namespace Just.PreciseMath.Benchmarks;
+
+/// Measures bounded dependent chains; each invocation resets the seed.
+[MemoryDiagnoser]
+[CategoriesColumn]
+public class ArithmeticLatencyBenchmarks
+{
+ private DoubleDouble[] _left = [];
+ private DoubleDouble[] _right = [];
+ private double[] _scalarRight = [];
+
+ [Params("BothResidual", "Mixed")]
+ public string Scenario { get; set; } = "BothResidual";
+
+ [GlobalSetup]
+ public void Setup()
+ {
+ ArithmeticInputs inputs = ArithmeticInputs.Create(Scenario, true);
+ _left = inputs.Left;
+ _right = inputs.Right;
+ _scalarRight = Array.ConvertAll(_right, value => value.High);
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Addition")]
+ public DoubleDouble DDAdd()
+ {
+ DoubleDouble value = _left[0];
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ value = value + _right[i];
+ }
+ return value;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Addition")]
+ public DoubleDouble DDScalarAdd()
+ {
+ DoubleDouble value = _left[0];
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ value = value + _scalarRight[i];
+ }
+ return value;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Addition")]
+ public DoubleDouble ScalarDDAdd()
+ {
+ DoubleDouble value = _left[0];
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ value = _scalarRight[i] + value;
+ }
+ return value;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Subtraction")]
+ public DoubleDouble DDSubtract()
+ {
+ DoubleDouble value = _left[0];
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ value = value - _right[i];
+ }
+ return value;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Subtraction")]
+ public DoubleDouble DDScalarSubtract()
+ {
+ DoubleDouble value = _left[0];
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ value = value - _scalarRight[i];
+ }
+ return value;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Subtraction")]
+ public DoubleDouble ScalarDDSubtract()
+ {
+ DoubleDouble value = _left[0];
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ value = _scalarRight[i] - value;
+ }
+ return value;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Multiplication")]
+ public DoubleDouble DDMultiply()
+ {
+ DoubleDouble value = _left[0];
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ value = value * _right[i];
+ }
+ return value;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Multiplication")]
+ public DoubleDouble DDScalarMultiply()
+ {
+ DoubleDouble value = _left[0];
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ value = value * _scalarRight[i];
+ }
+ return value;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Multiplication")]
+ public DoubleDouble ScalarDDMultiply()
+ {
+ DoubleDouble value = _left[0];
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ value = _scalarRight[i] * value;
+ }
+ return value;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Division")]
+ public DoubleDouble DDDivide()
+ {
+ DoubleDouble value = _left[0];
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ value = value / _right[i];
+ }
+ return value;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Division")]
+ public DoubleDouble DDScalarDivide()
+ {
+ DoubleDouble value = _left[0];
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ value = value / _scalarRight[i];
+ }
+ return value;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Division")]
+ public DoubleDouble ScalarDDDivide()
+ {
+ DoubleDouble value = _left[0];
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ value = _scalarRight[i] / value;
+ }
+ return value;
+ }
+}
diff --git a/2-benchmarks/Just.PreciseMath.Benchmarks/ArithmeticThroughputBenchmarks.cs b/2-benchmarks/Just.PreciseMath.Benchmarks/ArithmeticThroughputBenchmarks.cs
new file mode 100644
index 0000000..f4ad681
--- /dev/null
+++ b/2-benchmarks/Just.PreciseMath.Benchmarks/ArithmeticThroughputBenchmarks.cs
@@ -0,0 +1,189 @@
+using BenchmarkDotNet.Attributes;
+
+namespace Just.PreciseMath.Benchmarks;
+
+/// Measures independent array operations, including loads, stores, and loop overhead.
+[MemoryDiagnoser]
+[CategoriesColumn]
+public class ArithmeticThroughputBenchmarks
+{
+ private DoubleDouble[] _left = [];
+ private DoubleDouble[] _right = [];
+ private double[] _scalarLeft = [];
+ private double[] _scalarRight = [];
+ private readonly DoubleDouble[] _results = new DoubleDouble[ArithmeticInputs.Count];
+ private readonly double[] _doubleResults = new double[ArithmeticInputs.Count];
+
+ [Params("BinaryExact", "DecimalResidual", "BothResidual", "Mixed", "SignsAndScales", "Cancellation")]
+ public string Scenario { get; set; } = "BothResidual";
+
+ [GlobalSetup]
+ public void Setup()
+ {
+ ArithmeticInputs inputs = ArithmeticInputs.Create(Scenario, false);
+ _left = inputs.Left;
+ _right = inputs.Right;
+ _scalarLeft = Array.ConvertAll(_left, value => value.High);
+ _scalarRight = Array.ConvertAll(_right, value => value.High);
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Addition")]
+ public DoubleDouble[] DDAdd()
+ {
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ _results[i] = _left[i] + _right[i];
+ }
+ return _results;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Addition")]
+ public DoubleDouble[] DDScalarAdd()
+ {
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ _results[i] = _left[i] + _scalarRight[i];
+ }
+ return _results;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Addition")]
+ public DoubleDouble[] ScalarDDAdd()
+ {
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ _results[i] = _scalarLeft[i] + _right[i];
+ }
+ return _results;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Addition")]
+ public double[] DoubleAdd()
+ {
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ _doubleResults[i] = _scalarLeft[i] + _scalarRight[i];
+ }
+ return _doubleResults;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Subtraction")]
+ public DoubleDouble[] DDSubtract()
+ {
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ _results[i] = _left[i] - _right[i];
+ }
+ return _results;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Subtraction")]
+ public DoubleDouble[] DDScalarSubtract()
+ {
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ _results[i] = _left[i] - _scalarRight[i];
+ }
+ return _results;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Subtraction")]
+ public DoubleDouble[] ScalarDDSubtract()
+ {
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ _results[i] = _scalarLeft[i] - _right[i];
+ }
+ return _results;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Subtraction")]
+ public double[] DoubleSubtract()
+ {
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ _doubleResults[i] = _scalarLeft[i] - _scalarRight[i];
+ }
+ return _doubleResults;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Multiplication")]
+ public DoubleDouble[] DDMultiply()
+ {
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ _results[i] = _left[i] * _right[i];
+ }
+ return _results;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Multiplication")]
+ public DoubleDouble[] DDScalarMultiply()
+ {
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ _results[i] = _left[i] * _scalarRight[i];
+ }
+ return _results;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Multiplication")]
+ public DoubleDouble[] ScalarDDMultiply()
+ {
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ _results[i] = _scalarLeft[i] * _right[i];
+ }
+ return _results;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Multiplication")]
+ public double[] DoubleMultiply()
+ {
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ _doubleResults[i] = _scalarLeft[i] * _scalarRight[i];
+ }
+ return _doubleResults;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Division")]
+ public DoubleDouble[] DDDivide()
+ {
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ _results[i] = _left[i] / _right[i];
+ }
+ return _results;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Division")]
+ public DoubleDouble[] DDScalarDivide()
+ {
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ _results[i] = _left[i] / _scalarRight[i];
+ }
+ return _results;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Division")]
+ public DoubleDouble[] ScalarDDDivide()
+ {
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ _results[i] = _scalarLeft[i] / _right[i];
+ }
+ return _results;
+ }
+
+ [Benchmark(OperationsPerInvoke = ArithmeticInputs.Count), BenchmarkCategory("Division")]
+ public double[] DoubleDivide()
+ {
+ for (int i = 0; i < ArithmeticInputs.Count; i++)
+ {
+ _doubleResults[i] = _scalarLeft[i] / _scalarRight[i];
+ }
+ return _doubleResults;
+ }
+}
diff --git a/README.md b/README.md
index 4bc5576..05f9038 100644
--- a/README.md
+++ b/README.md
@@ -133,7 +133,7 @@ bool success = DoubleDouble.TryParse("1.25e-2".AsSpan(), CultureInfo.InvariantCu
The planned `PreciseMath` static class and its `Abs`, `Sqrt`, `Pow`, `Exp`, and
`Log` functions are not implemented. Generic-math interfaces beyond `ISignedNumber`,
-additional text formats/general round-trip formatting, and broader performance
+additional text formats/general round-trip formatting, and non-arithmetic performance
benchmarks remain deferred.
Replacing allocating arithmetic boundary fallbacks is also deferred; the current
`BigInteger` paths remain in place. That optimization does not require removing
@@ -156,18 +156,9 @@ on CI workflow runs.
## Benchmarks
-The BenchmarkDotNet suite compares same-type `+`, `-`, `*`, and `/` operations for
-`DoubleDouble`, `decimal`, and `double`: 12 methods with two input cases each.
-Each operation/input group uses `double` as its baseline and reports allocations.
-Operands are stored in fields and results are returned to the harness; conversions
-and construction happen in setup, outside the timed methods.
-
-The left inputs are `1.25` (binary-exact) and `1.1` (a nonzero low component in
-`DoubleDouble`); the right input is `0.75`. Decimal and double-double inputs originate
-from the same decimal values, while `double` rounds to binary64. These types have
-different precision and range contracts: this is a cost comparison, not an accuracy
-test. Exceptional values, exponent-boundary fallbacks, mixed-type operators,
-conversions, parsing, and formatting are not benchmarked yet.
+BenchmarkDotNet measures arithmetic throughput, dependent-chain latency, and allocations,
+including comparisons of `DoubleDouble`, `decimal`, and `double`, mixed scalar operations,
+and exponent-boundary paths. These are performance measurements, not accuracy tests.
After the Release build above, run from the repository root:
@@ -175,19 +166,16 @@ After the Release build above, run from the repository root:
# Discover benchmark methods without running them.
dotnet run --project 2-benchmarks/Just.PreciseMath.Benchmarks -c Release --no-build -- --list flat
-# Quick execution check (also run in CI); not useful for timing comparisons.
+# Smoke test; Dry timings are not performance measurements.
dotnet run --project 2-benchmarks/Just.PreciseMath.Benchmarks -c Release --no-build -- --job Dry --filter '*'
-# Full measurement run; optionally select an operation with --anyCategories Addition.
+# Full measurement run.
dotnet run --project 2-benchmarks/Just.PreciseMath.Benchmarks -c Release --no-build -- --filter '*'
```
Reports are written under the ignored `BenchmarkDotNet.Artifacts/` directory.
-Empty selections and failed benchmark runs return a nonzero exit code. CI only
-checks execution, with no performance gate. Use full runs on an idle, controlled
-machine for comparisons; a scalar `double` operation may approach harness overhead,
-so inspect BenchmarkDotNet warnings before interpreting ratios. Do not compare
-coverage-instrumented runs or treat Dry-job timings as performance measurements.
+Replace `'*'` with a benchmark-name pattern to select a subset; use `--artifacts `
+to keep runs separate. Run measurements on an idle machine and inspect BenchmarkDotNet warnings.
## Project structure