From beb9a084d088d484cc3908896aeefdac9e3b66b2 Mon Sep 17 00:00:00 2001 From: just Date: Mon, 14 Sep 2026 22:25:50 +0400 Subject: [PATCH] initial DDMath implementation --- 0-source/Just.PreciseMath/DDMath.cs | 57 +++++++ .../PreciseMathSqrtTests.cs | 161 ++++++++++++++++++ .../PreciseMathTests.cs | 35 ++++ README.md | 32 +++- 4 files changed, 281 insertions(+), 4 deletions(-) create mode 100644 0-source/Just.PreciseMath/DDMath.cs create mode 100644 1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs create mode 100644 1-tests/Just.PreciseMath.Tests/PreciseMathTests.cs diff --git a/0-source/Just.PreciseMath/DDMath.cs b/0-source/Just.PreciseMath/DDMath.cs new file mode 100644 index 0000000..99b3b7b --- /dev/null +++ b/0-source/Just.PreciseMath/DDMath.cs @@ -0,0 +1,57 @@ +namespace Just.PreciseMath; + +/// Provides mathematical functions for normalized double-double values. +/// Results use fixed-size double-double precision, not arbitrary precision. +public static class DDMath +{ + /// Returns the absolute value without discarding the low component. + /// Both signs of zero become positive zero; NaN is canonicalized and either infinity becomes positive infinity. + [Pure] + public static DoubleDouble Abs(DoubleDouble value) + { + return DoubleDouble.Abs(value); + } + + /// Returns the nonnegative square root with double-double precision. + /// + /// Uses power-of-two scaling and an FMA-based Newton correction. Results are + /// approximate, not guaranteed correctly rounded; finite positive inputs are + /// tested against a relative error bound of 2^-100 across the binary64 range. + /// Signed zero and positive infinity are preserved. Negative nonzero values + /// (including negative infinity) and NaN return canonical NaN. + /// + [Pure] + public static DoubleDouble Sqrt(DoubleDouble value) + { + if (double.IsNaN(value.High) || value.High < 0.0) + { + return DoubleDouble.NaN; + } + if (value.High == 0.0 || double.IsPositiveInfinity(value.High)) + { + return value; + } + + // Choose an even exponent, rounding negative odd exponents down too. + // For positive finite normalized inputs, exponent is in [-1074, 1022] + // and the scaled high is in [1, 4). Scaling the high is exact; any + // underflow in a very sparse low is far below the relative error bound. + int exponent = Math.ILogB(value.High) & ~1; + double high = Math.ScaleB(value.High, -exponent); + double low = Math.ScaleB(value.Low, -exponent); + double estimate = Math.Sqrt(high); + + // Estimate is in [1, 2]. FMA avoids rounding its square before + // cancellation. One Newton correction reduces O(u) error to O(u^2), + // where u=2^-53; both components of the input contribute to the residual. + double residual = Math.FusedMultiplyAdd(-estimate, estimate, high) + low; + double correction = residual / (2.0 * estimate); + DoubleDouble root = DoubleDouble.FromComponents(estimate, correction); + + // The root high stays normal, even if rounding reaches the next binade. + // Normalize again to canonicalize a zero low, including any underflow + // of an exceptionally sparse correction during rescaling. + int rootExponent = exponent / 2; + return DoubleDouble.FromComponents(Math.ScaleB(root.High, rootExponent), Math.ScaleB(root.Low, rootExponent)); + } +} diff --git a/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs b/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs new file mode 100644 index 0000000..a9cf5c4 --- /dev/null +++ b/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs @@ -0,0 +1,161 @@ +using System.Numerics; +using Shouldly; +using Xunit; + +namespace Just.PreciseMath.Tests; + +public class PreciseMathSqrtTests +{ + [Theory] + [InlineData(1.0, 1.0)] + [InlineData(4.0, 2.0)] + [InlineData(9.0, 3.0)] + [InlineData(0.25, 0.5)] + [InlineData(2.25, 1.5)] + public void ExactBinarySquaresHaveExactRoots(double square, double root) + { + DoubleDouble actual = DDMath.Sqrt(new DoubleDouble(square)); + actual.High.ShouldBe(root); + BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L); + } + + [Fact] + public void PowersOfFourHaveExactRootsAcrossTheFiniteRange() + { + // sqrt(2^(2k)) = 2^k exactly, including the minimum subnormal input. + for (int exponent = -1074; exponent <= 1022; exponent += 2) + { + DoubleDouble actual = DDMath.Sqrt(new DoubleDouble(Math.ScaleB(1.0, exponent))); + actual.High.ShouldBe(Math.ScaleB(1.0, exponent / 2)); + BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L); + } + } + + [Theory] + [InlineData(2.0)] + [InlineData(3.0)] + [InlineData(5.0)] + [InlineData(1e-308)] + [InlineData(1e308)] + public void IrrationalRootsRetainMoreThanBinary64Precision(double value) + { + DoubleDouble input = new(value); + DoubleDouble actual = DDMath.Sqrt(input); + actual.Low.ShouldNotBe(0.0); + AssertSqrtBound(input, actual); + } + + [Fact] + public void LowComponentsOnEitherSideOfAnExactSquareAffectTheRoot() + { + foreach (double low in new[] { Math.ScaleB(1.0, -54), -Math.ScaleB(1.0, -54) }) + { + DoubleDouble input = DoubleDouble.FromComponents(1.0, low); + DoubleDouble actual = DDMath.Sqrt(input); + actual.High.ShouldBe(1.0); + (Math.Sign(actual.Low) == Math.Sign(low)).ShouldBeTrue(); + AssertSqrtBound(input, actual); + } + } + + [Fact] + public void SpecialValuesMatchBinary64AndRemainCanonical() + { + double[] values = [0.0, -0.0, double.PositiveInfinity, double.NegativeInfinity, + double.NaN, -1.0, -double.Epsilon, double.MinValue]; + foreach (double value in values) + { + DoubleDouble actual = DDMath.Sqrt(new DoubleDouble(value)); + double expected = Math.Sqrt(value); + // Construction canonicalizes all NaNs to double.NaN, including + // domain errors; zeros must instead be compared by their sign bits. + BitConverter.DoubleToInt64Bits(actual.High).ShouldBe( + BitConverter.DoubleToInt64Bits(double.IsNaN(expected) ? double.NaN : expected)); + BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L); + DoubleDouble.IsCanonical(actual).ShouldBeTrue(); + } + + DoubleDouble negative = DoubleDouble.FromComponents(-1.0, Math.ScaleB(1.0, -54)); + DoubleDouble.IsNaN(DDMath.Sqrt(negative)).ShouldBeTrue(); + } + + [Fact] + public void PositiveFiniteInputsMeetExactRelativeErrorBound() + { + // Every binary64 exponent, including subnormal input exponents. Test zero, + // dense, and sparse lows of either sign, without relying on DD arithmetic + // to compute the expected square root or its square. + Random random = new(314159); + for (int exponent = -1074; exponent <= 1023; ++exponent) + { + double high = Math.ScaleB(1.0 + random.NextDouble(), exponent); + double low = Math.ScaleB(random.NextDouble(), exponent - 54); + double power = Math.ScaleB(1.0, exponent); + foreach (double boundary in new[] { Math.BitDecrement(power), power, Math.BitIncrement(power) }) + { + if (boundary > 0.0) + { + DoubleDouble input = new(boundary); + AssertSqrtBound(input, DDMath.Sqrt(input)); + } + } + foreach (double residual in new[] { 0.0, low, -low, double.Epsilon, -double.Epsilon }) + { + DoubleDouble input = DoubleDouble.FromComponents(high, residual); + if (input.High > 0.0) + { + AssertSqrtBound(input, DDMath.Sqrt(input)); + } + } + } + + // Pin the tiny-input precision loss described in the optional legacy review, + // and guard against intermediate squaring overflow near the finite limit. + DoubleDouble[] boundaries = [new(1e-308), new(double.Epsilon), + new(Math.BitDecrement(Math.ScaleB(1.0, -1022))), new(Math.ScaleB(1.0, -1022)), + new(Math.BitIncrement(Math.ScaleB(1.0, -1022))), new(double.MaxValue), + DoubleDouble.FromComponents(double.MaxValue, Math.ScaleB(1.0, 969)), + DoubleDouble.FromComponents(double.MaxValue, -Math.ScaleB(1.0, 969)), + DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -53)), + DoubleDouble.FromComponents(1.0, -Math.ScaleB(1.0, -54)), + DoubleDouble.FromComponents(4.0, Math.ScaleB(1.0, -51)), + DoubleDouble.FromComponents(4.0, -Math.ScaleB(1.0, -52))]; + foreach (DoubleDouble input in boundaries) + { + AssertSqrtBound(input, DDMath.Sqrt(input)); + } + } + + private static void AssertSqrtBound(DoubleDouble input, DoubleDouble actual) + { + DoubleDouble.IsFinite(actual).ShouldBeTrue(); + DoubleDouble.IsCanonical(actual).ShouldBeTrue(); + (actual.High > 0.0).ShouldBeTrue(); + + // Independent exact rational oracle: x = X*2^-1074, y = Y*2^-1074. + // For positive x and y, |y/sqrt(x) - 1| <= t iff + // x*(1-t)^2 <= y^2 <= x*(1+t)^2. With t=2^-100 we + // cross-multiply to integers: no rounded sqrt, product, or DD conversion. + BigInteger x = (Units(input.High) + Units(input.Low)) << 1074; + BigInteger y = Units(actual.High) + Units(actual.Low); + BigInteger scale = BigInteger.One << 100; + BigInteger squared = (y * y) << 200; + BigInteger lower = x * (scale - 1) * (scale - 1); + BigInteger upper = x * (scale + 1) * (scale + 1); + (squared >= lower && squared <= upper).ShouldBeTrue( + $"Sqrt bound failed for ({input.High:R}, {input.Low:R}): ({actual.High:R}, {actual.Low:R})"); + } + + 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; + } +} diff --git a/1-tests/Just.PreciseMath.Tests/PreciseMathTests.cs b/1-tests/Just.PreciseMath.Tests/PreciseMathTests.cs new file mode 100644 index 0000000..66ca109 --- /dev/null +++ b/1-tests/Just.PreciseMath.Tests/PreciseMathTests.cs @@ -0,0 +1,35 @@ +using Shouldly; +using Xunit; + +namespace Just.PreciseMath.Tests; + +public class PreciseMathTests +{ + [Fact] + public void AbsPreservesBothComponentsAndCanonicalSpecialValues() + { + double[] scalars = [0.0, -0.0, 1.0, -1.0, double.Epsilon, -double.Epsilon, + double.MaxValue, double.MinValue, double.PositiveInfinity, double.NegativeInfinity, double.NaN]; + foreach (double scalar in scalars) + { + // Binary64 is an independent reference when the input has no residual. + DoubleDouble actual = DDMath.Abs(new DoubleDouble(scalar)); + AssertBits(actual, double.IsNaN(scalar) ? double.NaN : Math.Abs(scalar), 0.0); + DoubleDouble.IsCanonical(actual).ShouldBeTrue(); + } + + foreach (double low in new[] { Math.ScaleB(1.0, -54), -Math.ScaleB(1.0, -54), double.Epsilon, -double.Epsilon }) + { + DoubleDouble positive = DoubleDouble.FromComponents(1.0, low); + DoubleDouble negative = DoubleDouble.FromComponents(-1.0, -low); + AssertBits(DDMath.Abs(positive), 1.0, low); + AssertBits(DDMath.Abs(negative), 1.0, low); + } + } + + private static void AssertBits(DoubleDouble actual, double high, double low) + { + BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(BitConverter.DoubleToInt64Bits(high)); + BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(BitConverter.DoubleToInt64Bits(low)); + } +} diff --git a/README.md b/README.md index 3fcb474..7df7270 100644 --- a/README.md +++ b/README.md @@ -75,8 +75,32 @@ 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. +constants. Precomputed logarithms do not imply a general `Log` function is implemented. + +## Mathematical functions + +The initial `DDMath` static class provides: + +- `Abs(DoubleDouble)`: preserves both components, maps either signed zero to + positive zero and either infinity to positive infinity, and returns canonical NaN. + It shares the existing `DoubleDouble.Abs` implementation. +- `Sqrt(DoubleDouble)`: uses power-of-two scaling and an FMA-based Newton correction + to retain extended precision, including for subnormal inputs, without squaring + an unscaled estimate near the exponent limits. Signed zero and positive infinity + are preserved; negative nonzero inputs and NaN return canonical NaN. + +```csharp +using Just.PreciseMath; + +DoubleDouble root = DDMath.Sqrt(new DoubleDouble(2.0)); +DoubleDouble magnitude = DDMath.Abs(-root); +``` + +Square-root tests compare the exact component sum against a `2^-100` relative +error bound using integer inequalities. They include samples at every binary64 +exponent, boundary neighbors, both signs of the low component, and exact binary +squares. This is a tested approximate-accuracy contract, not exhaustive coverage +of all component pairs or a guarantee of correctly rounded results. ## Conversions and formatting @@ -165,8 +189,8 @@ bool success = DoubleDouble.TryParse("1.25e-2".AsSpan(), CultureInfo.InvariantCu ## Deferred scope -The planned `PreciseMath` static class and its `Abs`, `Sqrt`, `Pow`, `Exp`, and -`Log` functions are not implemented. Generic-math interfaces beyond `ISignedNumber`, +The `DDMath.Pow`, `DDMath.Exp`, and `DDMath.Log` functions remain unimplemented. +Generic-math interfaces beyond `ISignedNumber`, additional text formats/general round-trip formatting, and non-arithmetic performance benchmarks remain deferred. Replacing allocating arithmetic boundary fallbacks is also deferred; the current