inv sqrt implementation

0-source/Just.PreciseMath/DDMath.cs

@@ -55,6 +55,59 @@ public static partial class DDMath
         return DoubleDouble.FromComponents(Math.ScaleB(root.High, rootExponent), Math.ScaleB(root.Low, rootExponent));
     }
 
+    /// <summary>Returns the reciprocal square root with double-double precision.</summary>
+    /// <remarks>
+    /// Uses power-of-two scaling and a compensated Newton step without double-double
+    /// division. Results are approximate, not guaranteed correctly rounded; positive
+    /// finite inputs are tested against a relative error bound of 2^-100 across the
+    /// binary64 range. Positive/negative zero maps to positive/negative infinity,
+    /// positive infinity to positive zero, and negative nonzero values or NaN to
+    /// canonical NaN.
+    /// </remarks>
+    [Pure]
+    public static DoubleDouble InvSqrt(DoubleDouble value)
+    {
+        if (double.IsNaN(value.High) || value.High < 0.0)
+        {
+            return DoubleDouble.NaN;
+        }
+        if (value.High == 0.0)
+        {
+            return new DoubleDouble(Math.CopySign(double.PositiveInfinity, value.High));
+        }
+        if (double.IsPositiveInfinity(value.High))
+        {
+            return DoubleDouble.Zero;
+        }
+
+        // As in Sqrt, the even exponent is in [-1074, 1022], high is in [1, 4),
+        // and any underflow in a sparse scaled low is below the 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 = 1.0 / Math.Sqrt(high);
+
+        // With u=2^-53, the seed has O(u) relative error. Preserve the square's
+        // FMA residual before cancellation in 1 - (high + low)*estimate^2:
+        // using only the rounded square would leave O(u) error after refinement.
+        // All leading products are bounded and normal. Omitting low*squareError
+        // contributes only O(u^2), as does one Newton step's remaining error.
+        double square = estimate * estimate;
+        double squareError = Math.FusedMultiplyAdd(estimate, estimate, -square);
+        double residual = Math.FusedMultiplyAdd(-high, square, 1.0);
+        residual = Math.FusedMultiplyAdd(-high, squareError, residual);
+        residual = Math.FusedMultiplyAdd(-low, square, residual);
+        double correction = (0.5 * estimate) * residual;
+        DoubleDouble inverseRoot = DoubleDouble.FromComponents(estimate, correction);
+
+        // The result high is normal and finite over the entire positive input
+        // range. Rescale in the opposite direction to Sqrt and canonicalize any
+        // zero low, including underflow of an exceptionally sparse correction.
+        int inverseExponent = -(exponent / 2);
+        return DoubleDouble.FromComponents(Math.ScaleB(inverseRoot.High, inverseExponent),
+            Math.ScaleB(inverseRoot.Low, inverseExponent));
+    }
+
     // For finite, nonzero normalized significands and exponents bounded by the
     // integer-power domain (|exponent| < 2^42). Keep the exponent separate until
     // the final result so an intermediate cannot overflow before reciprocation.

1-tests/Just.PreciseMath.Tests/PreciseMathInvSqrtTests.cs

@@ -0,0 +1,186 @@
+using System.Numerics;
+using Shouldly;
+using Xunit;
+
+namespace Just.PreciseMath.Tests;
+
+public class PreciseMathInvSqrtTests
+{
+    [Theory]
+    [InlineData(2.0)]
+    [InlineData(3.0)]
+    [InlineData(5.0)]
+    [InlineData(1e-308)]
+    [InlineData(1e308)]
+    public void IrrationalInverseRootsRetainMoreThanBinary64Precision(double value)
+    {
+        DoubleDouble input = new(value);
+        DoubleDouble actual = DDMath.InvSqrt(input);
+        actual.Low.ShouldNotBe(0.0);
+        AssertInvSqrtBound(input, actual);
+    }
+
+    [Fact]
+    public void SpecialValuesMatchReciprocalSquareRootAndRemainCanonical()
+    {
+        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.InvSqrt(new DoubleDouble(value));
+            // Binary64 is an independent oracle only for these special values.
+            double expected = 1.0 / Math.Sqrt(value);
+            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.InvSqrt(negative)).ShouldBeTrue();
+    }
+
+    [Fact]
+    public void PowersOfFourHaveExactInverseRootsAcrossTheFiniteRange()
+    {
+        // 1/sqrt(2^(2k)) = 2^-k exactly, including the minimum subnormal input.
+        for (int exponent = -1074; exponent <= 1022; exponent += 2)
+        {
+            DoubleDouble actual = DDMath.InvSqrt(new DoubleDouble(Math.ScaleB(1.0, exponent)));
+            actual.High.ShouldBe(Math.ScaleB(1.0, -(exponent / 2)));
+            BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L);
+        }
+    }
+
+    [Fact]
+    public void LowComponentsOnEitherSideOfOneAffectTheInverseRoot()
+    {
+        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.InvSqrt(input);
+            actual.High.ShouldBe(1.0);
+            Math.Sign(actual.Low).ShouldBe(-Math.Sign(low));
+            AssertInvSqrtBound(input, actual);
+        }
+    }
+
+    [Fact]
+    public void PositiveFiniteInputsMeetExactRelativeErrorBound()
+    {
+        // Every input exponent, both parities, dense/sparse lows of either sign,
+        // and neighbors of binade transitions. The oracle uses the exact stored
+        // input sum, including any rounding during public input construction.
+        Random random = new(271828);
+        for (int exponent = -1074; exponent <= 1023; ++exponent)
+        {
+            double high = Math.ScaleB(1.0 + random.NextDouble(), exponent);
+            double low = Math.ScaleB(random.NextDouble(), exponent - 53);
+            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);
+                    AssertInvSqrtBound(input, DDMath.InvSqrt(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)
+                {
+                    AssertInvSqrtBound(input, DDMath.InvSqrt(input));
+                }
+            }
+        }
+
+        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.BitDecrement(Math.ScaleB(1.0, 970))),
+            DoubleDouble.FromComponents(double.MaxValue, -Math.ScaleB(1.0, 970)),
+            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)),
+            DoubleDouble.FromComponents(0.0, double.Epsilon)];
+        foreach (DoubleDouble input in boundaries)
+        {
+            AssertInvSqrtBound(input, DDMath.InvSqrt(input));
+        }
+    }
+
+    [Fact]
+    public void DenseLowNormalizationBoundariesMeetTheBound()
+    {
+        // Exercise adjacent lows at half an ulp, both signs, for each exponent.
+        Random random = new(161803);
+        for (int exponent = -1074; exponent <= 1023; ++exponent)
+        {
+            double high = Math.ScaleB(1.0 + random.NextDouble(), exponent);
+            double halfUlp = Math.ScaleB(1.0, exponent - 53);
+            foreach (double magnitude in new[] { Math.BitDecrement(halfUlp), halfUlp, Math.BitIncrement(halfUlp) })
+            {
+                foreach (double low in new[] { magnitude, -magnitude })
+                {
+                    DoubleDouble input = DoubleDouble.FromComponents(high, low);
+                    if (DoubleDouble.IsFinite(input) && input.High > 0.0)
+                    {
+                        AssertInvSqrtBound(input, DDMath.InvSqrt(input));
+                    }
+                }
+            }
+        }
+    }
+
+    [Fact]
+    public void RepresentableSparseCorrectionsAreNotDiscarded()
+    {
+        // (1+d)^(-1/2) = 1-d/2+O(d^2). For these dyadic d, the quadratic
+        // term is below half an ulp of d/2, even at the minimum subnormal.
+        foreach (int exponent in new[] { -100, -500, -1000, -1073 })
+        {
+            foreach (double sign in new[] { -1.0, 1.0 })
+            {
+                double low = Math.ScaleB(sign, exponent);
+                DoubleDouble actual = DDMath.InvSqrt(DoubleDouble.FromComponents(1.0, low));
+                actual.High.ShouldBe(1.0);
+                actual.Low.ShouldBe(Math.ScaleB(-sign, exponent - 1));
+                DoubleDouble.IsCanonical(actual).ShouldBeTrue();
+            }
+        }
+    }
+
+    private static void AssertInvSqrtBound(DoubleDouble input, DoubleDouble actual)
+    {
+        DoubleDouble.IsFinite(actual).ShouldBeTrue();
+        DoubleDouble.IsCanonical(actual).ShouldBeTrue();
+        (actual.High > 0.0).ShouldBeTrue();
+
+        // Exact dyadic oracle: x = X*2^-1074, y = Y*2^-1074. For positive
+        // x and y, |y/(1/sqrt(x)) - 1| <= t iff (1-t)^2 <= x*y^2 <= (1+t)^2.
+        // Cross-multiply with t=2^-100; no DD product, division or rounded root.
+        BigInteger x = Units(input.High) + Units(input.Low);
+        BigInteger y = Units(actual.High) + Units(actual.Low);
+        BigInteger scale = BigInteger.One << 100;
+        BigInteger product = (x * y * y) << 200;
+        BigInteger lower = ((scale - 1) * (scale - 1)) << 3222;
+        BigInteger upper = ((scale + 1) * (scale + 1)) << 3222;
+        (product >= lower && product <= upper).ShouldBeTrue(
+            $"InvSqrt 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/PreciseMathSqrtTests.cs

@@ -116,6 +116,8 @@ public class PreciseMathSqrtTests
             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(double.MaxValue, Math.BitDecrement(Math.ScaleB(1.0, 970))),
+            DoubleDouble.FromComponents(double.MaxValue, -Math.ScaleB(1.0, 970)),
             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)),
@@ -126,6 +128,84 @@ public class PreciseMathSqrtTests
         }
     }
 
+    [Fact]
+    public void RootRoundingMidpointsAndAdjacentLowsMeetTheBound()
+    {
+        // Exact squares of root midpoints 1 + 2^-53 and 2 - 2^-53:
+        // (1 + 2^-52) + 2^-106 and (4 - 2^-51) + 2^-106.
+        // Bracket each with adjacent low doubles, then exercise exponent parity
+        // and rescaling. At tiny exponents the input itself loses low precision;
+        // the oracle always checks the exact stored input, not the unscaled square.
+        double midpointLow = Math.ScaleB(1.0, -106);
+        for (int exponent = -1074; exponent <= 1022; ++exponent)
+        {
+            foreach (double high in new[] { Math.BitIncrement(1.0), Math.BitDecrement(4.0) })
+            {
+                foreach (double low in new[] { Math.BitDecrement(midpointLow), midpointLow, Math.BitIncrement(midpointLow) })
+                {
+                    DoubleDouble input = DoubleDouble.FromComponents(Math.ScaleB(high, exponent), Math.ScaleB(low, exponent));
+                    if (DoubleDouble.IsFinite(input) && input.High > 0.0)
+                    {
+                        AssertSqrtBound(input, DDMath.Sqrt(input));
+                    }
+                }
+            }
+        }
+    }
+
+    [Fact]
+    public void PublicFactoryNormalizesLegacyZeroHighInputsBeforeTakingTheRoot()
+    {
+        // Legacy raw (0, low) examples are normalized at the public boundary.
+        foreach (double low in new[] { double.Epsilon, 1e-308, 2.0, double.MaxValue })
+        {
+            DoubleDouble input = DoubleDouble.FromComponents(0.0, low);
+            AssertSqrtBound(input, DDMath.Sqrt(input));
+        }
+        DoubleDouble.IsNaN(DDMath.Sqrt(DoubleDouble.FromComponents(0.0, -1.0))).ShouldBeTrue();
+    }
+
+    [Fact]
+    public void DenseLowNormalizationBoundariesMeetTheBound()
+    {
+        // Half-ulp lows and their neighbors test both sides of input normalization.
+        Random random = new(161803);
+        for (int exponent = -1074; exponent <= 1023; ++exponent)
+        {
+            double high = Math.ScaleB(1.0 + random.NextDouble(), exponent);
+            double halfUlp = Math.ScaleB(1.0, exponent - 53);
+            foreach (double magnitude in new[] { Math.BitDecrement(halfUlp), halfUlp, Math.BitIncrement(halfUlp) })
+            {
+                foreach (double low in new[] { magnitude, -magnitude })
+                {
+                    DoubleDouble input = DoubleDouble.FromComponents(high, low);
+                    if (DoubleDouble.IsFinite(input) && input.High > 0.0)
+                    {
+                        AssertSqrtBound(input, DDMath.Sqrt(input));
+                    }
+                }
+            }
+        }
+    }
+
+    [Fact]
+    public void RepresentableSparseCorrectionsAreNotDiscarded()
+    {
+        // sqrt(1+d) = 1+d/2+O(d^2). For these dyadic d, the quadratic term
+        // is below half an ulp of d/2, including when d/2 is the minimum subnormal.
+        foreach (int exponent in new[] { -100, -500, -1000, -1073 })
+        {
+            foreach (double sign in new[] { -1.0, 1.0 })
+            {
+                double low = Math.ScaleB(sign, exponent);
+                DoubleDouble actual = DDMath.Sqrt(DoubleDouble.FromComponents(1.0, low));
+                actual.High.ShouldBe(1.0);
+                actual.Low.ShouldBe(Math.ScaleB(sign, exponent - 1));
+                DoubleDouble.IsCanonical(actual).ShouldBeTrue();
+            }
+        }
+    }
+
     private static void AssertSqrtBound(DoubleDouble input, DoubleDouble actual)
     {
         DoubleDouble.IsFinite(actual).ShouldBeTrue();

README.md

@@ -88,6 +88,12 @@ The `DDMath` static class provides:
   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.
+- `InvSqrt(DoubleDouble)`: computes the reciprocal square root with power-of-two
+  scaling and a compensated Newton step, avoiding double-double division and its
+  allocating boundary paths. The estimate's squared-product residual is retained
+  with FMA. Signed zeros map to correspondingly signed infinities; positive infinity
+  maps to positive zero; negative nonzero inputs and NaN return canonical NaN.
+  This is a dedicated algorithm, not a claim of measured speedup over `1.0 / Sqrt(x)`.
 - `Pow(DoubleDouble, int)`: exponentiation by squaring with a separately tracked
   binary exponent. Supports the full `int` domain, including `int.MinValue`, and
   reciprocates a bounded significand before final scaling for negative exponents.
@@ -115,6 +121,7 @@ The `DDMath` static class provides:
 using Just.PreciseMath;
 
 DoubleDouble root = DDMath.Sqrt(new DoubleDouble(2.0));
+DoubleDouble inverseRoot = DDMath.InvSqrt(new DoubleDouble(2.0));
 DoubleDouble magnitude = DDMath.Abs(-root);
 DoubleDouble smallPower = DDMath.Pow(new DoubleDouble(2.0), -1024);
 DoubleDouble exponential = DDMath.Exp(new DoubleDouble(1.0));
@@ -124,10 +131,12 @@ DoubleDouble preciseExponent = DoubleDouble.FromComponents(0.5, 1e-30);
 DoubleDouble precisePower = DDMath.Pow(new DoubleDouble(2.0), preciseExponent);
 ```
 
-Square-root tests compare the exact component sum against a `2^-100` relative
-error bound using integer inequalities. They include samples at every binary64
+Square-root and inverse-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
+squares/powers of four. Square-root tests also bracket exact root-rounding midpoints;
+both suites exercise half-ulp low-component normalization boundaries.
+This is a tested approximate-accuracy contract, not exhaustive coverage
 of all component pairs or a guarantee of correctly rounded results.
 
 Logarithm tests compare exact component sums with independently generated