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

namespace Just.PreciseMath.Tests;

public class PreciseMathReciprocalTests
{
    [Theory]
    [InlineData(3.0)]
    [InlineData(7.0)]
    [InlineData(-3.0)]
    [InlineData(-7.0)]
    public void OrdinaryReciprocalsRetainMoreThanBinary64Precision(double high)
    {
        foreach (double low in new[] { 0.0, Math.ScaleB(1.0, -54), -Math.ScaleB(1.0, -54) })
        {
            DoubleDouble input = DoubleDouble.FromComponents(high, low);
            DoubleDouble actual = DDMath.Reciprocal(input);
            actual.Low.ShouldNotBe(0.0);
            AssertReciprocalBound(input, actual);
        }
    }

    [Fact]
    public void SpecialValuesFollowBinary64ReciprocalAndRemainCanonical()
    {
        double[] values = [0.0, -0.0, double.PositiveInfinity, double.NegativeInfinity,
            double.NaN, BitConverter.Int64BitsToDouble(0x7ff0000000000001L)];
        foreach (double value in values)
        {
            DoubleDouble actual = DDMath.Reciprocal(new DoubleDouble(value));
            // Binary64 is an independent oracle for special values only.
            double expected = 1.0 / value;
            BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(
                BitConverter.DoubleToInt64Bits(double.IsNaN(expected) ? double.NaN : expected));
            BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L);
            DoubleDouble.IsCanonical(actual).ShouldBeTrue();
        }
    }

    [Fact]
    public void PowersOfTwoHaveExactReciprocalsOrSignedOverflowAcrossTheRange()
    {
        // 1/(s*2^k) = s*2^-k. All finite results here are exactly representable,
        // including subnormals. Exponents below -1023 overflow the result.
        for (int exponent = -1074; exponent <= 1023; ++exponent)
        {
            foreach (double sign in new[] { -1.0, 1.0 })
            {
                DoubleDouble actual = DDMath.Reciprocal(new DoubleDouble(Math.ScaleB(sign, exponent)));
                actual.High.ShouldBe(Math.ScaleB(sign, -exponent));
                BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L);
                DoubleDouble.IsCanonical(actual).ShouldBeTrue();
            }
        }
    }

    [Fact]
    public void ReciprocalsMatchScalarDivisionBitsAcrossTheRange()
    {
        // Compatibility is bitwise, not just the same accuracy tolerance. Sample
        // every exponent, both signs, dense/sparse lows and binade neighbors;
        // this includes both ends of the scalar division fast-path guard.
        Random random = new(314159);
        for (int exponent = -1074; exponent <= 1023; ++exponent)
        {
            double high = Math.ScaleB(1.0 + random.NextDouble(), exponent);
            double halfUlp = Math.ScaleB(1.0, exponent - 53);
            double power = Math.ScaleB(1.0, exponent);
            foreach (double sign in new[] { -1.0, 1.0 })
            {
                foreach (double low in new[] { 0.0, halfUlp, -halfUlp,
                    Math.BitDecrement(halfUlp), -Math.BitDecrement(halfUlp),
                    Math.BitIncrement(halfUlp), -Math.BitIncrement(halfUlp),
                    double.Epsilon, -double.Epsilon })
                {
                    AssertMatchesDivision(DoubleDouble.FromComponents(sign * high, low));
                }
                foreach (double boundary in new[] { Math.BitDecrement(power), power, Math.BitIncrement(power) })
                {
                    AssertMatchesDivision(new DoubleDouble(sign * boundary));
                }
            }
        }

        DoubleDouble[] specials = [DoubleDouble.Zero, new(-0.0), DoubleDouble.NaN,
            new(double.PositiveInfinity), new(double.NegativeInfinity),
            new(double.MaxValue), new(double.MinValue),
            DoubleDouble.FromComponents(double.MaxValue, Math.BitDecrement(Math.ScaleB(1.0, 970))),
            DoubleDouble.FromComponents(double.MinValue, -Math.BitDecrement(Math.ScaleB(1.0, 970)))];
        foreach (DoubleDouble input in specials)
        {
            AssertMatchesDivision(input);
        }
    }

    [Fact]
    public void RepresentableSparseCorrectionsAreRetained()
    {
        // 1/(1+d) = 1-d+O(d^2). The omitted tail is less than half an ulp
        // of d for these exact dyadics, including the minimum subnormal.
        foreach (int exponent in new[] { -100, -500, -1000, -1074 })
        {
            foreach (double sign in new[] { -1.0, 1.0 })
            {
                foreach (double lowSign in new[] { -1.0, 1.0 })
                {
                    double low = Math.ScaleB(lowSign, exponent);
                    DoubleDouble input = DoubleDouble.FromComponents(sign, low);
                    DoubleDouble actual = DDMath.Reciprocal(input);
                    actual.High.ShouldBe(sign);
                    actual.Low.ShouldBe(-low);
                    AssertReciprocalBound(input, actual);
                }
            }
        }
    }

    private static void AssertMatchesDivision(DoubleDouble input)
    {
        DoubleDouble expected = 1.0 / input;
        DoubleDouble actual = DDMath.Reciprocal(input);
        string context = $"Input: ({input.High:R}, {input.Low:R})";
        BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(BitConverter.DoubleToInt64Bits(expected.High), context);
        BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(BitConverter.DoubleToInt64Bits(expected.Low), context);
        DoubleDouble.IsCanonical(actual).ShouldBeTrue();
        if (DoubleDouble.IsFinite(input) && input.High != 0.0 && DoubleDouble.IsFinite(actual))
        {
            // Division establishes compatibility, not accuracy: check the exact
            // rational inequality independently for every finite matrix result.
            AssertReciprocalBound(input, actual);
        }
    }

    private static void AssertReciprocalBound(DoubleDouble input, DoubleDouble actual)
    {
        DoubleDouble.IsFinite(actual).ShouldBeTrue();
        DoubleDouble.IsCanonical(actual).ShouldBeTrue();
        Math.Sign(actual.High).ShouldBe(Math.Sign(input.High));

        // Exact dyadic oracle: x = X*2^-1074, y = Y*2^-1074.
        // |y - 1/x| <= |1/x|*2^-100 + 2^-1074 is equivalent to
        // |X*Y - 2^2148|*2^100 <= 2^2148 + |X|*2^100.
        BigInteger x = Units(input.High) + Units(input.Low);
        BigInteger y = Units(actual.High) + Units(actual.Low);
        BigInteger scale = BigInteger.One << 2148;
        BigInteger error = BigInteger.Abs((x * y) - scale) << 100;
        (error <= scale + (BigInteger.Abs(x) << 100)).ShouldBeTrue(
            $"Reciprocal 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;
    }
}