addoed exponential functions

0-source/Just.PreciseMath/DDMath.Exp.cs

@@ -2,63 +2,45 @@ 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>
+    /// <inheritdoc cref="DoubleDouble.Exp"/>
     [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;
-        }
+        return DoubleDouble.Exp(value);
+    }
 
-        // 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;
+    /// <inheritdoc cref="DoubleDouble.Exp2"/>
+    [Pure]
+    public static DoubleDouble Exp2(DoubleDouble value)
+    {
+        return DoubleDouble.Exp2(value);
+    }
 
-        // 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);
+    /// <inheritdoc cref="DoubleDouble.Exp10"/>
+    [Pure]
+    public static DoubleDouble Exp10(DoubleDouble value)
+    {
+        return DoubleDouble.Exp10(value);
+    }
+
+    /// <inheritdoc cref="DoubleDouble.ExpM1"/>
+    [Pure]
+    public static DoubleDouble ExpM1(DoubleDouble value)
+    {
+        return DoubleDouble.ExpM1(value);
+    }
+
+    /// <inheritdoc cref="DoubleDouble.Exp2M1"/>
+    [Pure]
+    public static DoubleDouble Exp2M1(DoubleDouble value)
+    {
+        return DoubleDouble.Exp2M1(value);
+    }
+
+    /// <inheritdoc cref="DoubleDouble.Exp10M1"/>
+    [Pure]
+    public static DoubleDouble Exp10M1(DoubleDouble value)
+    {
+        return DoubleDouble.Exp10M1(value);
     }
 }

0-source/Just.PreciseMath/DDMath.Pow.cs

@@ -70,7 +70,7 @@ public static partial class DDMath
             result = 1.0 / result;
             resultExponent = -resultExponent;
         }
-        return ScalePowerOfTwo(negative ? -result : result, resultExponent);
+        return PreciseMathHelper.ScalePowerOfTwo(negative ? -result : result, resultExponent);
     }
 
     /// <summary>Raises a double-double value to a binary64 power.</summary>

0-source/Just.PreciseMath/DDMath.cs

@@ -89,29 +89,4 @@ public static partial class DDMath
     {
         return DoubleDouble.InvSqrt(value);
     }
-
-    // For finite, nonzero normalized significands and exponents bounded by the
-    // integer-power domain (|exponent| < 2^42). Keep the exponent separate until
-    // the final result so an intermediate cannot overflow before reciprocation.
-    private static DoubleDouble ScalePowerOfTwo(DoubleDouble value, long exponent)
-    {
-        long resultExponent = Math.ILogB(value.High) + exponent;
-        if (resultExponent > 1024)
-        {
-            return new DoubleDouble(Math.CopySign(double.PositiveInfinity, value.High));
-        }
-        if (resultExponent < -1075)
-        {
-            return new DoubleDouble(Math.CopySign(0.0, value.High));
-        }
-        if (resultExponent >= -969 && resultExponent <= 1022)
-        {
-            // Normal high and room for a dense low; a sparse low may underflow
-            // by at most half a minimum subnormal. Normalize its signed zero.
-            return DoubleDouble.FromComponents(Math.ScaleB(value.High, (int)exponent),
-                Math.ScaleB(value.Low, (int)exponent));
-        }
-
-        return PreciseMathHelper.ScalePowerOfTwoBoundary(value, exponent);
-    }
 }

0-source/Just.PreciseMath/DoubleDouble.ExponentialFunctions.cs

@@ -0,0 +1,237 @@
+namespace Just.PreciseMath;
+
+public readonly partial struct DoubleDouble : IExponentialFunctions<DoubleDouble>
+{
+    /// <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;
+
+        return PreciseMathHelper.ScalePowerOfTwo(ExpReduced(reduced), exponent);
+    }
+
+    /// <summary>Returns two raised to the specified double-double value.</summary>
+    /// <remarks>
+    /// Reduces in base two before evaluating a bounded natural exponential, avoiding
+    /// amplification of ln(2) rounding error by a large input. Integer powers in the
+    /// finite binary64 range are exact. Special values and approximate accuracy
+    /// follow <see cref="Exp"/>; both components determine range boundaries.
+    /// </remarks>
+    [Pure]
+    public static DoubleDouble Exp2(DoubleDouble value)
+    {
+        if (IsNaN(value))
+        {
+            return NaN;
+        }
+        if (value.High > 1024.0)
+        {
+            return PositiveInfinity;
+        }
+        if (value <= new DoubleDouble(-1075.0))
+        {
+            return Zero;
+        }
+        if (value < new DoubleDouble(-1074.0))
+        {
+            // Include sparse lows just above the exact half-subnormal tie;
+            // rounding exp(r) to one must not erase which side the input is on.
+            return Epsilon;
+        }
+        int exponent = (int)Math.Round(value.High);
+        DoubleDouble fraction = value - exponent;
+        if (fraction.High != 0.0 && Math.Abs(fraction.High) < Math.ScaleB(1.0, -500) && exponent <= 1023)
+        {
+            // 2^(k+d) = 2^k + 2^k*d*ln(2) + O(2^k*d²). Scale d to
+            // a bounded mantissa before multiplying: d*ln(2) may otherwise
+            // round in the subnormal range before 2^k restores the correction.
+            // The omitted term is < 2^-500 relative to the correction itself.
+            int adjustment = Math.ILogB(fraction.High);
+            DoubleDouble mantissa = FromComponents(Math.ScaleB(fraction.High, -adjustment),
+                Math.ScaleB(fraction.Low, -adjustment));
+            DoubleDouble correction = PreciseMathHelper.ScalePowerOfTwo(mantissa * Ln2, exponent + adjustment);
+            return new DoubleDouble(Math.ScaleB(1.0, exponent)) + correction;
+        }
+        DoubleDouble reduced = fraction * Ln2;
+        return PreciseMathHelper.ScalePowerOfTwo(ExpReduced(reduced), exponent);
+    }
+
+    /// <summary>Returns ten raised to the specified double-double value.</summary>
+    /// <remarks>
+    /// Subtracts three split log10(2) products before conversion to a bounded
+    /// natural exponent. Special values and approximate accuracy follow
+    /// <see cref="Exp"/>; precision decreases near underflow.
+    /// </remarks>
+    [Pure]
+    public static DoubleDouble Exp10(DoubleDouble value)
+    {
+        if (IsNaN(value))
+        {
+            return NaN;
+        }
+        // Loose guards include infinities and keep the exponent conversion bounded.
+        if (value.High > 309.0)
+        {
+            return PositiveInfinity;
+        }
+        if (value.High < -324.0)
+        {
+            return Zero;
+        }
+        if (IsZero(value))
+        {
+            return One;
+        }
+        // Independently split ln(2)/ln(10) at 180 and 260 decimal digits.
+        const double log10Of2Tail = 5.471948402314639e-35;
+        int exponent = (int)Math.Round(value.High * Log2Of10.High);
+        DoubleDouble reduced = value - (new DoubleDouble(Log10Of2.High) * exponent);
+        reduced -= new DoubleDouble(Log10Of2.Low) * exponent;
+        reduced -= new DoubleDouble(log10Of2Tail) * exponent;
+        return PreciseMathHelper.ScalePowerOfTwo(ExpReduced(reduced * Ln10), exponent);
+    }
+
+    /// <summary>Returns e raised to the specified value, minus one.</summary>
+    /// <remarks>
+    /// Uses a direct series near zero rather than subtracting one from a rounded
+    /// exponential, preserving tiny results. Signed zeros are preserved; negative
+    /// infinity returns negative one, positive infinity returns positive infinity,
+    /// and NaN returns canonical NaN. The tested approximate error bound is
+    /// 2^-100 relative to exp(value)-1 plus one minimum binary64 subnormal;
+    /// results are not guaranteed correctly rounded.
+    /// </remarks>
+    [Pure]
+    public static DoubleDouble ExpM1(DoubleDouble value)
+    {
+        if (IsZero(value))
+        {
+            return value;
+        }
+        if (Math.Abs(value.High) <= 0.5)
+        {
+            return ExpM1Small(value);
+        }
+        return Exp(value) - 1.0;
+    }
+
+    /// <summary>Returns two raised to the specified value, minus one.</summary>
+    /// <remarks>
+    /// Uses a cancellation-safe series near zero. Special values, signed zeros,
+    /// and approximate accuracy follow <see cref="ExpM1"/>, with the error bound
+    /// relative to 2^value-1. Other inputs use the base-two range reduction of
+    /// <see cref="Exp2"/> before subtracting one.
+    /// </remarks>
+    [Pure]
+    public static DoubleDouble Exp2M1(DoubleDouble value)
+    {
+        if (IsZero(value))
+        {
+            return value;
+        }
+        if (Math.Abs(value.High) <= 0.5)
+        {
+            return ExpM1Small(value * Ln2);
+        }
+        return Exp2(value) - 1.0;
+    }
+
+    /// <summary>Returns ten raised to the specified value, minus one.</summary>
+    /// <remarks>
+    /// Uses a cancellation-safe series near zero. Special values, signed zeros,
+    /// and approximate accuracy follow <see cref="ExpM1"/>, with the error bound
+    /// relative to 10^value-1. Other inputs use the base-ten range reduction of
+    /// <see cref="Exp10"/> before subtracting one.
+    /// </remarks>
+    [Pure]
+    public static DoubleDouble Exp10M1(DoubleDouble value)
+    {
+        if (IsZero(value))
+        {
+            return value;
+        }
+        if (Math.Abs(value.High) <= 0.125)
+        {
+            return ExpM1Small(value * Ln10);
+        }
+        return Exp10(value) - 1.0;
+    }
+
+    // Requires a finite normalized argument with |value| <= 0.5 plus rounding.
+    private static DoubleDouble ExpM1Small(DoubleDouble value)
+    {
+        if (Math.Abs(value.High) <= Math.ScaleB(1.0, -54))
+        {
+            // expm1(x) = x + x²/2 + O(x³): omitted relative error < 2^-110.
+            // Retain x even when its square underflows; never halve x first.
+            return value + ((value * value) * 0.5);
+        }
+        DoubleDouble term = value;
+        DoubleDouble sum = value;
+        for (int denominator = 2; denominator <= 32; ++denominator)
+        {
+            term = (term * value) / denominator;
+            sum += term;
+        }
+        // Relative truncation error <= 2*(0.5)^32/33! < 5.4e-47;
+        // double-double rounding dominates, including for negative arguments.
+        return sum;
+    }
+
+    // Requires a finite normalized argument with |reduced| < 0.347.
+    private static DoubleDouble ExpReduced(DoubleDouble reduced)
+    {
+        // 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 numerator / denominator;
+    }
+}

0-source/Just.PreciseMath/PreciseMathHelper.cs

@@ -212,6 +212,31 @@ internal static class PreciseMathHelper
         return ArithmeticFromRatio(ArithmeticUnits(left), ArithmeticUnits(right));
     }
 
+    // For finite, nonzero normalized significands and exponents bounded by the
+    // integer-power domain (|exponent| < 2^42). Keep the exponent separate until
+    // the final result so an intermediate cannot overflow before reciprocation.
+    internal static DoubleDouble ScalePowerOfTwo(DoubleDouble value, long exponent)
+    {
+        long resultExponent = Math.ILogB(value.High) + exponent;
+        if (resultExponent > 1024)
+        {
+            return new DoubleDouble(Math.CopySign(double.PositiveInfinity, value.High));
+        }
+        if (resultExponent < -1075)
+        {
+            return new DoubleDouble(Math.CopySign(0.0, value.High));
+        }
+        if (resultExponent >= -969 && resultExponent <= 1022)
+        {
+            // Normal high and room for a dense low; a sparse low may underflow
+            // by at most half a minimum subnormal. Normalize its signed zero.
+            return DoubleDouble.FromComponents(Math.ScaleB(value.High, (int)exponent),
+                Math.ScaleB(value.Low, (int)exponent));
+        }
+
+        return ScalePowerOfTwoBoundary(value, exponent);
+    }
+
     // The caller supplies a finite, nonzero normalized value and has bounded
     // ILogB(value.High) + exponent to [-1075, 1024], keeping both shifts small.
     // Isolate all allocating setup from ordinary scaling, including conversion