namespace Just.PreciseMath;

public static partial class DDMath
{
    /// <summary>Returns e raised to the specified double-double value.</summary>
    /// <remarks>
    /// Uses binary range reduction and a [12/12] Padé approximation. Results are
    /// approximate, not guaranteed correctly rounded; tests check 2^-100 relative
    /// error plus one minimum binary64 subnormal against high-precision references.
    /// Precision decreases near underflow. NaN returns canonical NaN; positive
    /// infinity returns positive infinity, negative infinity returns positive zero,
    /// and either zero returns one. Both input components affect range boundaries.
    /// </remarks>
    [Pure]
    public static DoubleDouble Exp(DoubleDouble value)
    {
        if (double.IsNaN(value.High))
        {
            return DoubleDouble.NaN;
        }
        // These deliberately loose bounds only reject inputs safely outside the
        // result range, even with a normalized low of either sign. They also keep
        // infinities and huge finite inputs out of the floating-to-int conversion.
        if (value.High > 710.0)
        {
            return new DoubleDouble(double.PositiveInfinity);
        }
        if (value.High < -746.0)
        {
            return DoubleDouble.Zero;
        }
        if (value.High == 0.0)
        {
            return DoubleDouble.One;
        }

        // The binary64 estimate need not choose the nearest k on a tie: either
        // neighbor leaves |r| < 0.347. k is bounded by [-1076, 1024]. Subtract
        // separate products to retain cancellation residuals. A third ln(2)
        // component prevents k times the DD constant error from dominating r.
        // Generated independently at 180/260 decimal digits by generate_exp.py.
        const double ln2Tail = 5.707708438416212e-34;
        int exponent = (int)Math.Round(value.High / DoubleDouble.Ln2.High);
        DoubleDouble reduced = value - (new DoubleDouble(DoubleDouble.Ln2.High) * exponent);
        reduced -= new DoubleDouble(DoubleDouble.Ln2.Low) * exponent;
        reduced -= new DoubleDouble(ln2Tail) * exponent;

        // P(x)/P(-x), with exact binary64 integer coefficients. For n=12,
        // c_k = (24-k)!*12! / (24!*k!*(12-k)!), scaled by 1/c_12.
        // This retains the legacy approximation but not its unsafe E^k scaling.
        // Horner accumulators remain normal and finite on the reduced interval.
        ReadOnlySpan<double> coefficients = [156.0, 12012.0, 600600.0, 21621600.0,
            588107520.0, 12350257920.0, 201132771840.0, 2514159648000.0,
            23465490048000.0, 154872234316800.0, 647647525324800.0, 1295295050649600.0];
        DoubleDouble numerator = DoubleDouble.One;
        DoubleDouble denominator = DoubleDouble.One;
        foreach (double coefficient in coefficients)
        {
            numerator = (numerator * reduced) + coefficient;
            denominator = (denominator * -reduced) + coefficient;
        }
        return ScalePowerOfTwo(numerator / denominator, exponent);
    }
}