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

namespace Just.PreciseMath.Tests;

public class DoubleDoubleCbrtTests
{
    [Fact]
    public void ExactCubesAndSpecialValues()
    {
        foreach (double value in new[] { 0.0, -0.0, double.NaN, double.PositiveInfinity, double.NegativeInfinity })
        {
            DoubleDouble result = DoubleDouble.Cbrt(new DoubleDouble(value));
            AssertBits(new DoubleDouble(Math.Cbrt(value)), result);
        }
        foreach (double root in new[] { -3.0, -2.0, -1.0, 0.5, 1.0, 2.0, 3.0 })
        {
            AssertBits(new DoubleDouble(root), DoubleDouble.Cbrt(new DoubleDouble(root * root * root)));
        }
        for (int exponent = -1074; exponent <= 1023; exponent += 3)
        {
            DoubleDouble input = new(Math.ScaleB(1.0, exponent));
            AssertBits(new DoubleDouble(Math.ScaleB(1.0, exponent / 3)), DoubleDouble.Cbrt(input));
        }
    }

    [Fact]
    public void CubesMeetExactBoundAcrossEveryInputExponent()
    {
        Random random = new(271828);
        for (int exponent = -1074; exponent <= 1023; ++exponent)
        {
            double high = Math.ScaleB(1.0 + random.NextDouble(), exponent);
            double dense = Math.ScaleB(random.NextDouble(), exponent - 53);
            foreach (double low in new[] { 0.0, dense, -dense, double.Epsilon, -double.Epsilon })
            {
                DoubleDouble input = DoubleDouble.FromComponents(high, low);
                if (input.High > 0.0 && DoubleDouble.IsFinite(input))
                {
                    AssertBound(input, DoubleDouble.Cbrt(input));
                    AssertBound(-input, DoubleDouble.Cbrt(-input));
                    AssertBits(DoubleDouble.Cbrt(input), DDMath.Cbrt(input));
                }
            }
        }
        foreach (double low in new[] { 0.0, Math.BitDecrement(Math.ScaleB(1.0, 970)), -Math.ScaleB(1.0, 970) })
        {
            DoubleDouble input = DoubleDouble.FromComponents(double.MaxValue, low);
            AssertBound(input, DoubleDouble.Cbrt(input));
        }
    }

    [Fact]
    public void SparseNearOneCorrectionsSurvive()
    {
        // cbrt(1+d) = 1+d/3+O(d^2); for these d, the quadratic is
        // below half an ulp of the rounded binary64 correction d/3.
        foreach (int exponent in new[] { -100, -500, -1000, -1072 })
        {
            foreach (double sign in new[] { -1.0, 1.0 })
            {
                double low = Math.ScaleB(sign, exponent);
                DoubleDouble result = DoubleDouble.Cbrt(DoubleDouble.FromComponents(1.0, low));
                result.High.ShouldBe(1.0);
                result.Low.ShouldBe(low / 3.0);
            }
        }
    }

    [Fact]
    public void ScalingDoesNotEraseARepresentableSparseCorrection()
    {
        // cbrt(2^900 + d) = 2^300 + d/(3*2^600) + O(d²/2^1500).
        // Scaling d=2^-174 by 2^-900 leaves epsilon, whose division by 3
        // underflows, although the final correction 2^-774/3 is representable.
        foreach (double sign in new[] { -1.0, 1.0 })
        {
            DoubleDouble input = DoubleDouble.FromComponents(Math.ScaleB(1.0, 900), Math.ScaleB(sign, -174));
            DoubleDouble result = DoubleDouble.Cbrt(input);
            AssertBits(DoubleDouble.FromComponents(Math.ScaleB(1.0, 300), Math.ScaleB(sign / 3.0, -774)), result);
            AssertBits(-result, DoubleDouble.Cbrt(-input));
        }
    }

    private static void AssertBound(DoubleDouble input, DoubleDouble actual)
    {
        DoubleDouble.IsFinite(actual).ShouldBeTrue();
        DoubleDouble.IsCanonical(actual).ShouldBeTrue();
        Math.Sign(actual.High).ShouldBe(Math.Sign(input.High));
        // Exact dyadic oracle, independent of the implementation's arithmetic:
        // x*(1-t)^3 <= |y|^3 <= x*(1+t)^3, t=2^-100, x=|input|.
        BigInteger x = BigInteger.Abs(Units(input.High) + Units(input.Low)) << 2148;
        BigInteger y = BigInteger.Abs(Units(actual.High) + Units(actual.Low));
        BigInteger scale = BigInteger.One << 100;
        BigInteger cube = BigInteger.Pow(y, 3) << 300;
        (cube >= x * BigInteger.Pow(scale - 1, 3)
            && cube <= x * BigInteger.Pow(scale + 1, 3)).ShouldBeTrue(
                $"Cbrt 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;
    }

    private static void AssertBits(DoubleDouble expected, DoubleDouble actual)
    {
        BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(BitConverter.DoubleToInt64Bits(expected.High));
        BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(BitConverter.DoubleToInt64Bits(expected.Low));
    }
}