namespace Just.PreciseMath.Tests;

public class PreciseMathPowTests
{
    [Fact]
    public void NegativePowerDoesNotOverflowBeforeReciprocation()
    {
        // Legacy review: Pow(2, -1024) overflowed before taking the reciprocal.
        DoubleDouble actual = DDMath.Pow(new DoubleDouble(2.0), -1024);
        AssertBits(actual, Math.ScaleB(1.0, -1024), 0.0);
    }

    [Fact]
    public void MinimumIntegerExponentDoesNotOverflowItsMagnitude()
    {
        AssertBits(DDMath.Pow(DoubleDouble.One, int.MinValue), 1.0, 0.0);
        AssertBits(DDMath.Pow(DoubleDouble.NegativeOne, int.MinValue), 1.0, 0.0);
        AssertBits(DDMath.Pow(new DoubleDouble(2.0), int.MinValue), 0.0, 0.0);
        AssertBits(DDMath.Pow(new DoubleDouble(0.5), int.MinValue), double.PositiveInfinity, 0.0);
    }

    [Fact]
    public void PowerOnePreservesEvenAnExtremelySparseLow()
    {
        foreach (double sign in new[] { -1.0, 1.0 })
        {
            DoubleDouble input = DoubleDouble.FromComponents(sign * double.MaxValue, sign * double.Epsilon);
            AssertBits(DDMath.Pow(input, 1), input.High, input.Low);
        }
    }

    [Fact]
    public void SpecialValuesFollowIntegerPowerRulesWithCanonicalBits()
    {
        double[] inputs = [0.0, -0.0, double.PositiveInfinity, double.NegativeInfinity,
            double.NaN, 1.0, -1.0];
        int[] exponents = [int.MinValue, int.MinValue + 1, -1025, -2, -1, 0, 1, 2, 1025, int.MaxValue];
        foreach (double input in inputs)
        {
            foreach (int exponent in exponents)
            {
                // Every int converts exactly to double. BCL is an independent
                // special-value/sign oracle, not the finite precision oracle.
                double expected = Math.Pow(input, exponent);
                AssertBits(DDMath.Pow(new DoubleDouble(input), exponent),
                    double.IsNaN(expected) ? double.NaN : expected, 0.0);
            }
        }
    }

    [Fact]
    public void BinaryPowersAreExactAcrossInputAndOutputExponentRanges()
    {
        int[] powers = [int.MinValue, int.MinValue + 1, -1075, -1024, -3, -2, -1,
            0, 1, 2, 3, 1024, 1075, int.MaxValue];
        for (int inputExponent = -1074; inputExponent <= 1023; ++inputExponent)
        {
            foreach (int power in powers)
            {
                long outputExponent = (long)inputExponent * power;
                double magnitude = outputExponent > 1023 ? double.PositiveInfinity
                    : outputExponent < -1074 ? 0.0 : Math.ScaleB(1.0, (int)outputExponent);
                foreach (double sign in new[] { -1.0, 1.0 })
                {
                    double expected = sign < 0.0 && (power & 1) != 0 ? -magnitude : magnitude;
                    AssertBits(DDMath.Pow(new DoubleDouble(sign * Math.ScaleB(1.0, inputExponent)), power),
                        expected, 0.0);
                }
            }
        }
    }

    [Fact]
    public void FinitePowersMeetAnIndependentExactRationalErrorBound()
    {
        // Exact binary input sums, raised with BigInteger arithmetic; no DD
        // multiplication/division, decimal, or Math.Pow supplies the oracle.
        Random random = new(270718);
        for (int sample = 0; sample < 100; ++sample)
        {
            int scale = random.Next(-12, 13);
            double high = Math.ScaleB(1.0 + random.NextDouble(), scale);
            double low = Math.ScaleB(random.NextDouble(), scale - 54);
            foreach (double residual in new[] { 0.0, low, -low })
            {
                DoubleDouble input = DoubleDouble.FromComponents(sample % 2 == 0 ? high : -high, residual);
                foreach (int exponent in new[] { -63, -7, -2, -1, 1, 2, 7, 63 })
                {
                    AssertExactPowerBound(input, exponent);
                }
            }
        }
        foreach (DoubleDouble input in new[] {
            DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -53)),
            DoubleDouble.FromComponents(1.0, -Math.ScaleB(1.0, -54)),
            DoubleDouble.FromComponents(0.5, double.Epsilon) })
        {
            foreach (int exponent in new[] { -1024, -257, 257, 1024 })
            {
                if (input.High == 0.5 && exponent == -1024)
                {
                    continue; // Overflow is pinned by the binary-power matrix.
                }
                AssertExactPowerBound(input, exponent);
            }
        }
    }

    [Fact]
    public void NegativePowersRetainRepresentableLowsAtExtremeInputMagnitudes()
    {
        foreach (int scale in new[] { -1023, -1000, -500, 500, 1000, 1023 })
        {
            foreach (double sign in new[] { -1.0, 1.0 })
            {
                DoubleDouble input = DoubleDouble.FromComponents(Math.ScaleB(1.5, scale),
                    sign * Math.ScaleB(1.0, scale - 54));
                AssertExactPowerBound(input, -1);
                if (scale == -500 || scale == 500)
                {
                    DoubleDouble actual = DDMath.Pow(input, -2);
                    actual.Low.ShouldNotBe(0.0);
                    AssertExactPowerBound(input, -2);
                }
                if (scale < 0)
                {
                    DDMath.Pow(input, -1).Low.ShouldNotBe(0.0);
                }
            }
        }
        AssertExactPowerBound(new DoubleDouble(double.MaxValue), -1);
        AssertExactPowerBound(new DoubleDouble(1e-308), -1);
        AssertExactPowerBound(new DoubleDouble(1e308), -1);
    }

    [Fact]
    public void SparseReciprocalCorrectionsSurviveBothRangeDirections()
    {
        // 1/(h+l) = 1/h - l/h^2 + O(l^2/h^3). In these dyadic
        // cases the remainder is strictly below half an ulp of the stated low.
        foreach (double sign in new[] { -1.0, 1.0 })
        {
            DoubleDouble tiny = DoubleDouble.FromComponents(Math.ScaleB(1.0, -500), sign * double.Epsilon);
            AssertBits(DDMath.Pow(tiny, -1), Math.ScaleB(1.0, 500), -sign * Math.ScaleB(1.0, -74));
            DoubleDouble huge = DoubleDouble.FromComponents(Math.ScaleB(1.0, 500), sign * Math.ScaleB(1.0, -74));
            AssertBits(DDMath.Pow(huge, -1), Math.ScaleB(1.0, -500), -sign * double.Epsilon);
            DoubleDouble tinySquare = DoubleDouble.FromComponents(Math.ScaleB(1.0, -250), sign * double.Epsilon);
            AssertBits(DDMath.Pow(tinySquare, -2), Math.ScaleB(1.0, 500), -sign * Math.ScaleB(1.0, -323));
        }
    }

    [Fact]
    public void UnderflowMidpointAndAdjacentLowComponentsRespectBothSigns()
    {
        // (2^-215)^5 = 2^-1075, exactly the tie between zero and epsilon.
        // The low perturbations strictly bracket the tie; reciprocals of
        // the corresponding large inputs reverse the perturbation direction.
        foreach (int exponent in new[] { -5, 5 })
        {
            int scale = exponent < 0 ? 215 : -215;
            foreach (double lowSign in new[] { -1.0, 0.0, 1.0 })
            {
                bool aboveTie = lowSign * exponent > 0.0;
                foreach (double sign in new[] { -1.0, 1.0 })
                {
                    DoubleDouble input = DoubleDouble.FromComponents(sign * Math.ScaleB(1.0, scale),
                        sign * lowSign * Math.ScaleB(1.0, scale - 54));
                    AssertBits(DDMath.Pow(input, exponent), sign * (aboveTie ? double.Epsilon : 0.0), 0.0);
                }
            }
        }
    }

    [Fact]
    public void OverflowBoundaryUsesTheLowComponentBeforeFinalScaling()
    {
        foreach (int exponent in new[] { -4, 4 })
        {
            int scale = exponent < 0 ? -256 : 256;
            double lowSign = exponent < 0 ? 1.0 : -1.0;
            DoubleDouble finite = DoubleDouble.FromComponents(Math.ScaleB(1.0, scale),
                lowSign * Math.ScaleB(1.0, scale - 55));
            AssertExactPowerBound(finite, exponent);
            DoubleDouble overflowing = DoubleDouble.FromComponents(Math.ScaleB(1.0, scale),
                lowSign * Math.ScaleB(1.0, scale - 57));
            AssertBits(DDMath.Pow(overflowing, exponent), double.PositiveInfinity, 0.0);
        }
    }

    [Theory]
    // ReferenceData/generate_pow_reference.py: exact stored binary inputs,
    // Decimal integer powers at 160/240 digits, stable floor(value * 10^100).
    [InlineData(1.0, 5.551115123125783e-17, -2147483648, "9999998807907175546458285654165775755995358529592462953428551373115894598504056800747655101172369356")]
    [InlineData(1.0, 5.551115123125783e-17, 2147483647, "10000001192092966006977287879000612242425055148009778363959373202290565261791488930669591094667385585")]
    [InlineData(1.0, -5.551115123125783e-17, -2147483648, "10000001192092966562088932540488693037989493813706046557380379233423355139596032355926918800247828081")]
    [InlineData(1.0, -5.551115123125783e-17, 2147483647, "9999998807907176101569665617857901211208198096053270077732977033266503983279402161561236649189617033")]
    [InlineData(1.0000000009313226, 5.169878828456423e-26, -2147483648, "1353352833626534812694377803744507919692475898909062344519276068330639808800815294171180036409012889")]
    [InlineData(1.0000000009313226, 5.169878828456423e-26, 2147483647, "73890560851674615678069694043352842890214309016451211229831428411840370567020651340427147012473005377")]
    [InlineData(0.9999999990686774, -5.169878828456423e-26, -2147483648, "73890561058122458050918958513594285934973647528626440535855911314733978878721028315311791064891371241")]
    [InlineData(0.9999999990686774, -5.169878828456423e-26, 2147483647, "1353352832366126768492006308656122646859542550999041350022637326283256718420076811410419433559508500")]
    [InlineData(1.00000003, 2.7755575615628914e-17, -2147483648, "1049038452726052018566777293807106599923174085199532650373208349817256939")]
    [InlineData(1.00000003, 2.7755575615628914e-17, 2147483647, "95325387491886611267179809406372425766774105569877060927978408348321633401745440609674655439897287007094693160193846818457460816")]
    public void ExtremeIntegerExponentsMeetHighPrecisionReferenceIntervals(
        double high, double low, int exponent, string lowerNumerator)
    {
        BigInteger numerator = BigInteger.Parse(lowerNumerator, System.Globalization.CultureInfo.InvariantCulture);
        BigInteger denominator = BigInteger.Pow(10, 100);
        foreach (double sign in new[] { -1.0, 1.0 })
        {
            DoubleDouble input = DoubleDouble.FromComponents(sign * high, sign * low);
            DoubleDouble actual = DDMath.Pow(input, exponent);
            actual.Low.ShouldNotBe(0.0);
            int resultSign = sign < 0.0 && (exponent & 1) != 0 ? -1 : 1;
            // Check both interval endpoints so the reference uncertainty is not
            // silently omitted from the claimed exponent-dependent tolerance.
            AssertRationalBound(actual, resultSign * numerator, denominator, exponent, lowerNumerator);
            AssertRationalBound(actual, resultSign * (numerator + 1), denominator, exponent, lowerNumerator);
        }
    }

    private static void AssertExactPowerBound(DoubleDouble input, int exponent)
    {
        BigInteger numerator = Units(input.High) + Units(input.Low);
        int trailing = (int)BigInteger.TrailingZeroCount(BigInteger.Abs(numerator));
        numerator >>= trailing;
        BigInteger denominator = BigInteger.One;
        if (trailing <= 1074)
        {
            denominator <<= 1074 - trailing;
        }
        else
        {
            numerator <<= trailing - 1074;
        }
        int magnitude = Math.Abs(exponent); // This oracle is used only for small exponents.
        numerator = BigInteger.Pow(numerator, magnitude);
        denominator = BigInteger.Pow(denominator, magnitude);
        if (exponent < 0)
        {
            (numerator, denominator) = (denominator * numerator.Sign, BigInteger.Abs(numerator));
        }
        AssertRationalBound(DDMath.Pow(input, exponent), numerator, denominator, exponent,
            $"({input.High:R}, {input.Low:R})^{exponent}");
    }

    private static void AssertRationalBound(DoubleDouble actual, BigInteger numerator,
        BigInteger denominator, int exponent, string context)
    {
        DoubleDouble.IsFinite(actual).ShouldBeTrue(context);
        DoubleDouble.IsCanonical(actual).ShouldBeTrue(context);
        BigInteger actualUnits = Units(actual.High) + Units(actual.Low);
        BigInteger error = BigInteger.Abs(actualUnits * denominator - (numerator << 1074));
        // |actual-exact| <= |exact|*(|n|+1)*2^-100 + 2^-1074,
        // cross-multiplied exactly in units of the smallest subnormal.
        BigInteger tolerance = (BigInteger.Abs(numerator) * (Math.Abs((long)exponent) + 1) << 974)
            + denominator;
        (error <= tolerance).ShouldBeTrue($"Power bound failed for {context}: ({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 actual, double high, double low)
    {
        BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(BitConverter.DoubleToInt64Bits(high));
        BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(BitConverter.DoubleToInt64Bits(low));
        DoubleDouble.IsCanonical(actual).ShouldBeTrue();
    }
}