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

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