using System.Numerics;
using Shouldly;
using Xunit;

namespace Just.PreciseMath.Tests;

public class DoubleDoubleHypotTests
{
    [Fact]
    public void SpecialValuesAndZeroOperandsFollowTheContract()
    {
        double[] values = [0.0, -0.0, 1.0, -1.0, double.PositiveInfinity, double.NegativeInfinity, double.NaN];
        foreach (double x in values)
        {
            foreach (double y in values)
            {
                if (!double.IsFinite(x) || !double.IsFinite(y) || x == 0.0 || y == 0.0)
                {
                    AssertBits(new DoubleDouble(double.Hypot(x, y)), DoubleDouble.Hypot(new DoubleDouble(x), new DoubleDouble(y)));
                }
            }
        }
        DoubleDouble sparse = DoubleDouble.FromComponents(double.MaxValue, double.Epsilon);
        AssertBits(sparse, DoubleDouble.Hypot(-sparse, DoubleDouble.NegativeZero));
        AssertBits(sparse, DoubleDouble.Hypot(DoubleDouble.Zero, -sparse));
        AssertBits(DoubleDouble.PositiveInfinity, DoubleDouble.Hypot(new DoubleDouble(double.MaxValue), new DoubleDouble(double.MaxValue)));
    }

    [Fact]
    public void ScaledPythagoreanTriplesAreExact()
    {
        for (int exponent = -1074; exponent <= 1021; ++exponent)
        {
            DoubleDouble x = new(Math.ScaleB(3.0, exponent));
            DoubleDouble y = new(Math.ScaleB(4.0, exponent));
            DoubleDouble expected = new(Math.ScaleB(5.0, exponent));
            AssertBits(expected, DoubleDouble.Hypot(x, y));
            AssertBits(expected, DoubleDouble.Hypot(-y, -x));
        }
    }

    [Fact]
    public void FiniteResultsMeetExactSquaredBoundAcrossTheRange()
    {
        Random random = new(123456);
        for (int exponent = -1074; exponent <= 1023; ++exponent)
        {
            double high = Math.ScaleB(1.0 + (random.NextDouble() * 0.1), exponent);
            foreach (int gap in new[] { 0, 1, 26, 54, 55, 500, 1100 })
            {
                double second = Math.ScaleB(1.0, exponent - gap);
                foreach (double sign in new[] { -1.0, 1.0 })
                {
                    DoubleDouble x = DoubleDouble.FromComponents(high, Math.ScaleB(sign, exponent - 54));
                    DoubleDouble y = DoubleDouble.FromComponents(second, Math.ScaleB(-sign, exponent - gap - 54));
                    DoubleDouble result = DoubleDouble.Hypot(x, y);
                    AssertBound(x, y, result);
                    AssertBits(result, DoubleDouble.Hypot(-y, x));
                    AssertBits(result, DDMath.Hypot(x, y));
                }
            }
        }
        DoubleDouble maximum = new(double.MaxValue);
        AssertBound(maximum, new DoubleDouble(Math.ScaleB(1.0, 995)), DoubleDouble.Hypot(maximum, new DoubleDouble(Math.ScaleB(1.0, 995))));
    }

    [Fact]
    public void WidelySeparatedOperandsRetainRepresentableCorrections()
    {
        // sqrt(a^2+b^2) = a+b^2/(2a)+O(b^4/a^3). These exact dyadics
        // put the omitted term below half an ulp of the expected low.
        foreach ((int large, int small) in new[] { (0, -100), (1000, 450), (-500, -600) })
        {
            DoubleDouble result = DoubleDouble.Hypot(new DoubleDouble(Math.ScaleB(1.0, large)), new DoubleDouble(Math.ScaleB(1.0, small)));
            AssertBits(DoubleDouble.FromComponents(Math.ScaleB(1.0, large), Math.ScaleB(1.0, (2 * small) - large - 1)), result);
        }
    }

    [Fact]
    public void OverflowClassificationUsesTheExactInputSquares()
    {
        DoubleDouble x = DoubleDouble.FromComponents(double.MaxValue, Math.BitDecrement(Math.ScaleB(1.0, 970)));
        double boundaryY = Math.ScaleB(1.0, 971);
        // T=2^1024-2^970 is the binary64 overflow midpoint. For the lower y,
        // T²-x²-y² = 3*2^1888-5*2^1834 > 0, though the scaled root can round
        // up to T before rescaling. The adjacent boundaryY is above the threshold.
        foreach (double y in new[] { Math.BitDecrement(boundaryY), boundaryY, Math.BitIncrement(boundaryY) })
        {
            AssertOverflowBoundary(x, new DoubleDouble(y));
        }
    }

    [Fact]
    public void ExactOverflowMidpointAndAdjacentLowsAreClassified()
    {
        // 55²+48²=73² and 73 divides 2^54-1, so these scaled integer legs
        // have hypotenuse exactly T=(2^54-1)*2^970. Split integer components
        // before scaling; no DD arithmetic supplies the expected threshold.
        const long factor = ((1L << 54) - 1) / 73;
        const long first = 55 * factor;
        const long second = 48 * factor;
        DoubleDouble x = DoubleDouble.FromComponents(Math.ScaleB((double)first, 970), Math.ScaleB(first - (long)(double)first, 970));
        DoubleDouble y = DoubleDouble.FromComponents(Math.ScaleB((double)second, 970), Math.ScaleB(second - (long)(double)second, 970));
        BigInteger threshold = Units(double.MaxValue) + Units(Math.ScaleB(1.0, 970));
        BigInteger a = Units(x.High) + Units(x.Low);
        BigInteger b = Units(y.High) + Units(y.Low);
        ((a * a) + (b * b)).ShouldBe(threshold * threshold);
        foreach (double low in new[] { Math.BitDecrement(x.Low), x.Low, Math.BitIncrement(x.Low) })
        {
            AssertOverflowBoundary(DoubleDouble.FromComponents(x.High, low), y);
        }
    }

    [Fact]
    public void ThreeDimensionalSpecialValuesAndZeroReductionFollowTheContract()
    {
        double[] values = [0.0, -0.0, 1.0, -1.0, double.PositiveInfinity, double.NegativeInfinity, double.NaN];
        foreach (double x in values)
        {
            foreach (double y in values)
            {
                foreach (double z in values)
                {
                    DoubleDouble result = DoubleDouble.Hypot(new DoubleDouble(x), new DoubleDouble(y), new DoubleDouble(z));
                    if (double.IsInfinity(x) || double.IsInfinity(y) || double.IsInfinity(z))
                    {
                        AssertBits(DoubleDouble.PositiveInfinity, result);
                    }
                    else if (double.IsNaN(x) || double.IsNaN(y) || double.IsNaN(z))
                    {
                        AssertBits(DoubleDouble.NaN, result);
                    }
                    else if (x == 0.0 && y == 0.0 && z == 0.0)
                    {
                        AssertBits(DoubleDouble.Zero, result);
                    }
                    else
                    {
                        AssertBound(new DoubleDouble(x), new DoubleDouble(y), new DoubleDouble(z), result);
                    }
                    AssertBits(result, DDMath.Hypot(new DoubleDouble(x), new DoubleDouble(y), new DoubleDouble(z)));
                }
            }
        }
        DoubleDouble[] finite = [new(double.Epsilon), new(double.MaxValue), new(2.0),
            DoubleDouble.FromComponents(1.0, double.Epsilon)];
        foreach (DoubleDouble x in finite)
        {
            foreach (DoubleDouble y in finite)
            {
                DoubleDouble expected = DoubleDouble.Hypot(x, y);
                AssertBits(expected, DoubleDouble.Hypot(x, y, DoubleDouble.NegativeZero));
                AssertBits(expected, DoubleDouble.Hypot(x, DoubleDouble.Zero, y));
                AssertBits(expected, DoubleDouble.Hypot(DoubleDouble.NegativeZero, x, y));
            }
        }
    }

    [Fact]
    public void ThreeDimensionalExactBinaryNormsCoverTheExponentRange()
    {
        // 1²+2²+2²=3²; all scaled inputs and outputs here are exact dyadics.
        for (int exponent = -1074; exponent <= 1021; ++exponent)
        {
            DoubleDouble x = new(Math.ScaleB(1.0, exponent));
            DoubleDouble y = new(Math.ScaleB(2.0, exponent));
            DoubleDouble expected = new(Math.ScaleB(3.0, exponent));
            AssertBits(expected, DoubleDouble.Hypot(x, y, y));
            AssertBits(expected, DoubleDouble.Hypot(-y, x, -y));
        }
        AssertBits(new DoubleDouble(7.0), DoubleDouble.Hypot(new DoubleDouble(2.0), new DoubleDouble(3.0), new DoubleDouble(6.0)));
        AssertBits(new DoubleDouble(13.0), DoubleDouble.Hypot(new DoubleDouble(3.0), new DoubleDouble(4.0), new DoubleDouble(12.0)));
    }

    [Fact]
    public void ThreeDimensionalFiniteResultsMeetExactSquaredBound()
    {
        Random random = new(314265);
        for (int exponent = -1074; exponent <= 1023; ++exponent)
        {
            double high = Math.ScaleB(1.0 + (random.NextDouble() * 0.1), exponent);
            foreach (int gap in new[] { 0, 1, 26, 54, 55, 500, 1100 })
            {
                foreach (double sign in new[] { -1.0, 1.0 })
                {
                    DoubleDouble x = DoubleDouble.FromComponents(high, Math.ScaleB(sign, exponent - 54));
                    DoubleDouble y = DoubleDouble.FromComponents(Math.ScaleB(1.0, exponent - gap), Math.ScaleB(-sign, exponent - gap - 54));
                    DoubleDouble z = DoubleDouble.FromComponents(Math.ScaleB(0.75, exponent - gap), Math.ScaleB(sign, exponent - gap - 55));
                    DoubleDouble result = DoubleDouble.Hypot(x, y, z);
                    AssertBound(x, y, z, result);
                    AssertBits(result, DoubleDouble.Hypot(-z, x, -y));
                    AssertBits(result, DoubleDouble.Hypot(y, -z, x));
                    AssertBits(result, DDMath.Hypot(x, y, z));
                }
            }
        }
    }

    [Fact]
    public void ThreeDimensionalSparseCorrectionsAreCombinedBeforeFinalRounding()
    {
        // sqrt(a²+2b²) = a+b²/a+O(b^4/a^3). The omitted term is below
        // half an ulp of the expected low in these exact dyadic cases.
        foreach ((int large, int small) in new[] { (0, -100), (1000, 450), (-500, -600), (0, -537), (-500, -787) })
        {
            DoubleDouble x = new(Math.ScaleB(1.0, large));
            DoubleDouble y = new(Math.ScaleB(1.0, small));
            DoubleDouble expected = DoubleDouble.FromComponents(x.High, Math.ScaleB(1.0, (2 * small) - large));
            AssertThreeDimensionalPermutations(x, y, y, expected);
        }
        DoubleDouble sparse = DoubleDouble.FromComponents(Math.ScaleB(1.0, 1000), double.Epsilon);
        AssertThreeDimensionalPermutations(sparse, DoubleDouble.One, DoubleDouble.One,
            DoubleDouble.FromComponents(sparse.High, Math.ScaleB(1.0, -1000)));
    }

    [Fact]
    public void ThreeDimensionalOverflowUsesAllExactInputSquares()
    {
        // The third coordinate changes a finite two-coordinate norm to overflow.
        DoubleDouble value = new(Math.ScaleB(1.25, 1023));
        DoubleDouble.IsFinite(DoubleDouble.Hypot(value, value)).ShouldBeTrue();
        AssertThreeDimensionalOverflowBoundary(value, value, value);
        AssertThreeDimensionalOverflowBoundary(new DoubleDouble(double.MaxValue), value, DoubleDouble.One);

        DoubleDouble x = DoubleDouble.FromComponents(double.MaxValue, Math.BitDecrement(Math.ScaleB(1.0, 970)));
        double boundaryY = Math.ScaleB(1.0, 971);
        foreach (double y in new[] { Math.BitDecrement(boundaryY), boundaryY, Math.BitIncrement(boundaryY) })
        {
            AssertThreeDimensionalOverflowBoundary(x, new DoubleDouble(y), DoubleDouble.One);
        }

        // 1²+2²+2²=3², and 3 divides 2^54-1. These exact scaled integers
        // give norm T=(2^54-1)*2^970. Even an epsilon low on the first
        // coordinate decides which side of the exact overflow midpoint we are on.
        const long factor = ((1L << 54) - 1) / 3;
        DoubleDouble first = new(Math.ScaleB((double)factor, 970));
        DoubleDouble second = new(Math.ScaleB((double)(2 * factor), 970));
        BigInteger a = Units(first.High);
        BigInteger b = Units(second.High);
        BigInteger threshold = Units(double.MaxValue) + Units(Math.ScaleB(1.0, 970));
        ((a * a) + (2 * b * b)).ShouldBe(threshold * threshold);
        foreach (double low in new[] { -double.Epsilon, 0.0, double.Epsilon })
        {
            AssertThreeDimensionalOverflowBoundary(DoubleDouble.FromComponents(first.High, low), second, second);
        }
    }

    private static void AssertThreeDimensionalOverflowBoundary(DoubleDouble x, DoubleDouble y, DoubleDouble z)
    {
        BigInteger a = Units(x.High) + Units(x.Low);
        BigInteger b = Units(y.High) + Units(y.Low);
        BigInteger c = Units(z.High) + Units(z.Low);
        BigInteger threshold = Units(double.MaxValue) + Units(Math.ScaleB(1.0, 970));
        bool overflow = ((a * a) + (b * b) + (c * c)) >= threshold * threshold;
        DoubleDouble result = DoubleDouble.Hypot(x, y, z);
        DoubleDouble.IsPositiveInfinity(result).ShouldBe(overflow);
        if (!overflow)
        {
            AssertBound(x, y, z, result);
        }
        AssertThreeDimensionalPermutations(x, y, z, result);
    }

    private static void AssertThreeDimensionalPermutations(DoubleDouble x, DoubleDouble y, DoubleDouble z, DoubleDouble expected)
    {
        foreach (DoubleDouble a in new[] { x, -x })
        {
            foreach (DoubleDouble b in new[] { y, -y })
            {
                foreach (DoubleDouble c in new[] { z, -z })
                {
                    AssertBits(expected, DoubleDouble.Hypot(a, b, c));
                    AssertBits(expected, DoubleDouble.Hypot(a, c, b));
                    AssertBits(expected, DoubleDouble.Hypot(b, a, c));
                    AssertBits(expected, DoubleDouble.Hypot(b, c, a));
                    AssertBits(expected, DoubleDouble.Hypot(c, a, b));
                    AssertBits(expected, DoubleDouble.Hypot(c, b, a));
                    AssertBits(expected, DDMath.Hypot(a, b, c));
                }
            }
        }
    }

    private static void AssertOverflowBoundary(DoubleDouble x, DoubleDouble y)
    {
        BigInteger a = Units(x.High) + Units(x.Low);
        BigInteger b = Units(y.High) + Units(y.Low);
        BigInteger threshold = Units(double.MaxValue) + Units(Math.ScaleB(1.0, 970));
        bool overflow = ((a * a) + (b * b)) >= threshold * threshold;
        foreach (DoubleDouble first in new[] { x, -x })
        {
            foreach (DoubleDouble second in new[] { y, -y })
            {
                DoubleDouble result = DoubleDouble.Hypot(first, second);
                DoubleDouble.IsPositiveInfinity(result).ShouldBe(overflow);
                if (!overflow)
                {
                    AssertBound(first, second, result);
                }
                AssertBits(result, DoubleDouble.Hypot(second, first));
                AssertBits(result, DDMath.Hypot(first, second));
            }
        }
    }

    private static void AssertBound(DoubleDouble x, DoubleDouble y, DoubleDouble actual)
    {
        AssertBound(x, y, DoubleDouble.Zero, actual);
    }

    private static void AssertBound(DoubleDouble x, DoubleDouble y, DoubleDouble z, DoubleDouble actual)
    {
        DoubleDouble.IsFinite(actual).ShouldBeTrue();
        DoubleDouble.IsCanonical(actual).ShouldBeTrue();
        double.IsNegative(actual.High).ShouldBeFalse();
        // Exact integer units of 2^-1074. The squared inequalities are equivalent
        // to |actual/sqrt(x*x+y*y+z*z)-1| <= 2^-100, allowing one subnormal unit.
        BigInteger a = Units(x.High) + Units(x.Low);
        BigInteger b = Units(y.High) + Units(y.Low);
        BigInteger c = Units(z.High) + Units(z.Low);
        BigInteger r = Units(actual.High) + Units(actual.Low);
        BigInteger sum = (a * a) + (b * b) + (c * c);
        BigInteger scale = BigInteger.One << 100;
        BigInteger lower = BigInteger.Max(BigInteger.Zero, r - 1) * scale;
        BigInteger upper = (r + 1) * scale;
        ((lower * lower <= sum * (scale + 1) * (scale + 1))
            && (upper * upper >= sum * (scale - 1) * (scale - 1))).ShouldBeTrue(
                $"Hypot bound failed for ({x.High:R}, {x.Low:R}), ({y.High:R}, {y.Low:R}), ({z.High:R}, {z.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));
    }
}