using System.Numerics; using Shouldly; using Xunit; namespace Just.PreciseMath.Tests; public class DoubleDoubleArithmeticTests { [Theory] [InlineData(6.0, 2.0, 8.0, 4.0, 12.0, 3.0)] [InlineData(-6.0, 2.0, -4.0, -8.0, -12.0, -3.0)] [InlineData(6.0, -2.0, 4.0, 8.0, -12.0, -3.0)] [InlineData(-6.0, -2.0, -8.0, -4.0, 12.0, 3.0)] [InlineData(1.5, 0.5, 2.0, 1.0, 0.75, 3.0)] [InlineData(-1.5, 0.5, -1.0, -2.0, -0.75, -3.0)] [InlineData(0.75, 1.5, 2.25, -0.75, 1.125, 0.5)] [InlineData(0.0, 2.0, 2.0, -2.0, 0.0, 0.0)] [InlineData(2.0, 2.0, 4.0, 0.0, 4.0, 1.0)] public void BasicArithmeticHasExactDyadicResults(double left, double right, double sum, double difference, double product, double quotient) { // Small integers and binary fractions: each expected result is an exact // rational representable in binary64, so no accuracy tolerance is needed. DoubleDouble a = new(left); DoubleDouble b = new(right); CheckBoundary(a + b, sum, 0.0); CheckBoundary(a + right, sum, 0.0); CheckBoundary(left + b, sum, 0.0); CheckBoundary(a - b, difference, 0.0); CheckBoundary(a - right, difference, 0.0); CheckBoundary(left - b, difference, 0.0); CheckBoundary(a * b, product, 0.0); CheckBoundary(a * right, product, 0.0); CheckBoundary(left * b, product, 0.0); CheckBoundary(a / b, quotient, 0.0); CheckBoundary(a / right, quotient, 0.0); CheckBoundary(left / b, quotient, 0.0); } [Theory] [InlineData(30.0, 0.5235987755982989, -5.360408832255455e-17)] [InlineData(45.0, 0.7853981633974483, 3.061616997868383e-17)] [InlineData(60.0, 1.0471975511965979, -1.072081766451091e-16)] [InlineData(90.0, 1.5707963267948966, 6.123233995736766e-17)] [InlineData(180.0, 3.141592653589793, 1.2246467991473532e-16)] [InlineData(360.0, 6.283185307179586, 2.4492935982947064e-16)] public void AngleConversionConstantsRetainExtendedPrecision(double degrees, double radiansHigh, double radiansLow) { // Independent precomputed pi fractions from ReferenceData/generate_constants.py. // Check each direction against its own reference, not a computed round trip. foreach (double sign in new[] { -1.0, 1.0 }) { double signedDegrees = sign * degrees; BigInteger radiansUnits = Units(sign * radiansHigh) + Units(sign * radiansLow); AssertRelative(signedDegrees * DoubleDouble.DegToRad, radiansUnits, BigInteger.One); AssertRelative(DoubleDouble.DegToRad * signedDegrees, radiansUnits, BigInteger.One); AssertRelative(new DoubleDouble(signedDegrees) * DoubleDouble.DegToRad, radiansUnits, BigInteger.One); DoubleDouble radians = DoubleDouble.FromComponents(sign * radiansHigh, sign * radiansLow); AssertRelative(radians * DoubleDouble.RadToDeg, Units(signedDegrees), BigInteger.One); AssertRelative(DoubleDouble.RadToDeg * radians, Units(signedDegrees), BigInteger.One); } } // Generated by: python3 1-tests/Just.PreciseMath.Tests/ReferenceData/generate_irrational_arithmetic.py // Constants: Python Decimal at 160 digits, checked again at 240 digits; // pi uses Machin's formula, e/roots/ln use Decimal's exp/sqrt/ln, phi=(1+sqrt(5))/2. // Expected operations use exact Fraction sums of the stored input components, // then round high and residual separately to binary64. A scalar input has no low. // Column pairs: left, right, sum, difference, product, quotient. [Theory] // pi, e: DD/DD [InlineData("DD/DD", 3.141592653589793, 1.2246467991473532e-16, 2.718281828459045, 1.4456468917292502e-16, 5.859874482048839, -1.7705984076240228e-16, 0.423310825130748, -2.2100009258189695e-17, 8.539734222673568, -6.773815290502424e-16, 1.1557273497909217, -1.3998972600526045e-17)] // pi, e: DD/double [InlineData("DD/double", 3.141592653589793, 1.2246467991473532e-16, 2.718281828459045, 0.0, 5.859874482048839, -3.2162452993532727e-16, 0.42331082513074814, 1.1442377452219667e-17, 8.539734222673566, 6.44811944875855e-16, 1.1557273497909217, 4.746535510161172e-17)] // pi, e: double/DD [InlineData("double/DD", 3.141592653589793, 0.0, 2.718281828459045, 1.4456468917292502e-16, 5.859874482048839, -2.995245206771376e-16, 0.42331082513074786, 2.1968764520848465e-17, 8.539734222673566, 7.660817963097297e-16, 1.1557273497909217, -5.905121061079843e-17)] // sqrt2, sqrt3: DD/DD [InlineData("DD/DD", 1.4142135623730951, -9.667293313452913e-17, 1.7320508075688772, 1.0035084221806903e-16, 3.1462643699419726, -2.1836669584149143e-16, -0.31783724519578227, 2.5020829572433146e-17, 2.449489742783178, 2.168616518103246e-16, 0.816496580927726, -1.7276510382355668e-18)] // sqrt2, sqrt3: DD/double [InlineData("DD/double", 1.4142135623730951, -9.667293313452913e-17, 1.7320508075688772, 0.0, 3.146264369941972, 1.2537167179050217e-16, -0.31783724519578216, 1.434936932798652e-17, 2.449489742783178, 7.494412974996883e-17, 0.816496580927726, 4.5578189648549696e-17)] // sqrt2, sqrt3: double/DD [InlineData("double/DD", 1.4142135623730951, 0.0, 1.7320508075688772, 1.0035084221806903e-16, 3.1462643699419726, -1.2169376270696227e-16, -0.31783724519578216, 1.0671460244446626e-17, 2.4494897427831783, -5.978512613402474e-17, 0.816496580927726, 5.408649293033552e-17)] // ln2, phi: DD/DD [InlineData("DD/DD", 0.6931471805599453, 2.3190468138462996e-17, 1.618033988749895, -5.432115203682506e-17, 2.3111811693098403, -1.421529863608777e-16, -0.9248868081899495, -3.35106822872276e-17, 1.121535697352152, -9.053143373999594e-17, 0.4283885167922066, -2.699599415943276e-18)] // ln2, phi: DD/double [InlineData("DD/double", 0.6931471805599453, 2.3190468138462996e-17, 1.618033988749895, 0.0, 2.3111811693098403, -8.783183432405266e-17, -0.9248868081899496, 2.3190468138462996e-17, 1.121535697352152, -5.2878880360902515e-17, 0.4283885167922066, -1.708159504353726e-17)] // ln2, phi: double/DD [InlineData("double/DD", 0.6931471805599453, 0.0, 1.618033988749895, -5.432115203682506e-17, 2.3111811693098403, -1.653434544993407e-16, -0.9248868081899496, 5.432115203682506e-17, 1.1215356973521518, 9.399020552198074e-17, 0.4283885167922066, -1.703209694053491e-17)] public void BasicArithmeticMatchesPrecomputedIrrationalResults(string overload, double leftHigh, double leftLow, double rightHigh, double rightLow, double sumHigh, double sumLow, double differenceHigh, double differenceLow, double productHigh, double productLow, double quotientHigh, double quotientLow) { foreach (double sign in new[] { -1.0, 1.0 }) { DoubleDouble left = DoubleDouble.FromComponents(sign * leftHigh, sign * leftLow); DoubleDouble right = DoubleDouble.FromComponents(sign * rightHigh, sign * rightLow); (DoubleDouble sum, DoubleDouble difference, DoubleDouble product, DoubleDouble quotient) = overload switch { "DD/DD" => (left + right, left - right, left * right, left / right), "DD/double" => (left + right.High, left - right.High, left * right.High, left / right.High), "double/DD" => (left.High + right, left.High - right, left.High * right, left.High / right), _ => throw new ArgumentOutOfRangeException(nameof(overload)) }; // Negating both inputs negates sum/difference but not product/quotient. // Compare complete expansions at the established 2^-100 relative + epsilon bound, // not individual component equality: correct rounding is not guaranteed. // The generator checks reference-rounding error is below 2^-105 relative. AssertRelative(sum, Units(sign * sumHigh) + Units(sign * sumLow), BigInteger.One); AssertRelative(difference, Units(sign * differenceHigh) + Units(sign * differenceLow), BigInteger.One); AssertRelative(product, Units(productHigh) + Units(productLow), BigInteger.One); AssertRelative(quotient, Units(quotientHigh) + Units(quotientLow), BigInteger.One); } } [Theory] [InlineData("+", "DD/DD")] [InlineData("-", "DD/DD")] [InlineData("*", "DD/DD")] [InlineData("/", "DD/DD")] [InlineData("+", "DD/double")] [InlineData("-", "DD/double")] [InlineData("*", "DD/double")] [InlineData("/", "DD/double")] [InlineData("+", "double/DD")] [InlineData("-", "double/DD")] [InlineData("*", "double/DD")] [InlineData("/", "double/DD")] public void GeneralArithmeticMeetsExactRationalBound(string operation, string overload) { // Fixed seed, with dense, sparse and zero residuals independent of the // high sign. Include ordinary exponents as well as dispatch/range edges. // Every pair is checked in both orders; overflow cases are asserted, not skipped. Random random = new(65537); int[] exponents = [-1074, -1022, -969, -451, -450, -1, 0, 1, 450, 451, 900, 1020, 1021, 1023]; foreach (int leftExponent in exponents) { foreach (int rightExponent in exponents) { for (int sample = 0; sample < 4; ++sample) { DoubleDouble left = GeneralArithmeticSample(random, leftExponent, sample); DoubleDouble right = GeneralArithmeticSample(random, rightExponent, (sample + 1) % 4); AssertGeneralArithmetic(left, right, operation, overload); AssertGeneralArithmetic(right, left, operation, overload); } } } for (int sample = 0; sample < 512; ++sample) { int exponent = random.Next(-450, 451); DoubleDouble left = GeneralArithmeticSample(random, exponent, sample % 4); DoubleDouble right = GeneralArithmeticSample(random, random.Next(-450, 451), (sample + 1) % 4); AssertGeneralArithmetic(left, right, operation, overload); AssertGeneralArithmetic(right, left, operation, overload); // Correlate highs to force cancellation rather than hoping random // independent values happen to exercise it. Keep different low terms. DoubleDouble neighbor = DoubleDouble.FromComponents(left.High, -left.Low); AssertGeneralArithmetic(left, neighbor, operation, overload); AssertGeneralArithmetic(left, -neighbor, operation, overload); } } private static DoubleDouble GeneralArithmeticSample(Random random, int exponent, int residualKind) { double sign = random.Next(2) == 0 ? -1.0 : 1.0; double high = Math.ScaleB(sign * (1.0 + (0.75 * random.NextDouble())), exponent); double low = residualKind switch { 0 => 0.0, 1 => Math.ScaleB(random.NextDouble() - 0.5, exponent - 53), 2 => Math.ScaleB(random.NextDouble() - 0.5, exponent - 106), 3 => random.Next(2) == 0 ? -double.Epsilon : double.Epsilon, _ => throw new ArgumentOutOfRangeException(nameof(residualKind)) }; // At the subnormal floor a residual can cancel the high completely. // Use a scalar there so this finite-input matrix never divides by zero; // the separate special-value matrix covers zero denominators and NaNs. if (exponent < -1022) { low = 0.0; } DoubleDouble result = DoubleDouble.FromComponents(high, low); double.IsFinite(result.High).ShouldBeTrue(); result.High.ShouldNotBe(0.0); return result; } private static void AssertGeneralArithmetic(DoubleDouble left, DoubleDouble right, string operation, string overload) { DoubleDouble actual = (operation, overload) switch { ("+", "DD/DD") => left + right, ("-", "DD/DD") => left - right, ("*", "DD/DD") => left * right, ("/", "DD/DD") => left / right, ("+", "DD/double") => left + right.High, ("-", "DD/double") => left - right.High, ("*", "DD/double") => left * right.High, ("/", "DD/double") => left / right.High, ("+", "double/DD") => left.High + right, ("-", "double/DD") => left.High - right, ("*", "double/DD") => left.High * right, ("/", "double/DD") => left.High / right, _ => throw new ArgumentOutOfRangeException(nameof(operation)) }; // 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); BigInteger y = overload == "DD/double" ? Units(right.High) : Units(right); (BigInteger numerator, BigInteger denominator) = operation switch { "+" => (x + y, BigInteger.One), "-" => (x - y, BigInteger.One), "*" => (x * y, BigInteger.One << 1074), "/" => (x << 1074, y), _ => throw new ArgumentOutOfRangeException(nameof(operation)) }; if (denominator.Sign < 0) { numerator = -numerator; denominator = -denominator; } string context = $"{overload}: ({left.High:R}, {left.Low:R}) {operation} ({right.High:R}, {right.Low:R}); " + $"actual ({actual.High:R}, {actual.Low:R})"; // Nearest-even binary64 overflow begins at 2^1024 - 2^970. BigInteger overflowUnits = Units(double.MaxValue) + Units(Math.ScaleB(1.0, 970)); if (BigInteger.Abs(numerator) >= overflowUnits * denominator) { actual.High.ShouldBe(numerator.Sign < 0 ? double.NegativeInfinity : double.PositiveInfinity, context); BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L, context); return; } double.IsFinite(actual.High).ShouldBeTrue(context); double.IsFinite(actual.Low).ShouldBeTrue(context); BigInteger error = BigInteger.Abs((Units(actual) * denominator) - numerator); // Same conservative 2^-100 relative + epsilon contract as the existing // suite, cross-multiplied exactly to retain all low-component information. (error <= (BigInteger.Abs(numerator) >> 100) + denominator).ShouldBeTrue(context); (actual.High + actual.Low).ShouldBe(actual.High, context); if (actual.Low == 0.0) { BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L, context); } if (actual.High == 0.0) { // Exact cancellation is +0; a nonzero underflow keeps its sign. BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(numerator.Sign < 0 ? long.MinValue : 0L, context); } } [Fact] public void CancellationRetainsBothLowSumTerms() { double small = Math.ScaleB(1.0, -54); double tiny = Math.ScaleB(1.0, -108); DoubleDouble left = DoubleDouble.FromComponents(1.0, small); DoubleDouble right = DoubleDouble.FromComponents(-1.0, tiny); Check(left + right, small, tiny); Check(right + left, small, tiny); Check(left - (-right), small, tiny); Check(+left, 1.0, small); Check(-left, -1.0, -small); Check(left - left, 0.0, 0.0); } [Fact] public void ScalarOverloadsPreserveOperandOrderAndResiduals() { DoubleDouble value = DoubleDouble.FromComponents(2.0, Math.ScaleB(1.0, -80)); Check(3.0 - value, 1.0, -value.Low); Check(value - 3.0, -1.0, value.Low); Check(value + 3.0, 5.0, value.Low); Check(3.0 + value, 5.0, value.Low); Check(value * 2.0, 4.0, 2.0 * value.Low); Check(2.0 * value, 4.0, 2.0 * value.Low); Check(value / 2.0, 1.0, value.Low / 2.0); 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() { // Review-1 §1: the former operator -(double, DoubleDouble) returned arg - lvalue, // so 3.0 - DD(2.0) produced -1 instead of 1. Check(3.0 - new DoubleDouble(2.0), 1.0, 0.0); Check(new DoubleDouble(2.0) - 3.0, -1.0, 0.0); Check(3.0 - DoubleDouble.FromComponents(2.0, Math.ScaleB(1.0, -80)), 1.0, -Math.ScaleB(1.0, -80)); Check(DoubleDouble.FromComponents(2.0, Math.ScaleB(1.0, -80)) - 3.0, -1.0, Math.ScaleB(1.0, -80)); } [Fact] public void DivisionByAValueCarriedOnlyInTheLowComponentIsFinite() { // Review-2 §5: with high == 0 and low != 0 the former division returned // Infinity because it divided by the zero high component. The public factory // now folds such a pair into its high component, so the quotient is finite. DoubleDouble denominator = DoubleDouble.FromComponents(0.0, 1.0); Check(denominator, 1.0, 0.0); Check(new DoubleDouble(4.0) / denominator, 4.0, 0.0); Check(4.0 / denominator, 4.0, 0.0); Check(new DoubleDouble(4.0) / DoubleDouble.FromComponents(0.0, -2.0), -2.0, 0.0); } [Fact] public void ScalarCancellationPreservesTheRemainingExpansionInBothOrders() { // At a normal binade boundary, 2^e - BitDecrement(2^e) = 2^(e-53). // The low input becomes the representable residual of that exact difference. foreach (int exponent in new[] { -967, -900, -450, 0, 450, 969, 1020, 1023 }) { foreach (double sign in new[] { -1.0, 1.0 }) { double high = Math.ScaleB(1.0, exponent); double low = sign * Math.ScaleB(1.0, exponent - 107); double scalar = sign * Math.BitDecrement(high); double difference = sign * Math.ScaleB(1.0, exponent - 53); DoubleDouble value = DoubleDouble.FromComponents(sign * high, low); Check(value - scalar, difference, low); Check(scalar - value, -difference, -low); Check(value + (-scalar), difference, low); Check((-scalar) + value, difference, low); } } } [Fact] public void MixedSpecialValuesIgnoreFiniteResidualsButPreserveResultSigns() { foreach (double sign in new[] { -1.0, 1.0 }) { DoubleDouble value = DoubleDouble.FromComponents(sign, sign * Math.ScaleB(1.0, -54)); foreach (double scalar in new[] { double.NaN, double.NegativeInfinity, double.PositiveInfinity }) { CheckBits(value + scalar, sign + scalar); CheckBits(scalar + value, scalar + sign); CheckBits(value - scalar, sign - scalar); CheckBits(scalar - value, scalar - sign); CheckBits(value * scalar, sign * scalar); CheckBits(scalar * value, scalar * sign); CheckBits(value / scalar, sign / scalar); CheckBits(scalar / value, scalar / sign); } foreach (double zero in new[] { 0.0, -0.0 }) { CheckBits(value * zero, sign * zero); CheckBits(zero * value, zero * sign); CheckBits(value / zero, sign / zero); CheckBits(zero / value, zero / sign); Check(value + zero, value.High, value.Low); Check(zero + value, value.High, value.Low); Check(value - zero, value.High, value.Low); Check(zero - value, -value.High, -value.Low); } } } [Fact] public void ProductAndQuotientRetainExtraPrecision() { // (1 + 2^-52)(1 - 2^-52) = 1 - 2^-104 exactly. Check(new DoubleDouble(1.0 + Math.ScaleB(1.0, -52)) * new DoubleDouble(1.0 - Math.ScaleB(1.0, -52)), 1.0, -Math.ScaleB(1.0, -104)); // Binary expansion of 1/3, rounding high then residual ties-to-even. Check(DoubleDouble.One / 3.0, 0.3333333333333333, 1.850371707708594e-17); } [Fact] public void ScalarProductNormalizesACorrectionBeyondTheHighMidpoint() { // (1 + 2^-53)(1 + 2^-52) = 1 + 3*2^-53 + 2^-105, exactly. // The rounded high advances twice above 1; its residual is still exact. DoubleDouble value = DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -53)); double scalar = Math.BitIncrement(1.0); double high = 1.0 + Math.ScaleB(1.0, -51); double low = -Math.ScaleB(1.0, -53) + Math.ScaleB(1.0, -105); Check(value * scalar, high, low); 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] public void ScalarDivisionRetainsNumeratorAndDenominatorResiduals() { double low = Math.ScaleB(1.0, -80); DoubleDouble value = DoubleDouble.FromComponents(1.0, low); Check(value / 1.0, 1.0, low); Check(value / (-1.0), -1.0, -low); // Compare the reciprocal to its exact rational, not another DD operator. BigInteger numerator = BigInteger.One << 2148; AssertRelative(1.0 / value, numerator, Units(value)); AssertRelative(-1.0 / value, -numerator, Units(value)); } [Fact] public void DivisionNormalizesTheFirstCorrectionBeforeAddingTheLast() { // Exact rational: (1 + 2^-54) / (17/16 - 2^-54). // Independently round that rational to binary64, then round its exact // residual: high bits 3FEE1E1E1E1E1E1F, low bits 3C4FE3A76B2EF2C4. // Adding the last correction to the unnormalized pair instead loses // four low-component ULPs. Common power-of-two scaling preserves the ratio. double expectedHigh = BitConverter.UInt64BitsToDouble(0x3fee_1e1e_1e1e_1e1f); double expectedLow = BitConverter.UInt64BitsToDouble(0x3c4f_e3a7_6b2e_f2c4); foreach (int exponent 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, exponent), leftSign * Math.ScaleB(1.0, exponent - 54)); DoubleDouble right = DoubleDouble.FromComponents(rightSign * Math.ScaleB(1.0625, exponent), -rightSign * Math.ScaleB(1.0, exponent - 54)); double sign = leftSign * rightSign; DoubleDouble actual = left / right; CheckBoundary(actual, sign * expectedHigh, sign * expectedLow); AssertRelative(actual, Units(left) << 1074, Units(right)); } } } } [Fact] public void DivisionSecondResidualPreservesBitsAndRationalAccuracy() { Random random = new(104729); for (int i = 0; i < 4096; ++i) { int leftExponent = random.Next(-450, 451); int rightExponent = random.Next(-450, 451); double leftHigh = Math.ScaleB(1.0 + (0.75 * random.NextDouble()), leftExponent); double rightHigh = Math.ScaleB(1.0 + (0.75 * random.NextDouble()), rightExponent); double leftLow = i % 3 == 0 ? 0.0 : i % 3 == 1 ? double.Epsilon : Math.ScaleB(random.NextDouble() - 0.5, leftExponent - 53); double rightLow = i % 3 == 1 ? 0.0 : i % 3 == 2 ? -double.Epsilon : Math.ScaleB(random.NextDouble() - 0.5, rightExponent - 53); DoubleDouble left = DoubleDouble.FromComponents(i % 2 == 0 ? leftHigh : -leftHigh, leftLow); DoubleDouble right = DoubleDouble.FromComponents(i % 4 < 2 ? rightHigh : -rightHigh, rightLow); AssertDivisionMatchesPrevious(left, right); AssertDivisionMatchesPrevious(right, left); } } [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() { DoubleDouble maximum = new(double.MaxValue); DoubleDouble third = maximum / 3.0; DoubleDouble thirdPair = maximum / new DoubleDouble(3.0); DoubleDouble thirdScalar = double.MaxValue / new DoubleDouble(3.0); AssertRelative(third, Units(maximum), 3); thirdPair.ShouldBe(third); thirdScalar.ShouldBe(third); Check(new DoubleDouble(double.Epsilon) / new DoubleDouble(double.Epsilon), 1.0, 0.0); Check(new DoubleDouble(double.Epsilon) * new DoubleDouble(Math.ScaleB(1.0, 1023)), Math.ScaleB(1.0, -51), 0.0); Check(new DoubleDouble(Math.ScaleB(1.0, -1022)) / 2.0, Math.ScaleB(1.0, -1023), 0.0); Check(maximum * 2.0, double.PositiveInfinity, 0.0); Check(maximum + maximum, double.PositiveInfinity, 0.0); // High-only addition overflows, but the complete sum is exactly MaxValue. DoubleDouble below = DoubleDouble.FromComponents(double.MaxValue, -Math.ScaleB(1.0, 969)); 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)] public void ExactOverflowMidpointStillOverflows(double sign) { DoubleDouble maximum = new(sign * double.MaxValue); double halfUlp = sign * Math.ScaleB(1.0, 970); Check(maximum + halfUlp, sign * double.PositiveInfinity, 0.0); Check(maximum - (-halfUlp), sign * double.PositiveInfinity, 0.0); 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() { double[] values = [0.0, -0.0, 1.0, -1.0, double.PositiveInfinity, double.NegativeInfinity, double.NaN]; foreach (double left in values) { foreach (double right in values) { DoubleDouble a = new(left); DoubleDouble b = new(right); CheckBits(a + b, left + right); CheckBits(a - b, left - right); CheckBits(a * b, left * right); CheckBits(a / b, left / right); CheckBits(a + right, left + right); CheckBits(left + b, left + right); CheckBits(a - right, left - right); CheckBits(left - b, left - right); CheckBits(a * right, left * right); CheckBits(left * b, left * right); CheckBits(a / right, left / right); CheckBits(left / b, left / right); } } } [Fact] public void DeterministicArithmeticMeetsConservativeErrorBound() { Random random = new(1729); for (int i = 0; i < 250; ++i) { DoubleDouble a = DoubleDouble.FromComponents(Math.ScaleB((random.NextDouble() * 2.0) - 1.0, random.Next(-400, 401)), Math.ScaleB(random.NextDouble(), random.Next(-500, -450))); DoubleDouble b = DoubleDouble.FromComponents(Math.ScaleB((random.NextDouble() * 2.0) - 1.0, random.Next(-400, 401)), Math.ScaleB(random.NextDouble(), random.Next(-500, -450))); BigInteger x = Units(a); BigInteger y = Units(b); AssertRelative(a + b, x + y, BigInteger.One); AssertRelative(a - b, x - y, BigInteger.One); AssertRelative(a * b, x * y, BigInteger.One << 1074); AssertRelative(a / b, x << 1074, y); } } private static void AssertDivisionMatchesPrevious(DoubleDouble left, DoubleDouble right) { // Freeze the pre-specialization expression. The unchanged public operators // form the bitwise reference; Units supplies the independent rational oracle. double quotient = left.High / right.High; DoubleDouble remainder = left - (right * quotient); double correction = remainder.High / right.High; remainder -= right * correction; double finalCorrection = remainder.High / right.High; DoubleDouble expected = DoubleDouble.FromComponents(quotient, correction) + finalCorrection; DoubleDouble actual = left / right; CheckBoundary(actual, expected.High, expected.Low); AssertRelative(actual, Units(left) << 1074, Units(right)); } private static void Check(DoubleDouble value, double high, double low) { value.High.ShouldBe(high); value.Low.ShouldBe(low); } private static void CheckBits(DoubleDouble value, double expected) { if (double.IsNaN(expected)) { DoubleDouble.IsNaN(value).ShouldBeTrue(); } else { BitConverter.DoubleToInt64Bits(value.High).ShouldBe(BitConverter.DoubleToInt64Bits(expected)); } BitConverter.DoubleToInt64Bits(value.Low).ShouldBe(0L); } // Independent oracle: every finite binary64 is an integer multiple of 2^-1074. private static BigInteger Units(DoubleDouble value) { return Units(value.High) + Units(value.Low); } private static BigInteger Units(double value) { long bits = BitConverter.DoubleToInt64Bits(value); int exponent = (int)((bits >> 52) & 0x7ff); BigInteger significand = bits & 0xfffffffffffffL; if (exponent != 0) { significand += BigInteger.One << 52; significand <<= exponent - 1; } return bits < 0 ? -significand : significand; } private static void AssertRelative(DoubleDouble actual, BigInteger numerator, BigInteger denominator) { DoubleDouble.IsFinite(actual).ShouldBeTrue(); BigInteger error = BigInteger.Abs((Units(actual) * denominator) - numerator); // <= 2^-100 relative error plus one minimum subnormal (rounding floor). (error <= (BigInteger.Abs(numerator) >> 100) + BigInteger.Abs(denominator)).ShouldBeTrue(); if (actual.High != 0.0) { (Math.Abs(actual.Low) <= Math.ScaleB(1.0, Math.ILogB(actual.High) - 53)).ShouldBeTrue(); } } }