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

using System.Globalization;
using System.Numerics;
using Just.PreciseMath.Tests.ReferenceData;
using Shouldly;
using Xunit;

namespace Just.PreciseMath.Tests;

public class PreciseMathRealPowTests
{
    public static IEnumerable<TheoryDataRow<string, double, double, double, double, bool, string>> ReferenceCases => RealPowReferenceData.Cases();

    [Theory]
    [MemberData(nameof(ReferenceCases))]
    public void FinitePowersMeetIndependentComponentAwareReferenceBound(string label,
        double high, double low, double exponentHigh, double exponentLow, bool scalar, string reference)
    {
        // The stdlib generator constructs exact binary64 sums at 2200 digits,
        // checks Fraction equality, evaluates Decimal ln/exp at 450/650 digits,
        // and requires stable 120-digit references. Scalar rows have independent
        // references formed AFTER discarding exponentLow, not the DD reference.
        DoubleDouble input = DoubleDouble.FromComponents(high, low);
        DoubleDouble exponent = DoubleDouble.FromComponents(exponentHigh, scalar ? 0.0 : exponentLow);
        (Units(input.High) + Units(input.Low)).ShouldBe(Units(high) + Units(low));
        BigInteger exponentUnits = Units(exponentHigh) + Units(scalar ? 0.0 : exponentLow);
        (Units(exponent.High) + Units(exponent.Low)).ShouldBe(exponentUnits);
        DoubleDouble actual = scalar ? DDMath.Pow(input, exponentHigh) : DDMath.Pow(input, exponent);
        string context = $"{label}, scalar={scalar}: ({high:R}, {low:R})^({exponentHigh:R}, {exponentLow:R})";
        AssertReferenceBound(actual, reference, exponentUnits, context);
    }

    [Theory]
    [InlineData(4.0, 0.5, 2.0)]
    [InlineData(16.0, 0.25, 2.0)]
    [InlineData(16.0, -0.25, 0.5)]
    [InlineData(0.0625, -0.5, 4.0)]
    public void ExactBinaryReferencesMeetTheRealPathBound(double input, double exponent, double expected)
    {
        // These references are exact rational identities. The general ln/exp
        // implementation is approximate, so do not require universal bit equality.
        string reference = expected.ToString("e17", CultureInfo.InvariantCulture);
        AssertReferenceBound(DDMath.Pow(new DoubleDouble(input), exponent), reference,
            Units(exponent), "scalar exact binary reference");
        AssertReferenceBound(DDMath.Pow(new DoubleDouble(input), new DoubleDouble(exponent)), reference,
            Units(exponent), "DD exact binary reference");
    }

    [Fact]
    public void FiniteIntegerDispatchPreservesExistingPowerBits()
    {
        // Consistency is distinct from accuracy: ReferenceCases independently
        // check the exponent-dependent bound for int-range dispatch too.
        DoubleDouble[] inputs = [DoubleDouble.FromComponents(1.25, -Math.ScaleB(1.0, -55)),
            DoubleDouble.FromComponents(-1.25, Math.ScaleB(1.0, -55)),
            DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -54)),
            DoubleDouble.FromComponents(1.0, -Math.ScaleB(1.0, -54))];
        foreach (DoubleDouble input in inputs)
        {
            foreach (int exponent in new[] { -63, -7, -1, 0, 1, 7, 63 })
            {
                DoubleDouble expected = DDMath.Pow(input, exponent);
                AssertBits(DDMath.Pow(input, (double)exponent), expected);
                AssertBits(DDMath.Pow(input, new DoubleDouble(exponent)), expected);
            }
        }
        foreach (int exponent in new[] { int.MinValue, int.MaxValue })
        {
            DoubleDouble input = DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -54));
            DoubleDouble expected = DDMath.Pow(input, exponent);
            AssertBits(DDMath.Pow(input, (double)exponent), expected);
            AssertBits(DDMath.Pow(input, new DoubleDouble(exponent)), expected);
        }
    }

    [Fact]
    public void ExponentOnePreservesSparseFiniteComponentsExactly()
    {
        foreach (double high in new[] { double.MaxValue, -double.MaxValue, 1.0, -1.0 })
        {
            foreach (double low in new[] { double.Epsilon, -double.Epsilon })
            {
                DoubleDouble input = DoubleDouble.FromComponents(high, low);
                AssertBits(DDMath.Pow(input, 1.0), input);
                AssertBits(DDMath.Pow(input, DoubleDouble.One), input);
            }
        }
    }

    private static void AssertReferenceBound(DoubleDouble actual, string reference,
        BigInteger exponentUnits, string context)
    {
        DoubleDouble.IsCanonical(actual).ShouldBeTrue(context);
        DoubleDouble.IsFinite(actual).ShouldBeTrue(context);
        string[] parts = reference.Split('e');
        int point = parts[0].IndexOf('.', StringComparison.Ordinal);
        int decimals = point < 0 ? 0 : parts[0].Length - point - 1;
        BigInteger numerator = BigInteger.Parse(parts[0].Replace(".", "", StringComparison.Ordinal), CultureInfo.InvariantCulture);
        int exponent = int.Parse(parts[1], CultureInfo.InvariantCulture) - decimals;
        BigInteger denominator = BigInteger.One;
        if (exponent >= 0)
        {
            numerator *= BigInteger.Pow(10, exponent);
        }
        else
        {
            denominator = BigInteger.Pow(10, -exponent);
        }
        BigInteger oneUnits = BigInteger.One << 1074;
        BigInteger integralExponent = BigInteger.DivRem(exponentUnits, oneUnits, out BigInteger remainder);
        BigInteger factor = 1024;
        if (remainder.IsZero && integralExponent >= int.MinValue && integralExponent <= int.MaxValue)
        {
            factor = BigInteger.Max(factor, BigInteger.Abs(integralExponent) + 1);
        }
        BigInteger magnitude = BigInteger.Abs(numerator);
        BigInteger actualUnits = Units(actual.High) + Units(actual.Low);
        BigInteger error = BigInteger.Abs((actualUnits * denominator) - (numerator << 1074));
        // max(2^-90, (|n|+1)*2^-100 for int dispatch) plus 2^-1074.
        // Add the independently checked reference uncertainty of 2^-350 relative.
        // All comparisons retain both result components as exact rational units.
        BigInteger bound = (magnitude * factor << 1324) + (denominator << 350) + (magnitude << 1074);
        ((error << 350) <= bound).ShouldBeTrue(
            $"Real Pow bound failed for {context}: ({actual.High:R}, {actual.Low:R}), reference={reference}");
    }

    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 actual, DoubleDouble expected)
    {
        DoubleDouble.IsCanonical(actual).ShouldBeTrue();
        BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(BitConverter.DoubleToInt64Bits(expected.High));
        BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(BitConverter.DoubleToInt64Bits(expected.Low));
    }
}