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

using System.Globalization;
using Just.PreciseMath.Tests.ReferenceData;

namespace Just.PreciseMath.Tests;

public class PreciseMathLogTests
{
    public static IEnumerable<TheoryDataRow<string, double, double, string>> ReferenceCases => LogReferenceData.Cases();

    [Fact]
    public void SpecialValuesHaveExplicitCanonicalResults()
    {
        // Legacy Log(0) reached Exp(-Infinity) and threw. Domain handling must
        // precede finite range reduction, with either zero mapping to -Infinity.
        double[] values = [0.0, -0.0, 1.0, -1.0, -double.Epsilon, double.MinValue,
            double.PositiveInfinity, double.NegativeInfinity, double.NaN];
        foreach (double value in values)
        {
            DoubleDouble expected = new(Math.Log(value));
            AssertBits(DDMath.Log(new DoubleDouble(value)), expected.High, 0.0);
        }
        DoubleDouble negative = DoubleDouble.FromComponents(-1.0, Math.ScaleB(1.0, -54));
        AssertBits(DDMath.Log(negative), double.NaN, 0.0);
    }

    [Theory]
    [MemberData(nameof(ReferenceCases))]
    public void FiniteInputsMeetIndependentReferenceBound(string label, double high, double low, string reference)
    {
        // Decimal.ln at 450/650 digits, from exact component sums formed at
        // 2200 digits and checked with Fraction. No binary64 log supplies the oracle.
        DoubleDouble input = DoubleDouble.FromComponents(high, low);
        (Units(input.High) + Units(input.Low)).ShouldBe(Units(high) + Units(low));
        DoubleDouble actual = DDMath.Log(input);
        AssertReferenceBound(actual, reference, 1, $"{label}: ({high:R}, {low:R})");
        if (input == DoubleDouble.One)
        {
            AssertBits(actual, 0.0, 0.0);
        }
        else
        {
            Math.Sign(actual.High).ShouldBe(input > DoubleDouble.One ? 1 : -1);
        }
    }

    [Fact]
    public void PowersOfTwoCoverEveryFiniteBinaryInputExponent()
    {
        // ln(2^k) = k*ln(2). Multiply the independent decimal reference by k
        // with exact integers, never using the library constant or DD arithmetic.
        for (int exponent = -1074; exponent <= 1023; ++exponent)
        {
            DoubleDouble input = new(Math.ScaleB(1.0, exponent));
            AssertReferenceBound(DDMath.Log(input), LogReferenceData.Ln2, exponent, $"2^{exponent}");
        }
    }

    [Theory]
    [InlineData(1e308)]
    [InlineData(double.MaxValue)]
    [InlineData(1e-308)]
    public void ExtremeFiniteInputsRetainExtendedPrecisionWithoutDenominatorOverflow(double input)
    {
        // The legacy Log(1e308) formed u+value and overflowed. The finite
        // kernel must reduce the mantissa before constructing its denominator.
        DoubleDouble actual = DDMath.Log(new DoubleDouble(input));
        DoubleDouble.IsFinite(actual).ShouldBeTrue();
        actual.Low.ShouldNotBe(0.0);
        DoubleDouble.IsCanonical(actual).ShouldBeTrue();
        // ReferenceCases independently pins the complete result for these inputs.
    }

    [Fact]
    public void NearOneRetainsSparseCorrectionsOfEitherSign()
    {
        foreach (double sign in new[] { -1.0, 1.0 })
        {
            AssertBits(DDMath.Log(DoubleDouble.FromComponents(1.0, sign * double.Epsilon)),
                sign * double.Epsilon, 0.0);
            double delta = sign * Math.ScaleB(1.0, -500);
            // ln(1+d) = d-d^2/2+O(d^3). Here the cubic tail is below half
            // the minimum subnormal, so these two binary components are exact.
            AssertBits(DDMath.Log(DoubleDouble.FromComponents(1.0, delta)),
                delta, -Math.ScaleB(1.0, -1001));
        }
    }

    [Fact]
    public void ArbitraryPairsAreNormalizedByThePublicFactoryBeforeLog()
    {
        DoubleDouble positive = DoubleDouble.FromComponents(0.0, 2.0);
        AssertReferenceBound(DDMath.Log(positive), LogReferenceData.Ln2, 1, "zero-high public construction");
        AssertBits(DDMath.Log(DoubleDouble.FromComponents(0.0, -2.0)), double.NaN, 0.0);
    }

    private static void AssertReferenceBound(DoubleDouble actual, string reference, int factor, string context)
    {
        DoubleDouble.IsFinite(actual).ShouldBeTrue(context);
        DoubleDouble.IsCanonical(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 = factor * 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 actualUnits = Units(actual.High) + Units(actual.Low);
        BigInteger error = BigInteger.Abs((actualUnits * denominator) - (numerator << 1074));
        BigInteger magnitude = BigInteger.Abs(numerator);
        // 2^-100 relative + one minimum subnormal, with independent reference
        // uncertainty <2^-350 relative explicitly included. No component collapse.
        BigInteger bound = (magnitude << 1324) + (denominator << 350) + (magnitude << 1074);
        ((error << 350) <= bound).ShouldBeTrue(
            $"Log 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();
    }
}