This commit is contained in:
@@ -182,6 +182,170 @@ public class DoubleDoubleArithmeticTests
|
||||
Check(DoubleDouble.FromComponents(sign * double.MaxValue, halfUlp), sign * double.PositiveInfinity, 0.0);
|
||||
}
|
||||
|
||||
[Theory]
|
||||
[InlineData(-1)]
|
||||
[InlineData(0)]
|
||||
[InlineData(1)]
|
||||
public void ProductAndQuotientCrossTheExactOverflowMidpoint(int side)
|
||||
{
|
||||
// M = 2^1024 - 2^970. M/2 is the normalized pair (2^1023, -2^969).
|
||||
// Its adjacent normalized pairs use different highs across this tie:
|
||||
// below uses (MaxValue/2, BitDecrement(2^969)), above increments -2^969.
|
||||
// Doubling gives M +/- 2^917; below has residual BitDecrement(2^970).
|
||||
// The tie rounds to the even significand at 2^1024, hence infinity.
|
||||
double high = side < 0 ? Math.ScaleB(double.MaxValue, -1) : Math.ScaleB(1.0, 1023);
|
||||
double low = side switch
|
||||
{
|
||||
-1 => Math.BitDecrement(Math.ScaleB(1.0, 969)),
|
||||
1 => Math.BitIncrement(-Math.ScaleB(1.0, 969)),
|
||||
_ => -Math.ScaleB(1.0, 969)
|
||||
};
|
||||
BigInteger midpointUnits = Units(double.MaxValue) + Units(Math.ScaleB(1.0, 970));
|
||||
foreach (double sign in new[] { -1.0, 1.0 })
|
||||
{
|
||||
DoubleDouble value = DoubleDouble.FromComponents(sign * high, sign * low);
|
||||
(BigInteger.Abs(Units(value)) * 2).ShouldBe(midpointUnits + (side * Units(Math.ScaleB(1.0, 917))));
|
||||
double expectedHigh = sign * (side < 0 ? double.MaxValue : double.PositiveInfinity);
|
||||
double expectedLow = side < 0 ? sign * Math.BitDecrement(Math.ScaleB(1.0, 970)) : 0.0;
|
||||
CheckBoundary(value * new DoubleDouble(2.0), expectedHigh, expectedLow);
|
||||
CheckBoundary(new DoubleDouble(2.0) * value, expectedHigh, expectedLow);
|
||||
CheckBoundary(value * 2.0, expectedHigh, expectedLow);
|
||||
CheckBoundary(2.0 * value, expectedHigh, expectedLow);
|
||||
CheckBoundary(value / new DoubleDouble(0.5), expectedHigh, expectedLow);
|
||||
CheckBoundary(value / 0.5, expectedHigh, expectedLow);
|
||||
}
|
||||
}
|
||||
|
||||
[Theory]
|
||||
[InlineData(false)]
|
||||
[InlineData(true)]
|
||||
public void ScalarLeftDivisionStraddlesTheOverflowMidpoint(bool below)
|
||||
{
|
||||
// An exact M = (2^54 - 1)*2^970 quotient is impossible with a finite
|
||||
// binary64 numerator and dyadic denominator: its odd numerator would
|
||||
// need all 54 bits. Instead use adjacent low components bracketing
|
||||
// 2^1023/M = 1/2 + 2^-55 + 2^-109 + ... .
|
||||
double low = Math.ScaleB(1.0, -55);
|
||||
DoubleDouble denominator = DoubleDouble.FromComponents(0.5, below ? Math.BitIncrement(low) : low);
|
||||
BigInteger midpointUnits = Units(double.MaxValue) + Units(Math.ScaleB(1.0, 970));
|
||||
BigInteger numeratorUnits = Units(Math.ScaleB(1.0, 1023));
|
||||
((numeratorUnits << 1074) < midpointUnits * Units(denominator)).ShouldBe(below);
|
||||
// Exact rational residual rounding gives BitDecrement(2^970) below M.
|
||||
foreach (double sign in new[] { -1.0, 1.0 })
|
||||
{
|
||||
double numerator = sign * Math.ScaleB(1.0, 1023);
|
||||
double high = sign * (below ? double.MaxValue : double.PositiveInfinity);
|
||||
double residual = below ? sign * Math.BitDecrement(Math.ScaleB(1.0, 970)) : 0.0;
|
||||
CheckBoundary(numerator / denominator, high, residual);
|
||||
CheckBoundary(new DoubleDouble(numerator) / denominator, high, residual);
|
||||
CheckBoundary((-numerator) / (-denominator), high, residual);
|
||||
}
|
||||
}
|
||||
|
||||
[Theory]
|
||||
[InlineData(-1)]
|
||||
[InlineData(0)]
|
||||
[InlineData(1)]
|
||||
public void ProductAndQuotientRoundUnderflowTiesToSignedZero(int side)
|
||||
{
|
||||
// (2^-1022 + side*2^-1074)*2^-53 = epsilon/2 + side*2^-1127.
|
||||
// The low component is stepped by its smallest possible increment.
|
||||
// Ties select even zero, retaining the exact nonzero result's sign.
|
||||
foreach (double sign in new[] { -1.0, 1.0 })
|
||||
{
|
||||
DoubleDouble value = DoubleDouble.FromComponents(sign * Math.ScaleB(1.0, -1022), sign * side * double.Epsilon);
|
||||
double multiplier = Math.ScaleB(1.0, -53);
|
||||
double divisor = Math.ScaleB(1.0, 53);
|
||||
double expected = Math.CopySign(side > 0 ? double.Epsilon : 0.0, sign);
|
||||
CheckBits(value * new DoubleDouble(multiplier), expected);
|
||||
CheckBits(new DoubleDouble(multiplier) * value, expected);
|
||||
CheckBits(value * multiplier, expected);
|
||||
CheckBits(multiplier * value, expected);
|
||||
CheckBits(value / new DoubleDouble(divisor), expected);
|
||||
CheckBits(value / divisor, expected);
|
||||
|
||||
// epsilon/(2 - side*epsilon) brackets the same tie; the nonzero
|
||||
// denominator residual is essential despite being invisible in double.
|
||||
DoubleDouble denominator = DoubleDouble.FromComponents(2.0, -side * double.Epsilon);
|
||||
CheckBits((sign * double.Epsilon) / denominator, expected);
|
||||
CheckBits(new DoubleDouble(sign * double.Epsilon) / denominator, expected);
|
||||
CheckBits((-sign * double.Epsilon) / (-denominator), expected);
|
||||
}
|
||||
}
|
||||
|
||||
[Theory]
|
||||
[InlineData(0)]
|
||||
[InlineData(1)]
|
||||
[InlineData(2)]
|
||||
[InlineData(3)]
|
||||
[InlineData(4)]
|
||||
[InlineData(5)]
|
||||
[InlineData(6)]
|
||||
[InlineData(7)]
|
||||
[InlineData(8)]
|
||||
[InlineData(9)]
|
||||
[InlineData(10)]
|
||||
[InlineData(11)]
|
||||
public void WarmedFiniteArithmeticAllocatesSubstantiallyLessThanBoundaryFallback(int operation)
|
||||
{
|
||||
DoubleDouble ordinary = DoubleDouble.FromComponents(1.25, Math.ScaleB(1.0, -70));
|
||||
DoubleDouble boundary = DoubleDouble.FromComponents(Math.ScaleB(1.0, 1022), Math.ScaleB(1.0, 968));
|
||||
DoubleDouble right = DoubleDouble.FromComponents(1.5, Math.ScaleB(1.0, -55));
|
||||
// Synchronous per-thread counters exclude other parallel tests. Warm both
|
||||
// branches, keep setup/assertions outside measurement, and consume results.
|
||||
_ = MeasureArithmeticAllocations(operation, ordinary, right, 128, out _);
|
||||
_ = MeasureArithmeticAllocations(operation, boundary, right, 128, out _);
|
||||
long ordinaryBytes = long.MaxValue;
|
||||
long boundaryBytes = long.MaxValue;
|
||||
for (int sample = 0; sample < 3; ++sample)
|
||||
{
|
||||
long finite = MeasureArithmeticAllocations(operation, ordinary, right, 256, out double finiteChecksum);
|
||||
long fallback = MeasureArithmeticAllocations(operation, boundary, right, 256, out double fallbackChecksum);
|
||||
double.IsFinite(finiteChecksum).ShouldBeTrue();
|
||||
double.IsFinite(fallbackChecksum).ShouldBeTrue();
|
||||
ordinaryBytes = Math.Min(ordinaryBytes, finite);
|
||||
boundaryBytes = Math.Min(boundaryBytes, fallback);
|
||||
}
|
||||
// A coarse relative distinction, not a runtime-dependent BigInteger byte
|
||||
// count or timing benchmark. Minima discard incidental warm-up allocation.
|
||||
boundaryBytes.ShouldBeGreaterThan(0L);
|
||||
ordinaryBytes.ShouldBeLessThan(boundaryBytes / 16);
|
||||
}
|
||||
|
||||
private static long MeasureArithmeticAllocations(int operation, DoubleDouble left, DoubleDouble right, int iterations, out double checksum)
|
||||
{
|
||||
checksum = 0.0;
|
||||
long before = GC.GetAllocatedBytesForCurrentThread();
|
||||
for (int i = 0; i < iterations; ++i)
|
||||
{
|
||||
DoubleDouble result = operation switch
|
||||
{
|
||||
0 => left + right,
|
||||
1 => left - right,
|
||||
2 => left * right,
|
||||
3 => left / right,
|
||||
4 => left + right.High,
|
||||
5 => right.High + left,
|
||||
6 => left - right.High,
|
||||
7 => right.High - left,
|
||||
8 => left * right.High,
|
||||
9 => right.High * left,
|
||||
10 => left / right.High,
|
||||
11 => right.High / left,
|
||||
_ => throw new ArgumentOutOfRangeException(nameof(operation))
|
||||
};
|
||||
// Scale before accumulation so boundary-sized results cannot overflow.
|
||||
checksum += Math.ScaleB(result.High, -1023) + Math.ScaleB(result.Low, -1023);
|
||||
}
|
||||
return GC.GetAllocatedBytesForCurrentThread() - before;
|
||||
}
|
||||
|
||||
private static void CheckBoundary(DoubleDouble value, double high, double low)
|
||||
{
|
||||
BitConverter.DoubleToInt64Bits(value.High).ShouldBe(BitConverter.DoubleToInt64Bits(high));
|
||||
BitConverter.DoubleToInt64Bits(value.Low).ShouldBe(BitConverter.DoubleToInt64Bits(low));
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void SpecialValueMatrixMatchesBinary64IncludingZeroSigns()
|
||||
{
|
||||
|
||||
Reference in New Issue
Block a user