From 4efecbb1bc254710bbeb1977f1e627885951df38 Mon Sep 17 00:00:00 2001 From: just Date: Mon, 14 Sep 2026 21:21:48 +0400 Subject: [PATCH] add precomputed mathematical constants and specialize the second division residual --- .../DoubleDouble.Arithmetic.cs | 10 +- 0-source/Just.PreciseMath/DoubleDouble.cs | 87 ++++- .../Just.PreciseMath/PreciseMathHelper.cs | 29 +- .../ArithmeticRangeTests.cs | 91 ++++- .../DoubleDoubleArithmeticTests.cs | 362 ++++++++++++++++++ .../DoubleDoubleFormattingTests.cs | 6 +- .../DoubleDoubleSpanFormattingTests.cs | 2 +- .../DoubleDoubleTests.cs | 121 +++++- .../ReferenceData/generate_constants.py | 76 ++++ .../generate_irrational_arithmetic.py | 77 ++++ README.md | 38 +- 11 files changed, 875 insertions(+), 24 deletions(-) create mode 100644 1-tests/Just.PreciseMath.Tests/ReferenceData/generate_constants.py create mode 100644 1-tests/Just.PreciseMath.Tests/ReferenceData/generate_irrational_arithmetic.py diff --git a/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs b/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs index 1afd017..90d091e 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs @@ -101,9 +101,13 @@ public readonly partial struct DoubleDouble : 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; - return FromComponents(quotient, correction) + finalCorrection; + double finalRemainder = PreciseMathHelper.SubtractDivisionCorrectionHigh(remainder, right * correction); + double finalCorrection = finalRemainder / right._high; + // The entry guard gives |quotient| <= 2^901. Conservatively bounding + // the first product and subtraction gives |remainder.High| <= 2^457 + // at the first correction, hence |correction| <= 2^908 and a finite sum + // below 2^909 in magnitude. Retain normalization before the final add. + return PreciseMathHelper.NormalizeFinite(quotient, correction) + finalCorrection; } /// Applies the expansion operation without discarding the low component. diff --git a/0-source/Just.PreciseMath/DoubleDouble.cs b/0-source/Just.PreciseMath/DoubleDouble.cs index 23bf0a3..4a955a2 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.cs @@ -98,10 +98,15 @@ public readonly partial struct DoubleDouble : /// Represents a zero value. /// public static DoubleDouble Zero => new(); + + // Mathematical constants store normalized, precomputed binary64 pairs: the + // nearest high and then the nearest residual of the high-precision value. + // No DD arithmetic, parsing, normalization or heap allocation occurs on access. + // Reproduce with the tests' ReferenceData/generate_constants.py. /// /// Represents the ratio of the circumference of a circle to its diameter, specified by the constant, π. /// - public static DoubleDouble PI => new(3.141592653589793, 1.2246467991473532e-16); + public static DoubleDouble Pi => new(3.141592653589793, 1.2246467991473532e-16); /// /// Represents the natural logarithmic base, specified by the constant, e. /// @@ -109,7 +114,85 @@ public readonly partial struct DoubleDouble : /// /// Represents the natural logarithm of value 2. /// - public static DoubleDouble LN2 => new(0.6931471805599453, 2.3190468138462996e-17); + public static DoubleDouble Ln2 => new(0.6931471805599453, 2.3190468138462996e-17); + + /// Gets τ = 2π, the angle of one full turn in radians. + public static DoubleDouble Tau => new(6.283185307179586, 2.4492935982947064e-16); + + /// Gets π/2, the angle of 90 degrees in radians. + public static DoubleDouble PiOver2 => new(1.5707963267948966, 6.123233995736766e-17); + + /// Gets π/3, the angle of 60 degrees in radians. + public static DoubleDouble PiOver3 => new(1.0471975511965979, -1.072081766451091e-16); + + /// Gets π/4, the angle of 45 degrees in radians. + public static DoubleDouble PiOver4 => new(0.7853981633974483, 3.061616997868383e-17); + + /// Gets π/6, the angle of 30 degrees in radians. + public static DoubleDouble PiOver6 => new(0.5235987755982989, -5.360408832255455e-17); + + /// Gets 1/π. + public static DoubleDouble InvPi => new(0.3183098861837907, -1.9678676675182486e-17); + + /// Gets 1/(2π), the factor for converting radians to turns. + public static DoubleDouble InvTau => new(0.15915494309189535, -9.839338337591243e-18); + + /// Gets π/180. Multiply an angle in degrees by this value to obtain radians. + public static DoubleDouble DegToRad => new(0.017453292519943295, 2.9486522708701687e-19); + + /// Gets 180/π. Multiply an angle in radians by this value to obtain degrees. + public static DoubleDouble RadToDeg => new(57.29577951308232, -1.9878495670576283e-15); + + /// Gets 1/e = exp(-1). + public static DoubleDouble InvE => new(0.36787944117144233, -1.2428753672788363e-17); + + /// Gets ln(10), the factor for converting base-10 logarithms to natural logarithms. + public static DoubleDouble Ln10 => new(2.302585092994046, -2.1707562233822494e-16); + + /// Gets log₂(e) = 1/ln(2), the factor for converting natural logarithms to base 2. + public static DoubleDouble Log2E => new(1.4426950408889634, 2.0355273740931033e-17); + + /// Gets log₁₀(e) = 1/ln(10), the factor for converting natural logarithms to base 10. + public static DoubleDouble Log10E => new(0.4342944819032518, 1.098319650216765e-17); + + /// Gets log₂(10), the factor for converting base-10 logarithms to base 2. + public static DoubleDouble Log2Of10 => new(3.321928094887362, 1.661617516973592e-16); + + /// Gets log₁₀(2), the factor for converting base-2 logarithms to base 10. + public static DoubleDouble Log10Of2 => new(0.3010299956639812, -2.8037281277851704e-18); + + /// Gets √2. + public static DoubleDouble Sqrt2 => new(1.4142135623730951, -9.667293313452913e-17); + + /// Gets √3. + public static DoubleDouble Sqrt3 => new(1.7320508075688772, 1.0035084221806903e-16); + + /// Gets √5. + public static DoubleDouble Sqrt5 => new(2.23606797749979, -1.0864230407365012e-16); + + /// Gets 1/√2, also the sine and cosine of π/4. + public static DoubleDouble InvSqrt2 => new(0.7071067811865476, -4.833646656726457e-17); + + /// Gets 1/√3, also the tangent of π/6. + public static DoubleDouble InvSqrt3 => new(0.5773502691896257, 3.3450280739356345e-17); + + /// Gets √π, the Gaussian integral over the real line for exp(-x²). + public static DoubleDouble SqrtPi => new(1.772453850905516, -7.666586499825799e-17); + + /// Gets 1/√π, a Gaussian normalization factor. + public static DoubleDouble InvSqrtPi => new(0.5641895835477563, 7.66772980658294e-18); + + /// Gets 2/√π, the normalization factor in the error-function integral. + public static DoubleDouble TwoInvSqrtPi => new(1.1283791670955126, 1.533545961316588e-17); + + /// Gets √(2π), used in Gaussian integrals and Stirling's approximation. + public static DoubleDouble SqrtTau => new(2.5066282746310007, -1.8328579980459167e-16); + + /// Gets 1/√(2π), the standard normal probability density's normalization factor. + public static DoubleDouble InvSqrtTau => new(0.3989422804014327, -2.49232720227773e-17); + + /// Gets the golden ratio φ = (1 + √5)/2. + public static DoubleDouble GoldenRatio => new(1.618033988749895, -5.432115203682506e-17); #endregion /// diff --git a/0-source/Just.PreciseMath/PreciseMathHelper.cs b/0-source/Just.PreciseMath/PreciseMathHelper.cs index 5532e03..9c5a7b1 100644 --- a/0-source/Just.PreciseMath/PreciseMathHelper.cs +++ b/0-source/Just.PreciseMath/PreciseMathHelper.cs @@ -104,10 +104,35 @@ internal static class PreciseMathHelper return NormalizeFinite(sum, sumError + (middleError + lowError)); } + // Only for DD division's second remainder: normalized finite inputs with highs + // bounded by 2^462, including the correction product's exact fallback results. + // Same sign and normal binade imply Sterbenz-exact high subtraction. Its error + // is +0; with canonical input lows, TwoAdd(+0, lowSum) is (lowSum, +0). + // Eliminate those two transforms, not the low-sum error or final normalization. + [MethodImpl(MethodImplOptions.AggressiveInlining)] + internal static double SubtractDivisionCorrectionHigh(DoubleDouble remainder, DoubleDouble product) + { + ulong remainderBits = BitConverter.DoubleToUInt64Bits(remainder._high); + ulong productBits = BitConverter.DoubleToUInt64Bits(product._high); + if ((remainderBits & 0x7ff0_0000_0000_0000UL) == 0 + || ((remainderBits ^ productBits) & 0xfff0_0000_0000_0000UL) != 0) + { + return (remainder - product)._high; + } + + double high = remainder._high - product._high; + (double low, double lowError) = TwoAdd(remainder._low, -product._low); + (double sum, double sumError) = TwoAdd(high, low); + double error = sumError + lowError; + // Preserve NormalizeFinite's zero-low shortcut, including the high zero's + // sign. The normalized low output is unused by the final quotient correction. + return error == 0.0 ? sum : sum + error; + } + // 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. + // AddFinite and DD division establish 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) diff --git a/1-tests/Just.PreciseMath.Tests/ArithmeticRangeTests.cs b/1-tests/Just.PreciseMath.Tests/ArithmeticRangeTests.cs index 7fedfba..59fe0a2 100644 --- a/1-tests/Just.PreciseMath.Tests/ArithmeticRangeTests.cs +++ b/1-tests/Just.PreciseMath.Tests/ArithmeticRangeTests.cs @@ -80,6 +80,7 @@ public class ArithmeticRangeTests [InlineData("+")] [InlineData("-")] [InlineData("*")] + [InlineData("/")] public void FiniteKernelsPreservePreviousComponentBits(string operation) { // Differential characterization, not an independent accuracy oracle. @@ -113,7 +114,39 @@ public class ArithmeticRangeTests { DoubleDouble expected; DoubleDouble actual; - if (operation == "*") + if (operation == "/") + { + if (left.High == 0.0 || right.High == 0.0 + || Math.Abs(Math.ILogB(left.High)) > 450 || Math.Abs(Math.ILogB(right.High)) > 450) + { + continue; + } + // Retain the original division expression, including its public + // normalization. Sparse corrections may use MultiplyBoundary. + double quotient = left.High / right.High; + DoubleDouble remainder = left - (right * quotient); + double correction = remainder.High / right.High; + double.IsFinite(quotient).ShouldBeTrue(); + double.IsFinite(correction).ShouldBeTrue(); + double.IsFinite(quotient + correction).ShouldBeTrue(); + // Conservative bounds from the entry domain and the first + // scalar product/four-TwoSum remainder, not an accuracy claim. + (Math.Abs(quotient) <= Math.ScaleB(1.0, 901)).ShouldBeTrue(); + (Math.Abs(remainder.High) <= Math.ScaleB(1.0, 457)).ShouldBeTrue(); + (Math.Abs(correction) <= Math.ScaleB(1.0, 908)).ShouldBeTrue(); + (Math.Abs(quotient + correction) < Math.ScaleB(1.0, 909)).ShouldBeTrue(); + DoubleDouble correctionProduct = right * correction; + AssertDivisionResidualCancellation(remainder, correctionProduct); + remainder -= correctionProduct; + double finalCorrection = remainder.High / right.High; + DoubleDouble normalized = DoubleDouble.FromComponents(quotient, correction); + DoubleDouble finiteNormalized = PreciseMathHelper.NormalizeFinite(quotient, correction); + BitConverter.DoubleToInt64Bits(finiteNormalized.High).ShouldBe(BitConverter.DoubleToInt64Bits(normalized.High)); + BitConverter.DoubleToInt64Bits(finiteNormalized.Low).ShouldBe(BitConverter.DoubleToInt64Bits(normalized.Low)); + expected = normalized + finalCorrection; + actual = left / right; + } + else if (operation == "*") { if (left.High == 0.0 || right.High == 0.0) { @@ -161,6 +194,62 @@ public class ArithmeticRangeTests } } + [Fact] + public void DivisionResidualCancellationPreservesHighBitsAtBinadeEdges() + { + // Include both sides of normal/subnormal and binade transitions, with + // canonical zero, dense and sparse lows. Cross-binade/sign/zero cases + // are intentionally ineligible and must retain the general subtraction. + List values = [new(0.0), new(-0.0)]; + foreach (int exponent in new[] { -1074, -1022, -900, -54, 0, 1, 457, 461 }) + { + foreach (double significand in new[] { 1.0, Math.BitIncrement(1.0), 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, -54), -Math.ScaleB(high, -54), + double.Epsilon, -double.Epsilon }) + { + values.Add(DoubleDouble.FromComponents(high, low)); + } + } + } + } + foreach (DoubleDouble remainder in values) + { + foreach (DoubleDouble product in values) + { + AssertDivisionResidualCancellation(remainder, product); + } + } + } + + private static void AssertDivisionResidualCancellation(DoubleDouble remainder, DoubleDouble product) + { + DoubleDouble expected = remainder - product; + double actual = PreciseMathHelper.SubtractDivisionCorrectionHigh(remainder, product); + BitConverter.DoubleToInt64Bits(actual).ShouldBe(BitConverter.DoubleToInt64Bits(expected.High)); + // Independent BCL classification, rather than the proposed exponent-bit + // guard. Within one normal binade, same-sign subtraction is exact. + if (!double.IsNormal(remainder.High) || !double.IsNormal(product.High) + || Math.Sign(remainder.High) != Math.Sign(product.High) + || Math.ILogB(remainder.High) != Math.ILogB(product.High)) + { + return; + } + (double high, double highError) = PreciseMathHelper.TwoAdd(remainder.High, -product.High); + BitConverter.DoubleToInt64Bits(highError).ShouldBe(0L); + (double low, double lowError) = PreciseMathHelper.TwoAdd(remainder.Low, -product.Low); + (double middle, double middleError) = PreciseMathHelper.TwoAdd(highError, low); + BitConverter.DoubleToInt64Bits(middle).ShouldBe(BitConverter.DoubleToInt64Bits(low)); + BitConverter.DoubleToInt64Bits(middleError).ShouldBe(0L); + (double sum, double sumError) = PreciseMathHelper.TwoAdd(high, low); + double error = sumError + lowError; + double simplifiedHigh = error == 0.0 ? sum : sum + error; + BitConverter.DoubleToInt64Bits(simplifiedHigh).ShouldBe(BitConverter.DoubleToInt64Bits(expected.High)); + } + private static void AssertMultiplicationRange(double left, double right) { int exponent = Math.ILogB(left) + Math.ILogB(right); diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs index 104afab..d93e405 100644 --- a/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs @@ -6,6 +6,303 @@ 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() { @@ -232,6 +529,56 @@ public class DoubleDoubleArithmeticTests 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() { @@ -537,6 +884,21 @@ public class DoubleDoubleArithmeticTests } } + 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); diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleFormattingTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleFormattingTests.cs index 914fa27..0efb243 100644 --- a/1-tests/Just.PreciseMath.Tests/DoubleDoubleFormattingTests.cs +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleFormattingTests.cs @@ -58,12 +58,12 @@ public class DoubleDoubleFormattingTests public void FormattingDoesNotRouteThroughDecimalOrRestrictExponentRange() { // Review-1 §6 verified the former decimal-based formatter threw for 1e100, - // printed "0" for 1e-100, threw for NaN/infinity, and emitted for PI digits + // printed "0" for 1e-100, threw for NaN/infinity, and emitted for pi digits // already wrong at binary64 precision. Expected digits are exact references: - // the PI pair rounded to 32 significant digits, ties to even + // the Pi pair rounded to 32 significant digits, ties to even // (Python Fraction/Decimal at precision 120), and the exact binary64 values of // the powers of ten. - DoubleDouble.PI.ToString("G32", CultureInfo.InvariantCulture).ShouldBe("3.1415926535897932384626433832795"); + DoubleDouble.Pi.ToString("G32", CultureInfo.InvariantCulture).ShouldBe("3.1415926535897932384626433832795"); new DoubleDouble(1e100).ToString("G", CultureInfo.InvariantCulture).ShouldBe("1.0000000000000000159028911097599E+100"); new DoubleDouble(1e-100).ToString("G", CultureInfo.InvariantCulture).ShouldBe("1.0000000000000000199918998026029E-100"); // The binary64 values are exactly representable as sums with a zero residual, diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleSpanFormattingTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleSpanFormattingTests.cs index a6ef0eb..4d6ecc0 100644 --- a/1-tests/Just.PreciseMath.Tests/DoubleDoubleSpanFormattingTests.cs +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleSpanFormattingTests.cs @@ -38,7 +38,7 @@ public class DoubleDoubleSpanFormattingTests public void ShortDestinationIsUnchangedAndReportsZero() { char[] buffer = ['!', '!']; - DoubleDouble.PI.TryFormat(buffer, out int written, "G32", CultureInfo.InvariantCulture).ShouldBeFalse(); + DoubleDouble.Pi.TryFormat(buffer, out int written, "G32", CultureInfo.InvariantCulture).ShouldBeFalse(); written.ShouldBe(0); new string(buffer).ShouldBe("!!"); DoubleDouble.One.TryFormat(Span.Empty, out written, provider: CultureInfo.InvariantCulture).ShouldBeFalse(); diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleTests.cs index 004ea1e..9691a30 100644 --- a/1-tests/Just.PreciseMath.Tests/DoubleDoubleTests.cs +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleTests.cs @@ -6,25 +6,81 @@ namespace Just.PreciseMath.Tests; public class DoubleDoubleTests { [Theory] - [InlineData("PI", 3.141592653589793, 1.2246467991473532e-16)] + [InlineData("Pi", 3.141592653589793, 1.2246467991473532e-16)] [InlineData("E", 2.718281828459045, 1.4456468917292502e-16)] - [InlineData("LN2", 0.6931471805599453, 2.3190468138462996e-17)] + [InlineData("Ln2", 0.6931471805599453, 2.3190468138462996e-17)] + [InlineData("Tau", 6.283185307179586, 2.4492935982947064e-16)] + [InlineData("PiOver2", 1.5707963267948966, 6.123233995736766e-17)] + [InlineData("PiOver3", 1.0471975511965979, -1.072081766451091e-16)] + [InlineData("PiOver4", 0.7853981633974483, 3.061616997868383e-17)] + [InlineData("PiOver6", 0.5235987755982989, -5.360408832255455e-17)] + [InlineData("InvPi", 0.3183098861837907, -1.9678676675182486e-17)] + [InlineData("InvTau", 0.15915494309189535, -9.839338337591243e-18)] + [InlineData("DegToRad", 0.017453292519943295, 2.9486522708701687e-19)] + [InlineData("RadToDeg", 57.29577951308232, -1.9878495670576283e-15)] + [InlineData("InvE", 0.36787944117144233, -1.2428753672788363e-17)] + [InlineData("Ln10", 2.302585092994046, -2.1707562233822494e-16)] + [InlineData("Log2E", 1.4426950408889634, 2.0355273740931033e-17)] + [InlineData("Log10E", 0.4342944819032518, 1.098319650216765e-17)] + [InlineData("Log2Of10", 3.321928094887362, 1.661617516973592e-16)] + [InlineData("Log10Of2", 0.3010299956639812, -2.8037281277851704e-18)] + [InlineData("Sqrt2", 1.4142135623730951, -9.667293313452913e-17)] + [InlineData("Sqrt3", 1.7320508075688772, 1.0035084221806903e-16)] + [InlineData("Sqrt5", 2.23606797749979, -1.0864230407365012e-16)] + [InlineData("InvSqrt2", 0.7071067811865476, -4.833646656726457e-17)] + [InlineData("InvSqrt3", 0.5773502691896257, 3.3450280739356345e-17)] + [InlineData("SqrtPi", 1.772453850905516, -7.666586499825799e-17)] + [InlineData("InvSqrtPi", 0.5641895835477563, 7.66772980658294e-18)] + [InlineData("TwoInvSqrtPi", 1.1283791670955126, 1.533545961316588e-17)] + [InlineData("SqrtTau", 2.5066282746310007, -1.8328579980459167e-16)] + [InlineData("InvSqrtTau", 0.3989422804014327, -2.49232720227773e-17)] + [InlineData("GoldenRatio", 1.618033988749895, -5.432115203682506e-17)] public void ConstantsHaveNearestBinary64Residuals(string name, double high, double low) { // Each residual is round_binary64(constant - exact_binary64(high)). - // Reproduced with Python decimal at precision 90: e = Decimal(1).exp(), - // ln(2) = Decimal(2).ln(), pi = 16*atan(1/5) - 4*atan(1/239), using - // atan(x) = sum((-1)^k*x^(2k+1)/(2k+1)) until |term| < 1e-95. - // low = float(reference - Decimal.from_float(float(reference))). + // Reproduce with ReferenceData/generate_constants.py: Python Decimal at + // 160 and 240 digits, Machin's pi, exp/ln/sqrt and exact Fraction splitting. + // The generator records every formula and checks normalized components. DoubleDouble value = name switch { - "PI" => DoubleDouble.PI, + "Pi" => DoubleDouble.Pi, "E" => DoubleDouble.E, - "LN2" => DoubleDouble.LN2, + "Ln2" => DoubleDouble.Ln2, + "Tau" => DoubleDouble.Tau, + "PiOver2" => DoubleDouble.PiOver2, + "PiOver3" => DoubleDouble.PiOver3, + "PiOver4" => DoubleDouble.PiOver4, + "PiOver6" => DoubleDouble.PiOver6, + "InvPi" => DoubleDouble.InvPi, + "InvTau" => DoubleDouble.InvTau, + "DegToRad" => DoubleDouble.DegToRad, + "RadToDeg" => DoubleDouble.RadToDeg, + "InvE" => DoubleDouble.InvE, + "Ln10" => DoubleDouble.Ln10, + "Log2E" => DoubleDouble.Log2E, + "Log10E" => DoubleDouble.Log10E, + "Log2Of10" => DoubleDouble.Log2Of10, + "Log10Of2" => DoubleDouble.Log10Of2, + "Sqrt2" => DoubleDouble.Sqrt2, + "Sqrt3" => DoubleDouble.Sqrt3, + "Sqrt5" => DoubleDouble.Sqrt5, + "InvSqrt2" => DoubleDouble.InvSqrt2, + "InvSqrt3" => DoubleDouble.InvSqrt3, + "SqrtPi" => DoubleDouble.SqrtPi, + "InvSqrtPi" => DoubleDouble.InvSqrtPi, + "TwoInvSqrtPi" => DoubleDouble.TwoInvSqrtPi, + "SqrtTau" => DoubleDouble.SqrtTau, + "InvSqrtTau" => DoubleDouble.InvSqrtTau, + "GoldenRatio" => DoubleDouble.GoldenRatio, _ => throw new ArgumentOutOfRangeException(nameof(name)), }; value.High.ShouldBe(high); value.Low.ShouldBe(low); + BitConverter.DoubleToInt64Bits(value.High).ShouldBe(BitConverter.DoubleToInt64Bits(high)); + BitConverter.DoubleToInt64Bits(value.Low).ShouldBe(BitConverter.DoubleToInt64Bits(low)); + DoubleDouble.IsCanonical(value).ShouldBeTrue(); + DoubleDouble.IsFinite(value).ShouldBeTrue(); + value.Low.ShouldNotBe(0.0); } [Fact] @@ -46,7 +102,7 @@ public class DoubleDoubleTests [InlineData(3.141592653589793, 1.2246467991473532e-16)] [InlineData(2.718281828459045, 1.4456468917292502e-16)] [InlineData(0.6931471805599453, 2.3190468138462996e-17)] - public void AdditiveIdentity(double high, double low) + public void AdditiveIdentityAdd(double high, double low) { DoubleDouble value = new(high, low); @@ -70,7 +126,32 @@ public class DoubleDoubleTests [InlineData(3.141592653589793, 1.2246467991473532e-16)] [InlineData(2.718281828459045, 1.4456468917292502e-16)] [InlineData(0.6931471805599453, 2.3190468138462996e-17)] - public void MultiplicativeIdentity(double high, double low) + public void AdditiveIdentitySubtract(double high, double low) + { + DoubleDouble value = new(high, low); + DoubleDouble negativeValue = -value; + + DoubleDouble result = value - DoubleDouble.AdditiveIdentity; + DoubleDouble resultInversedOrder = DoubleDouble.AdditiveIdentity - value; + + result.High.ShouldBe(high); + result.Low.ShouldBe(low); + + resultInversedOrder.High.ShouldBe(negativeValue.High); + resultInversedOrder.Low.ShouldBe(negativeValue.Low); + } + + [Theory] + [InlineData(1.0, 0.0)] + [InlineData(10.0, 0.0)] + [InlineData(100.0, 0.0)] + [InlineData(-1.0, 0.0)] + [InlineData(-10.0, 0.0)] + [InlineData(-100.0, 0.0)] + [InlineData(3.141592653589793, 1.2246467991473532e-16)] + [InlineData(2.718281828459045, 1.4456468917292502e-16)] + [InlineData(0.6931471805599453, 2.3190468138462996e-17)] + public void MultiplicativeIdentityMultiply(double high, double low) { DoubleDouble value = new(high, low); @@ -83,4 +164,24 @@ public class DoubleDoubleTests resultInversedOrder.High.ShouldBe(high); resultInversedOrder.Low.ShouldBe(low); } + + [Theory] + [InlineData(1.0, 0.0)] + [InlineData(10.0, 0.0)] + [InlineData(100.0, 0.0)] + [InlineData(-1.0, 0.0)] + [InlineData(-10.0, 0.0)] + [InlineData(-100.0, 0.0)] + [InlineData(3.141592653589793, 1.2246467991473532e-16)] + [InlineData(2.718281828459045, 1.4456468917292502e-16)] + [InlineData(0.6931471805599453, 2.3190468138462996e-17)] + public void MultiplicativeIdentityDivide(double high, double low) + { + DoubleDouble value = new(high, low); + + DoubleDouble result = value / DoubleDouble.MultiplicativeIdentity; + + result.High.ShouldBe(high); + result.Low.ShouldBe(low); + } } diff --git a/1-tests/Just.PreciseMath.Tests/ReferenceData/generate_constants.py b/1-tests/Just.PreciseMath.Tests/ReferenceData/generate_constants.py new file mode 100644 index 0000000..f514a21 --- /dev/null +++ b/1-tests/Just.PreciseMath.Tests/ReferenceData/generate_constants.py @@ -0,0 +1,76 @@ +"""Reproduce DoubleDouble's precomputed mathematical constants. + +Run: python3 -B 1-tests/Just.PreciseMath.Tests/ReferenceData/generate_constants.py +Standard library only; no dependency on the implementation under test. +Decimal evaluates formulas at 160 and 240 digits; Fraction-based splitting +rounds the high and then the exact residual to nearest-even binary64. +""" + +from decimal import Decimal, localcontext +from fractions import Fraction + +from generate_irrational_arithmetic import arctan_inverse, split + + +def constants(precision: int) -> dict[str, tuple[float, float]]: + with localcontext() as context: + context.prec = precision + one = Decimal(1) + pi = 16 * arctan_inverse(5) - 4 * arctan_inverse(239) + e = one.exp() + ln2 = Decimal(2).ln() + ln10 = Decimal(10).ln() + sqrt2 = Decimal(2).sqrt() + sqrt3 = Decimal(3).sqrt() + sqrt5 = Decimal(5).sqrt() + sqrt_pi = pi.sqrt() + sqrt_tau = (2 * pi).sqrt() + values = { + "Pi": pi, + "E": e, + "Ln2": ln2, + "Tau": 2 * pi, + "PiOver2": pi / 2, + "PiOver3": pi / 3, + "PiOver4": pi / 4, + "PiOver6": pi / 6, + "InvPi": one / pi, + "InvTau": one / (2 * pi), + "DegToRad": pi / 180, + "RadToDeg": 180 / pi, + "InvE": one / e, + "Ln10": ln10, + "Log2E": one / ln2, + "Log10E": one / ln10, + "Log2Of10": ln10 / ln2, + "Log10Of2": ln2 / ln10, + "Sqrt2": sqrt2, + "Sqrt3": sqrt3, + "Sqrt5": sqrt5, + "InvSqrt2": one / sqrt2, + "InvSqrt3": one / sqrt3, + "SqrtPi": sqrt_pi, + "InvSqrtPi": one / sqrt_pi, + "TwoInvSqrtPi": 2 / sqrt_pi, + "SqrtTau": sqrt_tau, + "InvSqrtTau": one / sqrt_tau, + "GoldenRatio": (one + sqrt5) / 2, + } + pairs = {name: split(Fraction(value)) for name, value in values.items()} + for name, value in values.items(): + high, low = pairs[name] + assert high + low == high, f"Unnormalized constant: {name}" + assert low != 0.0, f"Lost residual: {name}" + assert abs(Fraction(high) + Fraction(low) - Fraction(value)) <= abs(Fraction(value)) / (1 << 105) + return pairs + + +if __name__ == "__main__": + pairs = constants(160) + assert pairs == constants(240), "Increase precision: binary64 constants did not stabilize" + print("// Expected test components") + for name, (high, low) in pairs.items(): + print(f' [InlineData("{name}", {high!r}, {low!r})]') + print("\n// Precomputed properties") + for name, (high, low) in pairs.items(): + print(f" public static DoubleDouble {name} => new({high!r}, {low!r});") diff --git a/1-tests/Just.PreciseMath.Tests/ReferenceData/generate_irrational_arithmetic.py b/1-tests/Just.PreciseMath.Tests/ReferenceData/generate_irrational_arithmetic.py new file mode 100644 index 0000000..663aaa3 --- /dev/null +++ b/1-tests/Just.PreciseMath.Tests/ReferenceData/generate_irrational_arithmetic.py @@ -0,0 +1,77 @@ +"""Print BasicArithmeticMatchesPrecomputedIrrationalResults InlineData. + +Run with Python 3.11+; standard library only. Decimal supplies irrational +constants; Fraction supplies exact arithmetic on their stored binary64 pairs. +The test suite consumes the literals, not this script or Python at runtime. +""" + +from decimal import Decimal, localcontext +from fractions import Fraction + + +def arctan_inverse(inverse: int) -> Decimal: + """Evaluate atan(1/inverse) using its alternating power series.""" + x = Decimal(1) / inverse + power = x + total = x + index = 1 + while True: + power *= -(x * x) + updated = total + power / (2 * index + 1) + if updated == total: + return total + total = updated + index += 1 + + +def split(value: Fraction) -> tuple[float, float]: + """Round high, then the exact residual, to nearest-even binary64.""" + high = float(value) + return high, float(value - Fraction(high)) + + +def fixtures(precision: int) -> list[tuple[str, list[float]]]: + with localcontext() as context: + context.prec = precision + constants = { + "pi": 16 * arctan_inverse(5) - 4 * arctan_inverse(239), + "e": Decimal(1).exp(), + "sqrt2": Decimal(2).sqrt(), + "sqrt3": Decimal(3).sqrt(), + "ln2": Decimal(2).ln(), + "phi": (1 + Decimal(5).sqrt()) / 2, + } + pairs = {name: split(Fraction(value)) for name, value in constants.items()} + + rows = [] + for left_name, right_name in [("pi", "e"), ("sqrt2", "sqrt3"), ("ln2", "phi")]: + for overload in ["DD/DD", "DD/double", "double/DD"]: + left_high, left_low = pairs[left_name] + right_high, right_low = pairs[right_name] + if overload == "DD/double": + right_low = 0.0 + if overload == "double/DD": + left_low = 0.0 + left = Fraction(left_high) + Fraction(left_low) + right = Fraction(right_high) + Fraction(right_low) + values = [left_high, left_low, right_high, right_low] + for expected in [left + right, left - right, left * right, left / right]: + high, low = split(expected) + # These precomputed references are much closer than the test's + # 2^-100 bound; this does not assert library correct rounding. + assert abs(Fraction(high) + Fraction(low) - expected) <= abs(expected) / (1 << 105) + values.extend([high, low]) + rows.append((f"{left_name}, {right_name}: {overload}", values)) + return rows + + +if __name__ == "__main__": + rows = fixtures(160) + assert rows == fixtures(240), "Increase precision: binary64 fixtures did not stabilize" + for label, values in rows: + print(f" // {label}") + overload = label.split(": ")[1] + print(f' [InlineData("{overload}",') + for index in range(0, len(values), 2): + suffix = ")]" if index == len(values) - 2 else "," + print(f" {values[index]!r}, {values[index + 1]!r}{suffix}") diff --git a/README.md b/README.md index 05f9038..3fcb474 100644 --- a/README.md +++ b/README.md @@ -16,8 +16,9 @@ for a single `double`, or `DoubleDouble.FromComponents(high, low)` for arbitrary components. The factory normalizes finite sums and canonicalizes NaN/infinity with a positive-zero low component. The two-component constructor is internal and performs no normalization or validation; it is reserved for trusted, -already-normalized results. The constants `PI`, `E`, and `LN2` include binary64 -residuals checked against independently computed high-precision values. +already-normalized results. Mathematical constants use precomputed high/low pairs +checked against independently computed high-precision values; accessing them does +not perform double-double arithmetic or allocate on the heap. - Arithmetic: unary `+`/`-`, binary `+`, `-`, `*`, `/`, and both operand orders with a `double`. Addition retains residuals under cancellation; multiplication uses @@ -44,6 +45,39 @@ with exact component checks for selected representable cases. Near underflow, extended precision necessarily decreases; overflow produces infinity. Performance of the allocating exponent-boundary path is not covered by the basic benchmarks. +## Predefined mathematical constants + +All constants below are static `DoubleDouble` properties. Each stores the nearest +binary64 high component followed by the nearest binary64 residual, rather than +calculating a ratio, root, or logarithm on access. Names use PascalCase, including +`Pi`, `E`, and `Ln2`. + +| Group | Properties and values | +|---|---| +| Circle and common angles | `Pi` (π), `Tau` (2π), `PiOver2`, `PiOver3`, `PiOver4`, `PiOver6` | +| Angular conversion | `DegToRad` (π/180), `RadToDeg` (180/π), `InvPi` (1/π), `InvTau` (1/(2π), radians to turns) | +| Exponential and logarithmic | `E`, `InvE` (1/e), `Ln2` (ln 2), `Ln10` (ln 10) | +| Log-base conversion | `Log2E` (1/ln 2), `Log10E` (1/ln 10), `Log2Of10` (ln 10/ln 2), `Log10Of2` (ln 2/ln 10) | +| Roots and geometry | `Sqrt2`, `Sqrt3`, `Sqrt5`, `InvSqrt2`, `InvSqrt3`, `GoldenRatio` ((1+√5)/2) | +| Gaussian and error-function factors | `SqrtPi`, `InvSqrtPi`, `TwoInvSqrtPi` (2/√π), `SqrtTau` (√(2π)), `InvSqrtTau` (1/√(2π)) | + +Multiply by conversion factors instead of recomputing them: + +```csharp +using Just.PreciseMath; + +DoubleDouble degrees = new(180.0); +DoubleDouble radians = degrees * DoubleDouble.DegToRad; +DoubleDouble convertedDegrees = radians * DoubleDouble.RadToDeg; +DoubleDouble quarterTurn = DoubleDouble.PiOver2; +``` + +The factors avoid deriving constants at runtime; the multiplication itself remains +approximate DD arithmetic, so conversions are not guaranteed exact round trips. +The list is mathematical and dimensionless, not a table of unit-dependent physical +constants. Precomputed roots/logarithms do not imply general `Sqrt`/`Log` functions +are implemented. + ## Conversions and formatting - Explicit conversions support `double`, `float`, `int`, `long`, and `decimal`