initial DDMath implementation

0-source/Just.PreciseMath/DDMath.cs

@@ -0,0 +1,57 @@
+namespace Just.PreciseMath;
+
+/// <summary>Provides mathematical functions for normalized double-double values.</summary>
+/// <remarks>Results use fixed-size double-double precision, not arbitrary precision.</remarks>
+public static class DDMath
+{
+    /// <summary>Returns the absolute value without discarding the low component.</summary>
+    /// <remarks>Both signs of zero become positive zero; NaN is canonicalized and either infinity becomes positive infinity.</remarks>
+    [Pure]
+    public static DoubleDouble Abs(DoubleDouble value)
+    {
+        return DoubleDouble.Abs(value);
+    }
+
+    /// <summary>Returns the nonnegative square root with double-double precision.</summary>
+    /// <remarks>
+    /// 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.
+    /// </remarks>
+    [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));
+    }
+}

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;
+    }
+}

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));
+    }
+}

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
+constants. Precomputed logarithms do not imply a general `Log` function is implemented.
-are 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
+The `DDMath.Pow`, `DDMath.Exp`, and `DDMath.Log` functions remain unimplemented.
-`Log` functions are not implemented. Generic-math interfaces beyond `ISignedNumber`,
+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