485 lines
20 KiB
C#
485 lines
20 KiB
C#
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, "+")]
|
|
[InlineData(1.0, "-")]
|
|
[InlineData(-1.0, "-")]
|
|
public void AdditionBelowOverflowMidpointRemainsFinite(double sign, string operation)
|
|
{
|
|
// Exact magnitude = MaxValue + 2^970 - 2^916, strictly below
|
|
// the binary64 overflow midpoint. The right factory call reduces to
|
|
// (BitDecrement(2^969), 0), so scalar overloads share this case.
|
|
DoubleDouble left = DoubleDouble.FromComponents(sign * double.MaxValue, sign * Math.ScaleB(1.0, 969));
|
|
DoubleDouble right = DoubleDouble.FromComponents(sign * Math.ScaleB(1.0, 969), -sign * Math.ScaleB(1.0, 916));
|
|
if (operation == "-")
|
|
{
|
|
right = -right;
|
|
}
|
|
Rational expected = Expected(Exact(left), Exact(right), operation);
|
|
BelowOverflowMidpoint(expected).ShouldBeTrue();
|
|
// A canonical finite pair meets the requested accuracy; infinity is
|
|
// not forced by the representational limit or the error contract.
|
|
DoubleDouble finiteWitness = DoubleDouble.FromComponents(sign * double.MaxValue,
|
|
sign * Math.BitDecrement(Math.ScaleB(1.0, 970)));
|
|
AssertAccurate(finiteWitness, expected, "finite witness");
|
|
AssertOperation(left, right, operation);
|
|
}
|
|
|
|
[Theory]
|
|
[InlineData(1.0, "+", false)]
|
|
[InlineData(-1.0, "+", false)]
|
|
[InlineData(1.0, "+", true)]
|
|
[InlineData(-1.0, "+", true)]
|
|
[InlineData(1.0, "-", false)]
|
|
[InlineData(-1.0, "-", false)]
|
|
public void ScalarAdditionBelowOverflowMidpointRemainsFinite(double sign, string operation, bool scalarLeft)
|
|
{
|
|
DoubleDouble pair = DoubleDouble.FromComponents(sign * double.MaxValue, sign * Math.ScaleB(1.0, 969));
|
|
double scalar = sign * Math.BitDecrement(Math.ScaleB(1.0, 969));
|
|
if (operation == "-")
|
|
{
|
|
scalar = -scalar;
|
|
}
|
|
Rational expected = Expected(Exact(pair), Exact(scalar), operation);
|
|
DoubleDouble actual = operation == "-" ? pair - scalar : scalarLeft ? scalar + pair : pair + scalar;
|
|
AssertAccurate(actual, expected, Describe(pair, new DoubleDouble(scalar), operation));
|
|
}
|
|
|
|
[Theory]
|
|
[InlineData(1.0)]
|
|
[InlineData(-1.0)]
|
|
public void MultiplicationBelowOverflowMidpointRemainsFinite(double sign)
|
|
{
|
|
// Exact magnitude = MaxValue + 2^970 - 3*2^914.
|
|
DoubleDouble left = DoubleDouble.FromComponents(sign * double.MaxValue, sign * Math.ScaleB(1.0, 969));
|
|
DoubleDouble right = DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -55));
|
|
BelowOverflowMidpoint(Exact(left) * Exact(right)).ShouldBeTrue();
|
|
AssertOperation(left, right, "*");
|
|
}
|
|
|
|
[Theory]
|
|
[InlineData(1.0)]
|
|
[InlineData(-1.0)]
|
|
public void DivisionWithLowNumeratorBelowOverflowMidpointRemainsFinite(double sign)
|
|
{
|
|
DoubleDouble left = DoubleDouble.FromComponents(sign * double.MaxValue, sign * Math.ScaleB(1.0, 969));
|
|
DoubleDouble right = DoubleDouble.FromComponents(1.0, -Math.ScaleB(1.0, -55));
|
|
Rational expected = Exact(left) / Exact(right);
|
|
BelowOverflowMidpoint(expected).ShouldBeTrue();
|
|
AssertAccurate(left / right, expected, Describe(left, right, "/"));
|
|
}
|
|
|
|
[Theory]
|
|
[InlineData(1.0, false)]
|
|
[InlineData(-1.0, false)]
|
|
[InlineData(1.0, true)]
|
|
[InlineData(-1.0, true)]
|
|
public void DivisionBelowOverflowMidpointRemainsFinite(double sign, bool scalarLeft)
|
|
{
|
|
DoubleDouble left = new(sign * double.MaxValue);
|
|
DoubleDouble right = DoubleDouble.FromComponents(1.0, -Math.ScaleB(1.0, -54));
|
|
Rational expected = Exact(left) / Exact(right);
|
|
BelowOverflowMidpoint(expected).ShouldBeTrue();
|
|
DoubleDouble actual = scalarLeft ? (sign * double.MaxValue) / right : left / right;
|
|
AssertAccurate(actual, expected, Describe(left, right, "/"));
|
|
}
|
|
|
|
[Fact]
|
|
public void FactoryAvoidsIntermediateOverflowForFiniteOppositeSignSums()
|
|
{
|
|
// With U = 2^971 and M = MaxValue, M - 1.5U rounds to M - U.
|
|
// Smaller-first TwoSum computes (M - U) - (-1.5U) = M + 0.5U,
|
|
// which overflows even though the original exact sum is finite.
|
|
foreach (double sign in new[] { -1.0, 1.0 })
|
|
{
|
|
double small = -sign * Math.ScaleB(3.0, 970);
|
|
double large = sign * double.MaxValue;
|
|
foreach (DoubleDouble actual in new[] { DoubleDouble.FromComponents(small, large),
|
|
DoubleDouble.FromComponents(large, small) })
|
|
{
|
|
actual.High.ShouldBe(sign * Math.BitDecrement(double.MaxValue));
|
|
actual.Low.ShouldBe(-sign * Math.ScaleB(1.0, 970));
|
|
Exact(actual).CompareTo(Exact(small) + Exact(large)).ShouldBe(0);
|
|
AssertNormalized(actual);
|
|
}
|
|
}
|
|
}
|
|
|
|
[Fact]
|
|
public void FactoryPreservesExactFiniteSumsAcrossExponentBoundaries()
|
|
{
|
|
double[] components =
|
|
[
|
|
0.0, -0.0, double.Epsilon, -double.Epsilon,
|
|
Math.BitDecrement(Math.ScaleB(1.0, -1022)), Math.ScaleB(1.0, -1022),
|
|
Math.ScaleB(1.0, -969), Math.ScaleB(1.0, -53), 1.0,
|
|
Math.BitIncrement(1.0), Math.ScaleB(1.0, 970), double.MaxValue,
|
|
Math.ScaleB(3.0, 970), -Math.ScaleB(3.0, 970),
|
|
-Math.ScaleB(1.0, -1022), -1.0, -double.MaxValue
|
|
];
|
|
foreach (double high in components)
|
|
{
|
|
foreach (double low in components)
|
|
{
|
|
Rational expected = Exact(high) + Exact(low);
|
|
if (!BelowOverflowMidpoint(expected))
|
|
{
|
|
continue;
|
|
}
|
|
DoubleDouble actual = DoubleDouble.FromComponents(high, low);
|
|
Exact(actual).CompareTo(expected).ShouldBe(0, $"factory ({high:R}, {low:R})");
|
|
AssertNormalized(actual);
|
|
DoubleDouble repeated = DoubleDouble.FromComponents(actual.High, actual.Low);
|
|
repeated.Equals(actual).ShouldBeTrue();
|
|
repeated.GetHashCode().ShouldBe(actual.GetHashCode());
|
|
}
|
|
}
|
|
}
|
|
|
|
[Theory]
|
|
[InlineData("+")]
|
|
[InlineData("-")]
|
|
[InlineData("*")]
|
|
[InlineData("/")]
|
|
public void ArithmeticAcrossFastPathTransitionsMeetsExactRationalBound(string operation)
|
|
{
|
|
int[] exponents = [-1074, -1022, -970, -901, -900, -899, -451, -450, -449,
|
|
-54, -1, 0, 1, 54, 449, 450, 451, 899, 900, 901, 969, 1020, 1021, 1023];
|
|
Random random = new(0x5eed);
|
|
foreach (int leftExponent in exponents)
|
|
{
|
|
foreach (int rightExponent in exponents)
|
|
{
|
|
for (int sample = 0; sample < 4; ++sample)
|
|
{
|
|
DoubleDouble left = Sample(random, leftExponent);
|
|
DoubleDouble right = Sample(random, rightExponent);
|
|
AssertOperation(left, right, operation);
|
|
// Exercise scalar overloads independently, not via equality
|
|
// with the corresponding potentially faulty DD operation.
|
|
AssertScalarOperations(left, right.High, operation);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
[Fact]
|
|
public void CancellationAcrossBinadesRetainsSmallResiduals()
|
|
{
|
|
int[] exponents = [-1022, -969, -450, 0, 450, 969, 1020, 1023];
|
|
foreach (int exponent in exponents)
|
|
{
|
|
double high = Math.ScaleB(1.0, exponent);
|
|
foreach (int gap in new[] { 53, 54, 105, 106, 107, 200, 1000 })
|
|
{
|
|
double low = Math.ScaleB(1.0, exponent - gap);
|
|
DoubleDouble left = DoubleDouble.FromComponents(high, low);
|
|
DoubleDouble right = DoubleDouble.FromComponents(-high, Math.ScaleB(1.0, exponent - gap - 54));
|
|
Rational expected = Exact(left) + Exact(right);
|
|
AssertAccurate(left + right, expected, Describe(left, right, "+"));
|
|
AssertAccurate(right + left, expected, Describe(right, left, "+"));
|
|
AssertAccurate(left - (-right), expected, Describe(left, -right, "-"));
|
|
}
|
|
}
|
|
}
|
|
|
|
[Fact]
|
|
public void FiniteComparisonsAgreeWithExactValuesAtAdjacentHighMidpoints()
|
|
{
|
|
List<DoubleDouble> values = [new(0.0), new(-0.0)];
|
|
foreach (int exponent in new[] { -1022, -970, -450, 0, 450, 970, 1023 })
|
|
{
|
|
double high = Math.ScaleB(1.0, exponent);
|
|
double adjacent = Math.BitIncrement(high);
|
|
double midpointLow = Math.ScaleB(1.0, exponent - 53);
|
|
foreach (double sign in new[] { -1.0, 1.0 })
|
|
{
|
|
values.Add(new DoubleDouble(sign * high));
|
|
values.Add(DoubleDouble.FromComponents(sign * high, sign * midpointLow));
|
|
values.Add(DoubleDouble.FromComponents(sign * adjacent, -sign * midpointLow));
|
|
values.Add(DoubleDouble.FromComponents(sign * high, sign * Math.BitDecrement(midpointLow)));
|
|
values.Add(DoubleDouble.FromComponents(sign * high, sign * Math.BitIncrement(midpointLow)));
|
|
}
|
|
}
|
|
foreach (DoubleDouble left in values)
|
|
{
|
|
foreach (DoubleDouble right in values)
|
|
{
|
|
int order = Exact(left).CompareTo(Exact(right));
|
|
string context = Describe(left, right, "compare");
|
|
Math.Sign(left.CompareTo(right)).ShouldBe(Math.Sign(order), context);
|
|
(left < right).ShouldBe(order < 0, context);
|
|
(left > right).ShouldBe(order > 0, context);
|
|
(left <= right).ShouldBe(order <= 0, context);
|
|
(left >= right).ShouldBe(order >= 0, context);
|
|
(left == right).ShouldBe(order == 0, context);
|
|
(left != right).ShouldBe(order != 0, context);
|
|
left.Equals(right).ShouldBe(order == 0, context);
|
|
if (order == 0)
|
|
{
|
|
left.GetHashCode().ShouldBe(right.GetHashCode(), context);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
[Fact]
|
|
public void NonzeroUnderflowRetainsResultSign()
|
|
{
|
|
foreach (double sign in new[] { -1.0, 1.0 })
|
|
{
|
|
DoubleDouble tiny = new(sign * double.Epsilon);
|
|
DoubleDouble[] zeros = [tiny * 0.25, 0.25 * tiny, tiny / 4.0,
|
|
tiny * new DoubleDouble(0.25), tiny / new DoubleDouble(4.0),
|
|
(sign * double.Epsilon) / new DoubleDouble(4.0)];
|
|
foreach (DoubleDouble zero in zeros)
|
|
{
|
|
BitConverter.DoubleToInt64Bits(zero.High).ShouldBe(sign < 0.0 ? long.MinValue : 0L);
|
|
BitConverter.DoubleToInt64Bits(zero.Low).ShouldBe(0L);
|
|
}
|
|
}
|
|
}
|
|
|
|
private static DoubleDouble Sample(Random random, int exponent)
|
|
{
|
|
double sign = random.Next(2) == 0 ? -1.0 : 1.0;
|
|
double high = Math.ScaleB(sign * (1.0 + random.NextDouble()), exponent);
|
|
double low = Math.ScaleB((random.NextDouble() * 2.0) - 1.0, exponent - random.Next(53, 121));
|
|
return DoubleDouble.FromComponents(high, low);
|
|
}
|
|
|
|
private static void AssertOperation(DoubleDouble left, DoubleDouble right, string operation)
|
|
{
|
|
Rational expected = Expected(Exact(left), Exact(right), operation);
|
|
if (!BelowOverflowMidpoint(expected))
|
|
{
|
|
return;
|
|
}
|
|
DoubleDouble actual = operation switch
|
|
{
|
|
"+" => left + right,
|
|
"-" => left - right,
|
|
"*" => left * right,
|
|
"/" => left / right,
|
|
_ => throw new ArgumentOutOfRangeException(nameof(operation))
|
|
};
|
|
AssertAccurate(actual, expected, Describe(left, right, operation));
|
|
AssertNormalized(actual);
|
|
}
|
|
|
|
private static void AssertScalarOperations(DoubleDouble left, double right, string operation)
|
|
{
|
|
Rational forward = Expected(Exact(left), Exact(right), operation);
|
|
Rational reverse = Expected(Exact(right), Exact(left), operation);
|
|
if (BelowOverflowMidpoint(forward))
|
|
{
|
|
DoubleDouble actual = operation switch
|
|
{
|
|
"+" => left + right,
|
|
"-" => left - right,
|
|
"*" => left * right,
|
|
"/" => left / right,
|
|
_ => throw new ArgumentOutOfRangeException(nameof(operation))
|
|
};
|
|
AssertAccurate(actual, forward, Describe(left, new DoubleDouble(right), operation));
|
|
AssertNormalized(actual);
|
|
}
|
|
if (BelowOverflowMidpoint(reverse))
|
|
{
|
|
DoubleDouble actual = operation switch
|
|
{
|
|
"+" => right + left,
|
|
"-" => right - left,
|
|
"*" => right * left,
|
|
"/" => right / left,
|
|
_ => throw new ArgumentOutOfRangeException(nameof(operation))
|
|
};
|
|
AssertAccurate(actual, reverse, Describe(new DoubleDouble(right), left, operation));
|
|
AssertNormalized(actual);
|
|
}
|
|
}
|
|
|
|
private static Rational Expected(Rational left, Rational right, string operation)
|
|
{
|
|
return operation switch
|
|
{
|
|
"+" => left + right,
|
|
"-" => left - right,
|
|
"*" => left * right,
|
|
"/" => left / right,
|
|
_ => throw new ArgumentOutOfRangeException(nameof(operation))
|
|
};
|
|
}
|
|
|
|
private static bool BelowOverflowMidpoint(Rational value)
|
|
{
|
|
return value.Abs().CompareTo(Exact(double.MaxValue) + Exact(Math.ScaleB(1.0, 970))) < 0;
|
|
}
|
|
|
|
private static void AssertNormalized(DoubleDouble value)
|
|
{
|
|
double.IsFinite(value.High).ShouldBeTrue();
|
|
double.IsFinite(value.Low).ShouldBeTrue();
|
|
(value.High + value.Low).ShouldBe(value.High);
|
|
if (value.Low == 0.0)
|
|
{
|
|
BitConverter.DoubleToInt64Bits(value.Low).ShouldBe(0L);
|
|
}
|
|
}
|
|
|
|
private static string Describe(DoubleDouble left, DoubleDouble right, string operation)
|
|
{
|
|
return $"({left.High:R}, {left.Low:R}) {operation} ({right.High:R}, {right.Low:R})";
|
|
}
|
|
|
|
private static void AssertAccurate(DoubleDouble actual, Rational expected, string context)
|
|
{
|
|
string diagnostic = $"{context}: actual ({actual.High:R}, {actual.Low:R})";
|
|
double.IsFinite(actual.High).ShouldBeTrue(diagnostic);
|
|
double.IsFinite(actual.Low).ShouldBeTrue(diagnostic);
|
|
Rational error = (Exact(actual) - expected).Abs();
|
|
// Conservative contract, not a correct-rounding assertion. All
|
|
// comparisons, including the subnormal floor, use exact rationals.
|
|
Rational tolerance = (expected.Abs() * new Rational(1, BigInteger.One << 100))
|
|
+ Exact(double.Epsilon);
|
|
error.CompareTo(tolerance).ShouldBeLessThanOrEqualTo(0, diagnostic);
|
|
}
|
|
|
|
private static Rational Exact(DoubleDouble value)
|
|
{
|
|
return Exact(value.High) + Exact(value.Low);
|
|
}
|
|
|
|
private static Rational Exact(double value)
|
|
{
|
|
// Decode IEEE-754 directly; no production helpers or conversions.
|
|
double.IsFinite(value).ShouldBeTrue();
|
|
ulong bits = BitConverter.DoubleToUInt64Bits(value);
|
|
int biasedExponent = (int)((bits >> 52) & 0x7ff);
|
|
BigInteger significand = bits & 0x000f_ffff_ffff_ffffUL;
|
|
int exponent = -1074;
|
|
if (biasedExponent != 0)
|
|
{
|
|
significand += BigInteger.One << 52;
|
|
exponent = biasedExponent - 1075;
|
|
}
|
|
if ((bits >> 63) != 0)
|
|
{
|
|
significand = -significand;
|
|
}
|
|
return exponent >= 0
|
|
? new Rational(significand << exponent, BigInteger.One)
|
|
: new Rational(significand, BigInteger.One << -exponent);
|
|
}
|
|
|
|
private readonly struct Rational
|
|
{
|
|
private readonly BigInteger _numerator;
|
|
private readonly BigInteger _denominator;
|
|
|
|
public Rational(BigInteger numerator, BigInteger denominator)
|
|
{
|
|
if (denominator.IsZero)
|
|
{
|
|
throw new DivideByZeroException();
|
|
}
|
|
BigInteger divisor = BigInteger.GreatestCommonDivisor(numerator, denominator);
|
|
_numerator = numerator / divisor * denominator.Sign;
|
|
_denominator = BigInteger.Abs(denominator / divisor);
|
|
}
|
|
|
|
public Rational Abs()
|
|
{
|
|
return new Rational(BigInteger.Abs(_numerator), _denominator);
|
|
}
|
|
|
|
public int CompareTo(Rational other)
|
|
{
|
|
return (_numerator * other._denominator).CompareTo(other._numerator * _denominator);
|
|
}
|
|
|
|
public static Rational operator +(Rational left, Rational right)
|
|
{
|
|
return new Rational((left._numerator * right._denominator) + (right._numerator * left._denominator),
|
|
left._denominator * right._denominator);
|
|
}
|
|
|
|
public static Rational operator -(Rational value)
|
|
{
|
|
return new Rational(-value._numerator, value._denominator);
|
|
}
|
|
|
|
public static Rational operator -(Rational left, Rational right)
|
|
{
|
|
return left + (-right);
|
|
}
|
|
|
|
public static Rational operator *(Rational left, Rational right)
|
|
{
|
|
return new Rational(left._numerator * right._numerator, left._denominator * right._denominator);
|
|
}
|
|
|
|
public static Rational operator /(Rational left, Rational right)
|
|
{
|
|
return new Rational(left._numerator * right._denominator, left._denominator * right._numerator);
|
|
}
|
|
}
|
|
}
|