using System.Numerics;
using Shouldly;
using Xunit;

namespace Just.PreciseMath.Tests;

public class DoubleDoubleArithmeticTests
{
    [Fact]
    public void CancellationRetainsBothLowSumTerms()
    {
        double small = Math.ScaleB(1.0, -54);
        double tiny = Math.ScaleB(1.0, -108);
        DoubleDouble left = DoubleDouble.FromComponents(1.0, small);
        DoubleDouble right = DoubleDouble.FromComponents(-1.0, tiny);
        Check(left + right, small, tiny);
        Check(right + left, small, tiny);
        Check(left - (-right), small, tiny);
        Check(+left, 1.0, small);
        Check(-left, -1.0, -small);
        Check(left - left, 0.0, 0.0);
    }

    [Fact]
    public void ScalarOverloadsPreserveOperandOrderAndResiduals()
    {
        DoubleDouble value = DoubleDouble.FromComponents(2.0, Math.ScaleB(1.0, -80));
        Check(3.0 - value, 1.0, -value.Low);
        Check(value - 3.0, -1.0, value.Low);
        Check(value + 3.0, 5.0, value.Low);
        Check(3.0 + value, 5.0, value.Low);
        Check(value * 2.0, 4.0, 2.0 * value.Low);
        Check(2.0 * value, 4.0, 2.0 * value.Low);
        Check(value / 2.0, 1.0, value.Low / 2.0);
        Check(6.0 / new DoubleDouble(2.0), 3.0, 0.0);
    }

    [Fact]
    public void ScalarLeftSubtractionAppliesTheRequestedOperandOrder()
    {
        // Review-1 §1: the former operator -(double, DoubleDouble) returned arg - lvalue,
        // so 3.0 - DD(2.0) produced -1 instead of 1.
        Check(3.0 - new DoubleDouble(2.0), 1.0, 0.0);
        Check(new DoubleDouble(2.0) - 3.0, -1.0, 0.0);
        Check(3.0 - DoubleDouble.FromComponents(2.0, Math.ScaleB(1.0, -80)), 1.0, -Math.ScaleB(1.0, -80));
        Check(DoubleDouble.FromComponents(2.0, Math.ScaleB(1.0, -80)) - 3.0, -1.0, Math.ScaleB(1.0, -80));
    }

    [Fact]
    public void DivisionByAValueCarriedOnlyInTheLowComponentIsFinite()
    {
        // Review-2 §5: with high == 0 and low != 0 the former division returned
        // Infinity because it divided by the zero high component. The public factory
        // now folds such a pair into its high component, so the quotient is finite.
        DoubleDouble denominator = DoubleDouble.FromComponents(0.0, 1.0);
        Check(denominator, 1.0, 0.0);
        Check(new DoubleDouble(4.0) / denominator, 4.0, 0.0);
        Check(4.0 / denominator, 4.0, 0.0);
        Check(new DoubleDouble(4.0) / DoubleDouble.FromComponents(0.0, -2.0), -2.0, 0.0);
    }

    [Fact]
    public void ScalarCancellationPreservesTheRemainingExpansionInBothOrders()
    {
        // At a normal binade boundary, 2^e - BitDecrement(2^e) = 2^(e-53).
        // The low input becomes the representable residual of that exact difference.
        foreach (int exponent in new[] { -967, -900, -450, 0, 450, 969, 1020, 1023 })
        {
            foreach (double sign in new[] { -1.0, 1.0 })
            {
                double high = Math.ScaleB(1.0, exponent);
                double low = sign * Math.ScaleB(1.0, exponent - 107);
                double scalar = sign * Math.BitDecrement(high);
                double difference = sign * Math.ScaleB(1.0, exponent - 53);
                DoubleDouble value = DoubleDouble.FromComponents(sign * high, low);
                Check(value - scalar, difference, low);
                Check(scalar - value, -difference, -low);
                Check(value + (-scalar), difference, low);
                Check((-scalar) + value, difference, low);
            }
        }
    }

    [Fact]
    public void MixedSpecialValuesIgnoreFiniteResidualsButPreserveResultSigns()
    {
        foreach (double sign in new[] { -1.0, 1.0 })
        {
            DoubleDouble value = DoubleDouble.FromComponents(sign, sign * Math.ScaleB(1.0, -54));
            foreach (double scalar in new[] { double.NaN, double.NegativeInfinity, double.PositiveInfinity })
            {
                CheckBits(value + scalar, sign + scalar);
                CheckBits(scalar + value, scalar + sign);
                CheckBits(value - scalar, sign - scalar);
                CheckBits(scalar - value, scalar - sign);
                CheckBits(value * scalar, sign * scalar);
                CheckBits(scalar * value, scalar * sign);
                CheckBits(value / scalar, sign / scalar);
                CheckBits(scalar / value, scalar / sign);
            }
            foreach (double zero in new[] { 0.0, -0.0 })
            {
                CheckBits(value * zero, sign * zero);
                CheckBits(zero * value, zero * sign);
                CheckBits(value / zero, sign / zero);
                CheckBits(zero / value, zero / sign);
                Check(value + zero, value.High, value.Low);
                Check(zero + value, value.High, value.Low);
                Check(value - zero, value.High, value.Low);
                Check(zero - value, -value.High, -value.Low);
            }
        }
    }

    [Fact]
    public void ProductAndQuotientRetainExtraPrecision()
    {
        // (1 + 2^-52)(1 - 2^-52) = 1 - 2^-104 exactly.
        Check(new DoubleDouble(1.0 + Math.ScaleB(1.0, -52)) * new DoubleDouble(1.0 - Math.ScaleB(1.0, -52)),
            1.0, -Math.ScaleB(1.0, -104));
        // Binary expansion of 1/3, rounding high then residual ties-to-even.
        Check(DoubleDouble.One / 3.0, 0.3333333333333333, 1.850371707708594e-17);
    }

    [Fact]
    public void ScalarProductNormalizesACorrectionBeyondTheHighMidpoint()
    {
        // (1 + 2^-53)(1 + 2^-52) = 1 + 3*2^-53 + 2^-105, exactly.
        // The rounded high advances twice above 1; its residual is still exact.
        DoubleDouble value = DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -53));
        double scalar = Math.BitIncrement(1.0);
        double high = 1.0 + Math.ScaleB(1.0, -51);
        double low = -Math.ScaleB(1.0, -53) + Math.ScaleB(1.0, -105);
        Check(value * scalar, high, low);
        Check(scalar * value, high, low);
        Check(value * (-scalar), -high, -low);
        Check((-scalar) * value, -high, -low);
    }

    [Fact]
    public void ScalarDivisionRetainsNumeratorAndDenominatorResiduals()
    {
        double low = Math.ScaleB(1.0, -80);
        DoubleDouble value = DoubleDouble.FromComponents(1.0, low);
        Check(value / 1.0, 1.0, low);
        Check(value / (-1.0), -1.0, -low);
        // Compare the reciprocal to its exact rational, not another DD operator.
        BigInteger numerator = BigInteger.One << 2148;
        AssertRelative(1.0 / value, numerator, Units(value));
        AssertRelative(-1.0 / value, -numerator, Units(value));
    }

    [Fact]
    public void ExtremeFiniteOperationsDoNotOverflowIntermediates()
    {
        DoubleDouble maximum = new(double.MaxValue);
        DoubleDouble third = maximum / 3.0;
        DoubleDouble thirdPair = maximum / new DoubleDouble(3.0);
        DoubleDouble thirdScalar = double.MaxValue / new DoubleDouble(3.0);
        AssertRelative(third, Units(maximum), 3);
        thirdPair.ShouldBe(third);
        thirdScalar.ShouldBe(third);
        Check(new DoubleDouble(double.Epsilon) / new DoubleDouble(double.Epsilon), 1.0, 0.0);
        Check(new DoubleDouble(double.Epsilon) * new DoubleDouble(Math.ScaleB(1.0, 1023)), Math.ScaleB(1.0, -51), 0.0);
        Check(new DoubleDouble(Math.ScaleB(1.0, -1022)) / 2.0, Math.ScaleB(1.0, -1023), 0.0);
        Check(maximum * 2.0, double.PositiveInfinity, 0.0);
        Check(maximum + maximum, double.PositiveInfinity, 0.0);
        // High-only addition overflows, but the complete sum is exactly MaxValue.
        DoubleDouble below = DoubleDouble.FromComponents(double.MaxValue, -Math.ScaleB(1.0, 969));
        Check(below + Math.ScaleB(1.0, 969), double.MaxValue, 0.0);
    }

    [Theory]
    [InlineData(1.0)]
    [InlineData(-1.0)]
    public void ExactOverflowMidpointStillOverflows(double sign)
    {
        DoubleDouble maximum = new(sign * double.MaxValue);
        double halfUlp = sign * Math.ScaleB(1.0, 970);
        Check(maximum + halfUlp, sign * double.PositiveInfinity, 0.0);
        Check(maximum - (-halfUlp), sign * double.PositiveInfinity, 0.0);
        Check(DoubleDouble.FromComponents(sign * double.MaxValue, halfUlp), sign * double.PositiveInfinity, 0.0);
    }

    [Fact]
    public void SpecialValueMatrixMatchesBinary64IncludingZeroSigns()
    {
        double[] values = [0.0, -0.0, 1.0, -1.0, double.PositiveInfinity, double.NegativeInfinity, double.NaN];
        foreach (double left in values)
        {
            foreach (double right in values)
            {
                DoubleDouble a = new(left);
                DoubleDouble b = new(right);
                CheckBits(a + b, left + right);
                CheckBits(a - b, left - right);
                CheckBits(a * b, left * right);
                CheckBits(a / b, left / right);
                CheckBits(a + right, left + right);
                CheckBits(left + b, left + right);
                CheckBits(a - right, left - right);
                CheckBits(left - b, left - right);
                CheckBits(a * right, left * right);
                CheckBits(left * b, left * right);
                CheckBits(a / right, left / right);
                CheckBits(left / b, left / right);
            }
        }
    }

    [Fact]
    public void DeterministicArithmeticMeetsConservativeErrorBound()
    {
        Random random = new(1729);
        for (int i = 0; i < 250; ++i)
        {
            DoubleDouble a = DoubleDouble.FromComponents(Math.ScaleB((random.NextDouble() * 2.0) - 1.0, random.Next(-400, 401)),
                Math.ScaleB(random.NextDouble(), random.Next(-500, -450)));
            DoubleDouble b = DoubleDouble.FromComponents(Math.ScaleB((random.NextDouble() * 2.0) - 1.0, random.Next(-400, 401)),
                Math.ScaleB(random.NextDouble(), random.Next(-500, -450)));
            BigInteger x = Units(a);
            BigInteger y = Units(b);
            AssertRelative(a + b, x + y, BigInteger.One);
            AssertRelative(a - b, x - y, BigInteger.One);
            AssertRelative(a * b, x * y, BigInteger.One << 1074);
            AssertRelative(a / b, x << 1074, y);
        }
    }

    private static void Check(DoubleDouble value, double high, double low)
    {
        value.High.ShouldBe(high);
        value.Low.ShouldBe(low);
    }

    private static void CheckBits(DoubleDouble value, double expected)
    {
        if (double.IsNaN(expected))
        {
            DoubleDouble.IsNaN(value).ShouldBeTrue();
        }
        else
        {
            BitConverter.DoubleToInt64Bits(value.High).ShouldBe(BitConverter.DoubleToInt64Bits(expected));
        }
        BitConverter.DoubleToInt64Bits(value.Low).ShouldBe(0L);
    }

    // Independent oracle: every finite binary64 is an integer multiple of 2^-1074.
    private static BigInteger Units(DoubleDouble value)
    {
        return Units(value.High) + Units(value.Low);
    }

    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 AssertRelative(DoubleDouble actual, BigInteger numerator, BigInteger denominator)
    {
        DoubleDouble.IsFinite(actual).ShouldBeTrue();
        BigInteger error = BigInteger.Abs((Units(actual) * denominator) - numerator);
        // <= 2^-100 relative error plus one minimum subnormal (rounding floor).
        (error <= (BigInteger.Abs(numerator) >> 100) + BigInteger.Abs(denominator)).ShouldBeTrue();
        if (actual.High != 0.0)
        {
            (Math.Abs(actual.Low) <= Math.ScaleB(1.0, Math.ILogB(actual.High) - 53)).ShouldBeTrue();
        }
    }
}