diff --git a/0-source/Just.PreciseMath/DDMath.Exp.cs b/0-source/Just.PreciseMath/DDMath.Exp.cs index 9013076..e06c287 100644 --- a/0-source/Just.PreciseMath/DDMath.Exp.cs +++ b/0-source/Just.PreciseMath/DDMath.Exp.cs @@ -2,63 +2,45 @@ namespace Just.PreciseMath; public static partial class DDMath { - /// Returns e raised to the specified double-double value. - /// - /// 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. - /// + /// [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; + /// + [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 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); + /// + [Pure] + public static DoubleDouble Exp10(DoubleDouble value) + { + return DoubleDouble.Exp10(value); + } + + /// + [Pure] + public static DoubleDouble ExpM1(DoubleDouble value) + { + return DoubleDouble.ExpM1(value); + } + + /// + [Pure] + public static DoubleDouble Exp2M1(DoubleDouble value) + { + return DoubleDouble.Exp2M1(value); + } + + /// + [Pure] + public static DoubleDouble Exp10M1(DoubleDouble value) + { + return DoubleDouble.Exp10M1(value); } } diff --git a/0-source/Just.PreciseMath/DDMath.Pow.cs b/0-source/Just.PreciseMath/DDMath.Pow.cs index dc87fde..7ddfbdf 100644 --- a/0-source/Just.PreciseMath/DDMath.Pow.cs +++ b/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); } /// Raises a double-double value to a binary64 power. diff --git a/0-source/Just.PreciseMath/DDMath.cs b/0-source/Just.PreciseMath/DDMath.cs index fa58c36..8b0abc5 100644 --- a/0-source/Just.PreciseMath/DDMath.cs +++ b/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); - } } diff --git a/0-source/Just.PreciseMath/DoubleDouble.ExponentialFunctions.cs b/0-source/Just.PreciseMath/DoubleDouble.ExponentialFunctions.cs new file mode 100644 index 0000000..78593aa --- /dev/null +++ b/0-source/Just.PreciseMath/DoubleDouble.ExponentialFunctions.cs @@ -0,0 +1,237 @@ +namespace Just.PreciseMath; + +public readonly partial struct DoubleDouble : IExponentialFunctions +{ + /// Returns e raised to the specified double-double value. + /// + /// 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. + /// + [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); + } + + /// Returns two raised to the specified double-double value. + /// + /// 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 ; both components determine range boundaries. + /// + [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); + } + + /// Returns ten raised to the specified double-double value. + /// + /// Subtracts three split log10(2) products before conversion to a bounded + /// natural exponent. Special values and approximate accuracy follow + /// ; precision decreases near underflow. + /// + [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); + } + + /// Returns e raised to the specified value, minus one. + /// + /// 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. + /// + [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; + } + + /// Returns two raised to the specified value, minus one. + /// + /// Uses a cancellation-safe series near zero. Special values, signed zeros, + /// and approximate accuracy follow , with the error bound + /// relative to 2^value-1. Other inputs use the base-two range reduction of + /// before subtracting one. + /// + [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; + } + + /// Returns ten raised to the specified value, minus one. + /// + /// Uses a cancellation-safe series near zero. Special values, signed zeros, + /// and approximate accuracy follow , with the error bound + /// relative to 10^value-1. Other inputs use the base-ten range reduction of + /// before subtracting one. + /// + [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 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; + } +} diff --git a/0-source/Just.PreciseMath/PreciseMathHelper.cs b/0-source/Just.PreciseMath/PreciseMathHelper.cs index efd7927..28b4092 100644 --- a/0-source/Just.PreciseMath/PreciseMathHelper.cs +++ b/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 diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleExponentialFunctionsTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleExponentialFunctionsTests.cs new file mode 100644 index 0000000..fe45528 --- /dev/null +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleExponentialFunctionsTests.cs @@ -0,0 +1,227 @@ +using System.Globalization; +using System.Numerics; +using Just.PreciseMath.Tests.ReferenceData; +using Shouldly; +using Xunit; + +namespace Just.PreciseMath.Tests; + +public class DoubleDoubleExponentialFunctionsTests +{ + public static IEnumerable> ReferenceCases => + ExponentialFunctionsReferenceData.Cases(); + + [Theory] + [MemberData(nameof(ReferenceCases))] + public void FiniteInputsMeetIndependentReferenceBound(string operation, double high, double low, + string reference, bool overflow, bool underflow) + { + // Exact binary64 sums and Decimal ln/exp/expm1 at 450/650 digits; + // integer base-two/base-ten powers are independently rational-checked. + DoubleDouble input = DoubleDouble.FromComponents(high, low); + (Units(input.High) + Units(input.Low)).ShouldBe(Units(high) + Units(low)); + DoubleDouble actual = Evaluate(operation, input); + DoubleDouble.IsCanonical(actual).ShouldBeTrue(); + AssertBits(actual, EvaluateGeneric(operation, input)); + AssertBits(actual, EvaluateFacade(operation, input)); + if (overflow) + { + AssertBits(DoubleDouble.PositiveInfinity, actual); + return; + } + if (underflow) + { + AssertBits(new DoubleDouble(reference[0] == '-' ? -0.0 : 0.0), actual); + return; + } + DoubleDouble.IsFinite(actual).ShouldBeTrue(); + DoubleDouble.IsZero(actual).ShouldBeFalse(); + string[] parts = reference.Split('e'); + int point = parts[0].IndexOf('.', StringComparison.Ordinal); + int decimals = point < 0 ? 0 : parts[0].Length - point - 1; + BigInteger numerator = 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 + epsilon absolute + 2^-350 relative reference allowance. + BigInteger bound = (magnitude << 1324) + (denominator << 350) + (magnitude << 1074); + ((error << 350) <= bound).ShouldBeTrue( + $"{operation} bound failed for ({high:R}, {low:R}): ({actual.High:R}, {actual.Low:R})"); + } + + [Fact] + public void NaturalExponentialIsAvailableOnTheTypeAndMatchesFacadeBits() + { + DoubleDouble[] values = [DoubleDouble.Zero, DoubleDouble.NegativeZero, + DoubleDouble.NaN, DoubleDouble.PositiveInfinity, DoubleDouble.NegativeInfinity, + new(-1.0), new(double.Epsilon), new(1.0), new(709.5), new(-740.0), + new(double.MaxValue), DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -1000))]; + foreach (DoubleDouble value in values) + { + AssertBits(DDMath.Exp(value), DoubleDouble.Exp(value)); + } + } + + [Fact] + public void BinaryAndDecimalExponentialsHaveExactElementaryValues() + { + for (int exponent = -1074; exponent <= 1023; ++exponent) + { + AssertBits(new DoubleDouble(Math.ScaleB(1.0, exponent)), DoubleDouble.Exp2(new DoubleDouble(exponent))); + } + AssertBits(DoubleDouble.Zero, DoubleDouble.Exp2(new DoubleDouble(-1075.0))); + AssertBits(DoubleDouble.Epsilon, DoubleDouble.Exp2(DoubleDouble.FromComponents(-1075.0, double.Epsilon))); + AssertBits(DoubleDouble.Zero, DoubleDouble.Exp2(DoubleDouble.FromComponents(-1075.0, -double.Epsilon))); + AssertBits(DoubleDouble.PositiveInfinity, DoubleDouble.Exp2(new DoubleDouble(1024.0))); + DoubleDouble ten = DoubleDouble.Exp10(DoubleDouble.One); + ten.High.ShouldBe(10.0); + Math.Abs(ten.Low).ShouldBeLessThan(1e-30); + AssertBits(DoubleDouble.One, DoubleDouble.Exp10(DoubleDouble.Zero)); + } + + [Fact] + public void MinusOneFunctionsRetainTinyResultsAndSignedZeros() + { + double tiny = Math.ScaleB(1.0, -100); + // exp(±x)-1 = ±x + x²/2 + O(x³). At x=2^-100 the tail is + // below half an ulp of the low, so these particular splits are exact. + AssertBits(DoubleDouble.FromComponents(tiny, Math.ScaleB(1.0, -201)), DoubleDouble.ExpM1(new DoubleDouble(tiny))); + AssertBits(DoubleDouble.FromComponents(-tiny, Math.ScaleB(1.0, -201)), DoubleDouble.ExpM1(new DoubleDouble(-tiny))); + foreach (double sign in new double[] { -1.0, 1.0 }) + { + DoubleDouble input = new(sign * double.Epsilon); + AssertBits(input, DoubleDouble.ExpM1(input)); + AssertBits(input, DoubleDouble.Exp2M1(input)); + AssertBits(new DoubleDouble(sign * (2.0 * double.Epsilon)), DoubleDouble.Exp10M1(input)); + } + } + + [Theory] + [InlineData(1, double.Epsilon)] + [InlineData(500, 1.1210060331144859e-173)] + [InlineData(1000, 3.6694906201918696e-23)] + public void BinaryExponentialScalesSparseCorrectionsBeforeRounding(int exponent, double expectedLow) + { + // For delta=±2^-1074, 2^(n+delta)=2^n+delta*ln(2)*2^n+O(2^(n-2148)). + // Decimal ln(2)*2^(n-1074) at 450/650 digits gives the literals above. + // The quadratic term is far below their binary64 rounding intervals. + foreach (double sign in new double[] { -1.0, 1.0 }) + { + DoubleDouble input = DoubleDouble.FromComponents(exponent, sign * double.Epsilon); + DoubleDouble expected = DoubleDouble.FromComponents(Math.ScaleB(1.0, exponent), sign * expectedLow); + AssertBits(expected, DoubleDouble.Exp2(input)); + AssertBits(expected, DDMath.Exp2(input)); + } + } + + [Theory] + [InlineData("Exp")] + [InlineData("Exp2")] + [InlineData("Exp10")] + [InlineData("ExpM1")] + [InlineData("Exp2M1")] + [InlineData("Exp10M1")] + public void CompleteContractHasGenericDispatchAndThinFacade(string operation) + { + bool minusOne = operation.EndsWith("M1", StringComparison.Ordinal); + DoubleDouble[] values = [DoubleDouble.Zero, DoubleDouble.NegativeZero, + DoubleDouble.NaN, DoubleDouble.PositiveInfinity, DoubleDouble.NegativeInfinity, + new(double.MaxValue), new(double.MinValue), new(double.Epsilon), + new(0.125), new(-0.125), new(1.0), new(-1.0), + DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -1000))]; + foreach (DoubleDouble value in values) + { + DoubleDouble actual = Evaluate(operation, value); + DoubleDouble.IsCanonical(actual).ShouldBeTrue(); + AssertBits(actual, EvaluateGeneric(operation, value)); + AssertBits(actual, EvaluateFacade(operation, value)); + if (DoubleDouble.IsZero(value)) + { + AssertBits(minusOne ? value : DoubleDouble.One, actual); + } + else if (DoubleDouble.IsNaN(value)) + { + AssertBits(DoubleDouble.NaN, actual); + } + else if (value.High == double.PositiveInfinity || value.High == double.MaxValue) + { + AssertBits(DoubleDouble.PositiveInfinity, actual); + } + else if (value.High == double.NegativeInfinity || value.High == double.MinValue) + { + AssertBits(minusOne ? DoubleDouble.NegativeOne : DoubleDouble.Zero, actual); + } + } + } + + private static DoubleDouble Evaluate(string operation, DoubleDouble value) + { + return operation switch + { + "Exp" => DoubleDouble.Exp(value), + "Exp2" => DoubleDouble.Exp2(value), + "Exp10" => DoubleDouble.Exp10(value), + "ExpM1" => DoubleDouble.ExpM1(value), + "Exp2M1" => DoubleDouble.Exp2M1(value), + "Exp10M1" => DoubleDouble.Exp10M1(value), + _ => throw new ArgumentOutOfRangeException(nameof(operation)) + }; + } + + private static T EvaluateGeneric(string operation, T value) where T : IExponentialFunctions + { + return operation switch + { + "Exp" => T.Exp(value), + "Exp2" => T.Exp2(value), + "Exp10" => T.Exp10(value), + "ExpM1" => T.ExpM1(value), + "Exp2M1" => T.Exp2M1(value), + "Exp10M1" => T.Exp10M1(value), + _ => throw new ArgumentOutOfRangeException(nameof(operation)) + }; + } + + private static DoubleDouble EvaluateFacade(string operation, DoubleDouble value) + { + return operation switch + { + "Exp" => DDMath.Exp(value), + "Exp2" => DDMath.Exp2(value), + "Exp10" => DDMath.Exp10(value), + "ExpM1" => DDMath.ExpM1(value), + "Exp2M1" => DDMath.Exp2M1(value), + "Exp10M1" => DDMath.Exp10M1(value), + _ => throw new ArgumentOutOfRangeException(nameof(operation)) + }; + } + + 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)); + } +} diff --git a/1-tests/Just.PreciseMath.Tests/ReferenceData/ExponentialFunctionsReferenceData.cs b/1-tests/Just.PreciseMath.Tests/ReferenceData/ExponentialFunctionsReferenceData.cs new file mode 100644 index 0000000..1de9ca4 --- /dev/null +++ b/1-tests/Just.PreciseMath.Tests/ReferenceData/ExponentialFunctionsReferenceData.cs @@ -0,0 +1,427 @@ +// Generated by generate_exponential_functions.py; do not hand-edit. +// Exact binary64 sums verified with Fraction at 2200 digits; Decimal ln/exp and +// cancellation-safe expm1 at 450/650 digits; 120-digit references agree. +// Integer base-2/base-10 cases additionally use exact rational powers. +// Row: operation, high, low, reference, overflow, underflow. +// Overflow: exact result >= 2^1024 - 2^970 (binary64 overflow midpoint). +// Underflow: |exact result| <= 2^-1075, NOT exp(x) <= 2^-1075 for M1. +// Compare exact component sums with relative tolerance 2^-100 + double.Epsilon +// absolute; reference rounding uncertainty is bounded by 2^-350 relative. +// Tiny corrections to 1/-1 below 120 digits are not component-retention oracles. +// Signed zeros, infinities, NaNs, and enormous arguments are tested separately. +using Xunit; + +namespace Just.PreciseMath.Tests.ReferenceData; + +internal static class ExponentialFunctionsReferenceData +{ + internal static IEnumerable> Cases() + { + // Exp2: ordinary + yield return new("Exp2", -10.25, 0.0, "8.21187905521205608428833472883998920937533459332797373841041099584887456009670156068675785350871617356251536858905356524e-4", false, false); + yield return new("Exp2", -2.0, 0.0, "2.5e-1", false, false); + yield return new("Exp2", -1.0, 0.0, "5e-1", false, false); + yield return new("Exp2", -0.125, 0.0, "9.17004043204671231743541594794144428038655166436839749791662069353238831122347362841465489380660313672341059668949459679e-1", false, false); + yield return new("Exp2", 0.125, 0.0, "1.09050773266525765920701065576070797899270271854006712178566764768330053084884184033821114049420311989145161926291809001e+0", false, false); + yield return new("Exp2", 0.75, 0.0, "1.68179283050742908606225095246642979008006852471356902162645217194984950990780447962864800839858507234560314748703817016e+0", false, false); + yield return new("Exp2", 1.0, 0.0, "2e+0", false, false); + yield return new("Exp2", 3.25, 0.0, "9.51365692002176853373999976448380732234377673971053930415201779775573334581533727896625075630510138972829915819815370549e+0", false, false); + // Exp2: dense-and-sparse-low + yield return new("Exp2", -1.0, 5.551115123125783e-17, "5.00000000000000019238698982791550072604160854284538415449960340473288044078594786193838693052953287736809338580638344518e-1", false, false); + yield return new("Exp2", -1.0, -5.551115123125783e-17, "4.99999999999999980761301017208450667650916246644680054760727103703961478639310202113257847265570935275551812571858143798e-1", false, false); + yield return new("Exp2", -1.0, 5e-324, "5.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-1", false, false); + yield return new("Exp2", -1.0, -5e-324, "5.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-1", false, false); + yield return new("Exp2", 1.0, 5.551115123125783e-17, "2.00000000000000007695479593116620029041664341713815366179984136189315217631437914477535477221181315094723735432255337807e+0", false, false); + yield return new("Exp2", 1.0, -5.551115123125783e-17, "1.99999999999999992304520406883380267060366498657872021904290841481584591455724080845303138906228374110220725028743257519e+0", false, false); + yield return new("Exp2", 1.0, 5e-324, "2.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp2", 1.0, -5e-324, "2.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + // Exp2: tiny + yield return new("Exp2", -1e-20, 0.0, "9.99999999999999999993068528194400547286020253980627191054756284100556423022379868463187217573182659816366797509203417810e-1", false, false); + yield return new("Exp2", 1e-20, 0.0, "1.00000000000000000000693147180559945271402779132076462908244017094540055080820895837340445945043791570869517416287708474e+0", false, false); + yield return new("Exp2", -1e-100, 0.0, "9.99999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999930685281944005467673e-1", false, false); + yield return new("Exp2", 1e-100, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000006931471805599453233e+0", false, false); + yield return new("Exp2", -1e-300, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp2", 1e-300, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp2", -2.2250738585072014e-308, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp2", 2.2250738585072014e-308, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp2", -1e-323, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp2", 1e-323, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp2", -5e-324, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp2", 5e-324, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + // Exp2: absolute-log-half-switch + yield return new("Exp2", -0.7213475204444817, -1.0177636870465518e-17, "6.06530659712633423603799534991179582561900207458797347807571359489015133003532873390753100210381268501490009310808759694e-1", false, false); + yield return new("Exp2", -0.7213475204444817, -1.0177636870465517e-17, "6.06530659712633423603799534991180230313795998115094343092035618278013018112551514441716068365329629038114330858826160472e-1", false, false); + yield return new("Exp2", -0.7213475204444817, -1.0177636870465515e-17, "6.06530659712633423603799534991180878065691788771391338376499877067702677844635677291755121257605846666377380013169852043e-1", false, false); + yield return new("Exp2", 0.7213475204444817, 1.0177636870465515e-17, "1.64872127070012814684865078781416241740668824728100260241580595072860002906738406941721084254203061497902616864239766825e+0", false, false); + yield return new("Exp2", 0.7213475204444817, 1.0177636870465517e-17, "1.64872127070012814684865078781416417817889592491912837676279009285057271715417266380529283718974446456952552825625487127e+0", false, false); + yield return new("Exp2", 0.7213475204444817, 1.0177636870465518e-17, "1.64872127070012814684865078781416593895110360255725415110977423497442584362822937159780260677266950493781041819896683958e+0", false, false); + // Exp2: sparse-correction-switch + yield return new("Exp2", 1000.0, -3.0549363634996054e-151, "1.07150860718626732094842504906000181056140481170553360744375038837035105112493612249319837881569585812759467291755314683e+301", false, false); + yield return new("Exp2", 1000.0, -3.054936363499605e-151, "1.07150860718626732094842504906000181056140481170553360744375038837035105112493612249319837881569585812759467291755314683e+301", false, false); + yield return new("Exp2", 1000.0, -3.0549363634996043e-151, "1.07150860718626732094842504906000181056140481170553360744375038837035105112493612249319837881569585812759467291755314683e+301", false, false); + yield return new("Exp2", 1000.0, 3.0549363634996043e-151, "1.07150860718626732094842504906000181056140481170553360744375038837035105112493612249319837881569585812759467291755314683e+301", false, false); + yield return new("Exp2", 1000.0, 3.054936363499605e-151, "1.07150860718626732094842504906000181056140481170553360744375038837035105112493612249319837881569585812759467291755314683e+301", false, false); + yield return new("Exp2", 1000.0, 3.0549363634996054e-151, "1.07150860718626732094842504906000181056140481170553360744375038837035105112493612249319837881569585812759467291755314683e+301", false, false); + // Exp2: random-moderate + yield return new("Exp2", 8.78452800194328, -4.440892098500626e-16, "4.40967341476995196084861795619177036986833446346865068997727335172090686266840323292806805468322290875886544488139632805e+2", false, false); + yield return new("Exp2", 2.7202982684631807, 1.1102230246251565e-16, "6.59009045860470968459866968889737028605350979742874890837108301124392819835145507472654892784605100090392179416858372658e+0", false, false); + yield return new("Exp2", -3.4144589655515833, 1.1102230246251565e-16, "9.37876021346977504985476134528503635574385608741052067745862593845534384960660409900825318502561533036509957790741090194e-2", false, false); + yield return new("Exp2", 3.67643978981615, 1.1102230246251565e-16, "1.27855275975436677510222614327660438420381423323350572823327975657260641462810932407195792018571153800381071193261125222e+1", false, false); + yield return new("Exp2", 7.559099749783131, -2.220446049250313e-16, "1.88588742600978100098560391503204728439903232936850014376858946537478147694222110588671754338725081209043377353414161445e+2", false, false); + yield return new("Exp2", 1.1272538511782404, 5.551115123125783e-17, "2.18442542093619461179505632322739881400577647657120007022725830276626252983657181904591357904539459020575218758691298325e+0", false, false); + // Exp2: random-range + yield return new("Exp2", -816.2214079897565, -2.842170943040401e-14, "1.96278663442385783610897157419095408652661650100352671243288814392820295492839788615258722894068315168311070489220014454e-246", false, false); + yield return new("Exp2", -678.2241907878085, -2.842170943040401e-14, "6.82613366240855811298682466187049063105090237181643243484884807542082161846280596239346667442720742223177409347157282675e-205", false, false); + yield return new("Exp2", -370.71699147959754, -1.4210854715202004e-14, "2.52968043778232094121601606971876522062602966535407014705556941483082015232179403400791844339483992427445921351569246286e-112", false, false); + yield return new("Exp2", 325.3049427158436, 1.4210854715202004e-14, "8.44394693210093550736509533332028359812000552799434020131118049620518568591929165240041457560841465860244163260586337592e+97", false, false); + yield return new("Exp2", 956.9782395104039, 2.842170943040401e-14, "1.19992828387600858106093631065960712137870691174075367927710136390487281510229839467771780838917585057717098231260902141e+288", false, false); + yield return new("Exp2", 364.3581585496251, -1.4210854715202004e-14, "4.81653686139180239709139190167158242071136853423964187582462482485565537885773294376979018541476826443148951897942014775e+109", false, false); + // Exp2: integer-power + yield return new("Exp2", -1076.0, 0.0, "1.23516411460311636044142198217055343091264950653581191106396420625168876817552187966324959090408998094949141173861429433e-324", false, true); + yield return new("Exp2", -1075.0, 0.0, "2.47032822920623272088284396434110686182529901307162382212792841250337753635104375932649918180817996189898282347722858865e-324", false, true); + yield return new("Exp2", -1074.0, 0.0, "4.94065645841246544176568792868221372365059802614324764425585682500675507270208751865299836361635992379796564695445717731e-324", false, false); + yield return new("Exp2", -1022.0, 0.0, "2.22507385850720138309023271733240406421921598046233183055332741688720443481391819585428315901251102056406733973103581101e-308", false, false); + yield return new("Exp2", -100.0, 0.0, "7.888609052210118054117285652827862296732064351090230047702789306640625e-31", false, false); + yield return new("Exp2", 10.0, 0.0, "1.024e+3", false, false); + yield return new("Exp2", 53.0, 0.0, "9.007199254740992e+15", false, false); + yield return new("Exp2", 100.0, 0.0, "1.267650600228229401496703205376e+30", false, false); + yield return new("Exp2", 1023.0, 0.0, "8.98846567431157953864652595394512366808988489471153286367150405788663379027504815663542386612037680105600569399356966788e+307", false, false); + yield return new("Exp2", 1024.0, 0.0, "1.79769313486231590772930519078902473361797697894230657273430081157732675805500963132708477322407536021120113879871393358e+308", true, false); + // Exp2: overflow-adjacent-low + yield return new("Exp2", 1024.0, -8.008566259537295e-17, "1.79769313486231580793728971405302296962136452782486920563609826833849467898852832283391964137375521264887063753644165960e+308", false, false); + yield return new("Exp2", 1024.0, -8.008566259537294e-17, "1.79769313486231580793728971405303832856973767181831156142887561931500193969224300858151352073314340505642260157691456128e+308", true, false); + yield return new("Exp2", 1024.0, -8.008566259537293e-17, "1.79769313486231580793728971405305368751811081581175391722165297042273142656692749614371314702169582876170167447259768366e+308", true, false); + // Exp2: underflow-adjacent-low + yield return new("Exp2", -1075.0, -5e-324, "2.47032822920623272088284396434110686182529901307162382212792841250337753635104375932649918180817996189898282347722858865e-324", false, true); + yield return new("Exp2", -1075.0, 5e-324, "2.47032822920623272088284396434110686182529901307162382212792841250337753635104375932649918180817996189898282347722858865e-324", false, false); + // Exp2: min-normal-adjacent-low + yield return new("Exp2", -1022.0, -5e-324, "2.22507385850720138309023271733240406421921598046233183055332741688720443481391819585428315901251102056406733973103581101e-308", false, false); + yield return new("Exp2", -1022.0, 5e-324, "2.22507385850720138309023271733240406421921598046233183055332741688720443481391819585428315901251102056406733973103581101e-308", false, false); + // Exp10: ordinary + yield return new("Exp10", -10.25, 0.0, "5.62341325190349080394951039776481231468251043098691664081689423735883568643062848905857984526220305928676107320100325218e-11", false, false); + yield return new("Exp10", -2.0, 0.0, "1e-2", false, false); + yield return new("Exp10", -1.0, 0.0, "1e-1", false, false); + yield return new("Exp10", -0.125, 0.0, "7.49894209332455827302184275615136438441867918164971014620419005429827525167160627980673695983144556246592084007724058545e-1", false, false); + yield return new("Exp10", 0.125, 0.0, "1.33352143216332402567593171529533109241566796476437099332954998716275894318019581864901349800473255887744566135767837809e+0", false, false); + yield return new("Exp10", 0.75, 0.0, "5.62341325190349080394951039776481231468251043098691664081689423735883568643062848905857984526220305928676107320100325218e+0", false, false); + yield return new("Exp10", 1.0, 0.0, "1.0e+1", false, false); + yield return new("Exp10", 3.25, 0.0, "1.77827941003892280122542119519268484473579052640225535801183072277630188153949380490030039927870215508827904815953580779e+3", false, false); + // Exp10: dense-and-sparse-low + yield return new("Exp10", -1.0, 5.551115123125783e-17, "1.00000000000000012781914932003234537245680382409083392670440261696902368520426109154220545198530672142968125154509195453e-1", false, false); + yield return new("Exp10", -1.0, -5.551115123125783e-17, "9.99999999999999872180850679967670965278129072632256181650151054098060763085821567736360614288865296549471193611253084211e-2", false, false); + yield return new("Exp10", -1.0, 5e-324, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-1", false, false); + yield return new("Exp10", -1.0, -5e-324, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-1", false, false); + yield return new("Exp10", 1.0, 5.551115123125783e-17, "1.00000000000000012781914932003234537245680382409083392670440261696902368520426109154220545198530672142968125154509195453e+1", false, false); + yield return new("Exp10", 1.0, -5.551115123125783e-17, "9.99999999999999872180850679967670965278129072632256181650151054098060763085821567736360614288865296549471193611253084211e+0", false, false); + yield return new("Exp10", 1.0, 5e-324, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+1", false, false); + yield return new("Exp10", 1.0, -5e-324, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+1", false, false); + // Exp10: tiny + yield return new("Exp10", -1e-20, 0.0, "9.99999999999999999976974149070059544422977775848960612440147585677561893017435981776880676797244875417274786359562174097e-1", false, false); + yield return new("Exp10", 1e-20, 0.0, "1.00000000000000000002302585092994045557755241396208722730275022767272793307903996143344756701626889076812422686037317339e+0", false, false); + yield return new("Exp10", -1e-100, 0.0, "9.99999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999769741490700595426995e-1", false, false); + yield return new("Exp10", 1e-100, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000023025850929940457301e+0", false, false); + yield return new("Exp10", -1e-300, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp10", 1e-300, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp10", -2.2250738585072014e-308, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp10", 2.2250738585072014e-308, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp10", -1e-323, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp10", 1e-323, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp10", -5e-324, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp10", 5e-324, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + // Exp10: absolute-log-half-switch + yield return new("Exp10", -0.2171472409516259, -5.491598251083826e-18, "6.06530659712633423603799534991179637117927421578992542354041909882833935611531036427818598568463751669923025111585399002e-1", false, false); + yield return new("Exp10", -0.2171472409516259, -5.491598251083825e-18, "6.06530659712633423603799534991180713010537993345061928799291965018127819688823024941030203305141921320141807673626493593e-1", false, false); + yield return new("Exp10", -0.2171472409516259, -5.4915982510838246e-18, "6.06530659712633423603799534991181788903148565111131315244542020155330172680170174417944401325366675980091613028373016444e-1", false, false); + yield return new("Exp10", 0.2171472409516259, 5.4915982510838246e-18, "1.64872127070012814684865078781415994149378081230563234818972414648475614387142248427209749759453651737942752570661371215e+0", false, false); + yield return new("Exp10", 0.2171472409516259, 5.491598251083825e-18, "1.64872127070012814684865078781416286607311350290140527066142468471602253792023469276276331171381066284701617511531416884e+0", false, false); + yield return new("Exp10", 0.2171472409516259, 5.491598251083826e-18, "1.64872127070012814684865078781416579065244619349717819313312522295247668833830201265250709635897011641393592320394302761e+0", false, false); + // Exp10: random-moderate + yield return new("Exp10", -2.9677357165710085, -1.1102230246251565e-16, "1.07712047983496528579184014784924244267541657758050199455170643282791807045975507835023158082182189438408871295653951051e-3", false, false); + yield return new("Exp10", 0.5995899657426765, -2.7755575615628914e-17, "3.97731479519890352879887600772233924298609686175885392875139750363632580171965156242522363738652933349600506618119299140e+0", false, false); + yield return new("Exp10", 0.8836851964549527, -2.7755575615628914e-17, "7.65041856144971858829962647061748349456294701097773926744045197071139270643977619590118903694461196808465346045233138263e+0", false, false); + yield return new("Exp10", -0.89395009436352, -2.7755575615628914e-17, "1.27658549534319782190947058351173152520381158276033885330906540617729206451590520184715981551175953612070355802078326647e-1", false, false); + yield return new("Exp10", 2.4580147142578235, -1.1102230246251565e-16, "2.87087784810333896489677828661765156448804286447564220961070487382944745109858177125943722517534602403829765495652338124e+2", false, false); + yield return new("Exp10", 1.3118516264245281, 5.551115123125783e-17, "2.05046153376223160193906750616367118569485139159211538092010333804564633619337531743249262303443651659625057726925490890e+1", false, false); + // Exp10: random-range + yield return new("Exp10", 120.58969150830279, 3.552713678800501e-15, "3.88768893311729066468034240903516540867755441511370722931132770507875572390376800025883066808459259868657675470129694955e+120", false, false); + yield return new("Exp10", 128.81118142762747, -7.105427357601002e-15, "6.47413017703446171983916778641885806651744797264667468232608054450205762476305583291383873034783101822647556402661401770e+128", false, false); + yield return new("Exp10", 42.94992135127048, 1.7763568394002505e-15, "8.91089551309232685378712887186962248614334767700386757268709525041370952477415725129562497023716824045224543686275322633e+42", false, false); + yield return new("Exp10", 116.81408221381518, -3.552713678800501e-15, "6.51751761797928744146122574992738336073782813915351110967844224366907044327161884077729178445278842875234452038465777427e+116", false, false); + yield return new("Exp10", 228.96449953361966, -7.105427357601002e-15, "9.21508899144045996730372464171725339874292277767159122312745668742022577561608598769932754185213295353703683516203829072e+228", false, false); + yield return new("Exp10", 7.524609271856807, 2.220446049250313e-16, "3.34664211369735369097716632117090216287459757705675042241703792408934408612576386525573391905086740939515500212837287795e+7", false, false); + // Exp10: integer-power + yield return new("Exp10", -324.0, 0.0, "1e-324", false, true); + yield return new("Exp10", -323.0, 0.0, "1e-323", false, false); + yield return new("Exp10", -308.0, 0.0, "1e-308", false, false); + yield return new("Exp10", 22.0, 0.0, "1.0000000000000000000000e+22", false, false); + yield return new("Exp10", 308.0, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+308", false, false); + yield return new("Exp10", 309.0, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+309", true, false); + // Exp10: overflow-adjacent-low + yield return new("Exp10", 308.25471555991675, -2.895125789515847e-15, "1.79769313486231580793728971405135339792924520420229515147549215774814264965213204794170642362509327672779802464860968196e+308", false, false); + yield return new("Exp10", 308.25471555991675, -2.8951257895158466e-15, "1.79769313486231580793728971405298608023672269474350236922097481966929437654805569985823689453925527448616186188699371312e+308", false, false); + yield return new("Exp10", 308.25471555991675, -2.8951257895158462e-15, "1.79769313486231580793728971405461876254420018528470958696645896440832277858828718822634729443205898882595653371847623943e+308", true, false); + // Exp10: underflow-adjacent-low + yield return new("Exp10", -323.60724533877976, -2.2576632980240802e-14, "2.47032822920623272088284396433158418397152352476997073314847977407861399577439129056325691274224699350470920929775150941e-324", false, true); + yield return new("Exp10", -323.60724533877976, -2.25766329802408e-14, "2.47032822920623272088284396434953279206384344178901504121903931494189587104807927555924014551585189156808716237545881368e-324", false, false); + yield return new("Exp10", -323.60724533877976, -2.2576632980240796e-14, "2.47032822920623272088284396436748140015616335880805934928972926460213150426298975260692893446611258989980884791171466090e-324", false, false); + // Exp10: min-normal-adjacent-low + yield return new("Exp10", -307.6526555685888, 2.754387844133928e-15, "2.22507385850720138309023271732978958155884829152404782449780586939390649146037065767628694530848808276307239212102756320e-308", false, false); + yield return new("Exp10", -307.6526555685888, 2.7543878441339283e-15, "2.22507385850720138309023271733181041517675806770793310094027818310578520394298980281484276867693115447731727042755624095e-308", false, false); + yield return new("Exp10", -307.6526555685888, 2.7543878441339287e-15, "2.22507385850720138309023271733383124879466784389181837738275233215831818760979011696161815647039528165171309747797451182e-308", false, false); + // ExpM1: ordinary + yield return new("ExpM1", -10.25, 0.0, "-9.99964642499149590017595412360236722574820646856009376883960968787652399508581301416230500941667060422476262375531037864e-1", false, false); + yield return new("ExpM1", -2.0, 0.0, "-8.64664716763387308106000505027515596592368454090424118531841127345926625898512310062901877509342951244922712810366447788e-1", false, false); + yield return new("ExpM1", -1.0, 0.0, "-6.32120558828557678404476229838539132554188868968232165492163198302538504255100196642852725654080356253372674723156004792e-1", false, false); + yield return new("ExpM1", -0.125, 0.0, "-1.17503097415404597135107856770949263777995175009349258229690680791890704788691269900644852670056561354676484984318222954e-1", false, false); + yield return new("ExpM1", 0.125, 0.0, "1.33148453066826316829007227811793872565503131745181625912820036078823577880048386513939990794941728573231527015647307566e-1", false, false); + yield return new("ExpM1", 0.75, 0.0, "1.11700001661267466854536981983709561013449158470240342177913303081098453336401282000279156026661579821888590471901551426e+0", false, false); + yield return new("ExpM1", 1.0, 0.0, "1.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193200305992e+0", false, false); + yield return new("ExpM1", 3.25, 0.0, "2.47903399171930620890801076693772218766552688489539736340424977271388720447207274888484147322397605262439439122840412454e+1", false, false); + // ExpM1: dense-and-sparse-low + yield return new("ExpM1", -1.0, 5.551115123125783e-17, "-6.32120558828557657983064936099986792792127784222576303251939251170289698165054910720245453994927423599213941197314427262e-1", false, false); + yield return new("ExpM1", -1.0, -5.551115123125783e-17, "-6.32120558828557698825887523577090338700199271276585542169177230003595532046796634595998012127496903342004079921719212416e-1", false, false); + yield return new("ExpM1", -1.0, 5e-324, "-6.32120558828557678404476229838539132554188868968232165492163198302538504255100196642852725654080356253372674723156004792e-1", false, false); + yield return new("ExpM1", -1.0, -5e-324, "-6.32120558828557678404476229838539132554188868968232165492163198302538504255100196642852725654080356253372674723156004792e-1", false, false); + yield return new("ExpM1", 1.0, 5.551115123125783e-17, "1.71828182845904538625524114012277294245207378041324564716578593432480201210931242665209562398991569452359995574369903890e+0", false, false); + yield return new("ExpM1", 1.0, -5.551115123125783e-17, "1.71828182845904508446533380258256042941501354772712009605959080847470919809942371990293379281005694591784351622350569198e+0", false, false); + yield return new("ExpM1", 1.0, 5e-324, "1.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193200305992e+0", false, false); + yield return new("ExpM1", 1.0, -5e-324, "1.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193200305992e+0", false, false); + // ExpM1: tiny + yield return new("ExpM1", -1e-20, 0.0, "-9.99999999999999945148271454209571652277987654911963382213526866733506442994862496388124890915444431248853256190340604605e-21", false, false); + yield return new("ExpM1", 1e-20, 0.0, "9.99999999999999945158271454209571651181053083996154845329753345878813976563872245650392070698072643274472760230268014801e-21", false, false); + yield return new("ExpM1", -1e-100, 0.0, "-1.00000000000000001999189980260288361964776078853415942018260300593659569925554346761767628861329298953274607481091184880e-100", false, false); + yield return new("ExpM1", 1e-100, 0.0, "1.00000000000000001999189980260288361964776078853415942018260300593659569925554346761767628861329298963274607481091185280e-100", false, false); + yield return new("ExpM1", -1e-300, 0.0, "-1.00000000000000002505909183520875968569614680770370524992534231990046604318405148467630281218195010089496230627027825415e-300", false, false); + yield return new("ExpM1", 1e-300, 0.0, "1.00000000000000002505909183520875968569614680770370524992534231990046604318405148467630281218195010089496230627027825415e-300", false, false); + yield return new("ExpM1", -2.2250738585072014e-308, 0.0, "-2.22507385850720138309023271733240406421921598046233183055332741688720443481391819585428315901251102056406733973103581101e-308", false, false); + yield return new("ExpM1", 2.2250738585072014e-308, 0.0, "2.22507385850720138309023271733240406421921598046233183055332741688720443481391819585428315901251102056406733973103581101e-308", false, false); + yield return new("ExpM1", -1e-323, 0.0, "-9.88131291682493088353137585736442744730119605228649528851171365001351014540417503730599672723271984759593129390891435462e-324", false, false); + yield return new("ExpM1", 1e-323, 0.0, "9.88131291682493088353137585736442744730119605228649528851171365001351014540417503730599672723271984759593129390891435462e-324", false, false); + yield return new("ExpM1", -5e-324, 0.0, "-4.94065645841246544176568792868221372365059802614324764425585682500675507270208751865299836361635992379796564695445717731e-324", false, false); + yield return new("ExpM1", 5e-324, 0.0, "4.94065645841246544176568792868221372365059802614324764425585682500675507270208751865299836361635992379796564695445717731e-324", false, false); + // ExpM1: absolute-log-half-switch + yield return new("ExpM1", -0.5, -5e-324, "-3.93469340287366576396200465008819546558081864512813044317107841264943480586251576001352388492010543973576210205960474823e-1", false, false); + yield return new("ExpM1", -0.5, 0.0, "-3.93469340287366576396200465008819546558081864512813044317107841264943480586251576001352388492010543973576210205960474823e-1", false, false); + yield return new("ExpM1", -0.5, 5e-324, "-3.93469340287366576396200465008819546558081864512813044317107841264943480586251576001352388492010543973576210205960474823e-1", false, false); + yield return new("ExpM1", 0.5, -5e-324, "6.48721270700128146848650787814163571653776100710148011575079311640661021194215608632776520056366643002866637756307797005e-1", false, false); + yield return new("ExpM1", 0.5, 0.0, "6.48721270700128146848650787814163571653776100710148011575079311640661021194215608632776520056366643002866637756307797005e-1", false, false); + yield return new("ExpM1", 0.5, 5e-324, "6.48721270700128146848650787814163571653776100710148011575079311640661021194215608632776520056366643002866637756307797005e-1", false, false); + // ExpM1: input-high-series-switch + yield return new("ExpM1", -0.5000000000000001, 0.0, "-3.93469340287366643734630820413961752893910543944792323093955342455151126428014543140163447963889188936567569177072330420e-1", false, false); + yield return new("ExpM1", -0.49999999999999994, 0.0, "-3.93469340287366542726985287306245639864824173983524003239030957439240694784182493508171378245214271234180273050238838247e-1", false, false); + yield return new("ExpM1", -0.5, -2.7755575615628914e-17, "-3.93469340287366593230808053860105799023374872074113261068818573492447875650538514346733263349345635337465711156202652272e-1", false, false); + yield return new("ExpM1", -0.5, 2.7755575615628914e-17, "-3.93469340287366559561592876157532826838564965149287614677159619553512682806916609709611128939222384005591032211301092594e-1", false, false); + // ExpM1: quadratic-series-switch + yield return new("ExpM1", -5.551115123125783e-17, -5e-324, "-5.55111512312578254804376278947523616897078869154779750199182577185527323675307989703604544287548416109752829359640238674e-17", false, false); + yield return new("ExpM1", -5.551115123125783e-17, 0.0, "-5.55111512312578254804376278947523616897078869154779750199182577185527323675307989703604544287548416109752829359640238674e-17", false, false); + yield return new("ExpM1", -5.551115123125783e-17, 5e-324, "-5.55111512312578254804376278947523616897078869154779750199182577185527323675307989703604544287548416109752829359640238674e-17", false, false); + // ExpM1: input-high-series-switch + yield return new("ExpM1", 0.49999999999999994, 0.0, "6.48721270700128055326234991787778407132549255798122477493863918355113649772615017427283934825041206826640936896304389788e-1", false, false); + yield return new("ExpM1", 0.5000000000000001, 0.0, "6.48721270700128329893482379866949142240222700377506852855419426513165259729393218477113635478077785280838920026492973826e-1", false, false); + yield return new("ExpM1", 0.5, -2.7755575615628914e-17, "6.48721270700128101087442889800970354328829640344026798281419430748374005506759742356897159987521249118694442987495686887e-1", false, false); + yield return new("ExpM1", 0.5, 2.7755575615628914e-17, "6.48721270700128192609858685827358059107388636896521370527076406346895321998228741244042216578759957646849801680426913917e-1", false, false); + // ExpM1: quadratic-series-switch + yield return new("ExpM1", 5.551115123125783e-17, -5e-324, "5.55111512312578285619255389143297265792725950513624759781930415130694675216220020400320592376210494940133730810649727575e-17", false, false); + yield return new("ExpM1", 5.551115123125783e-17, 0.0, "5.55111512312578285619255389143297265792725950513624759781930415130694675216220020400320592376210494940133730810649727575e-17", false, false); + yield return new("ExpM1", 5.551115123125783e-17, 5e-324, "5.55111512312578285619255389143297265792725950513624759781930415130694675216220020400320592376210494940133730810649727575e-17", false, false); + // ExpM1: random-moderate + yield return new("ExpM1", 1.8665732301127385, 5.551115123125783e-17, "5.46610055231604164741021410411512794541063980685522084739073087713841387734128034155842973464987332959651169399776074928e+0", false, false); + yield return new("ExpM1", -3.836275153651057, -1.1102230246251565e-16, "-9.78426189000834536644301856044613215358744428231199467965735926574776462346821532341070704496171040555177164906750092116e-1", false, false); + yield return new("ExpM1", -7.926421267629822, 2.220446049250313e-16, "-9.99638923701085515458633712931479267616303538413149722750742992858886680422776234425246545391734902973472717556638544616e-1", false, false); + yield return new("ExpM1", 0.809076467549632, -2.7755575615628914e-17, "1.24583292904435095047625826229914076404174762338755219183658975528974302916271313844842172254318969775653936199802268813e+0", false, false); + yield return new("ExpM1", 2.395252849322306, 1.1102230246251565e-16, "9.97097171116982451435344657554439752704646874789848866365033063395867538506970643915392963325796240244200998541676010281e+0", false, false); + yield return new("ExpM1", 6.360537727726566, 2.220446049250313e-16, "5.77557379107619288839611291974544367870013144362286685855659151210210389821741943112835250334310901695746967316162020030e+2", false, false); + // ExpM1: random-range + yield return new("ExpM1", -717.471294689501, -2.842170943040401e-14, "-1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("ExpM1", -527.7157218253692, -2.842170943040401e-14, "-1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("ExpM1", -322.9587865258026, -1.4210854715202004e-14, "-1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("ExpM1", -24.19275074210509, -8.881784197001252e-16, "-9.99999999968866936988232338142807791500416465558403067189833971761249730793225364843678313704660239258534490554851101589e-1", false, false); + yield return new("ExpM1", -26.228830895869123, -8.881784197001252e-16, "-9.99999999995935908519211296503644770852506238245358308752490322883971257357333755851647526243130367097065447570813612841e-1", false, false); + yield return new("ExpM1", 407.0814289754503, 1.4210854715202004e-14, "6.21181179920291117447313962056204710364585787558821473816816695867214972444088839539628676739709899564766755487635174426e+176", false, false); + // ExpM1: negative-saturation + yield return new("ExpM1", -40.0, 0.0, "-9.99999999999999995751645744708411004670765217141341982120434445833553711949181081073966936073085345895610771405272219090e-1", false, false); + yield return new("ExpM1", -100.0, 0.0, "-9.99999999999999999999999999999999999999999962799240239791640370403041961368816626411077076232180328793861233367095241042e-1", false, false); + yield return new("ExpM1", -745.0, 0.0, "-1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("ExpM1", -746.0, 0.0, "-1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + // ExpM1: overflow-adjacent-low + yield return new("ExpM1", 709.782712893384, 2.3691528222554846e-14, "1.79769313486231580793728971404631370472204944707885473394255307689401631566863219557545655233948936948926513968471917798e+308", false, false); + yield return new("ExpM1", 709.782712893384, 2.369152822255485e-14, "1.79769313486231580793728971405198622405675778669468436196455240891719566132544008818951355194454679628152242602530089901e+308", false, false); + yield return new("ExpM1", 709.782712893384, 2.3691528222554853e-14, "1.79769313486231580793728971405765874339146612631051398998656964025532405383200713096445663417114297060553540946791444483e+308", true, false); + // Exp2M1: ordinary + yield return new("Exp2M1", -10.25, 0.0, "-9.99178812094478794391571166527116001079062466540667202626158958900415112543990329843931324214649128382643748463141094643e-1", false, false); + yield return new("Exp2M1", -2.0, 0.0, "-7.5e-1", false, false); + yield return new("Exp2M1", -1.0, 0.0, "-5e-1", false, false); + yield return new("Exp2M1", -0.125, 0.0, "-8.29959567953287682564584052058555719613448335631602502083379306467611688776526371585345106193396863276589403310505403211e-2", false, false); + yield return new("Exp2M1", 0.125, 0.0, "9.05077326652576592070106557607079789927027185400671217856676476833005308488418403382111404942031198914516192629180900103e-2", false, false); + yield return new("Exp2M1", 0.75, 0.0, "6.81792830507429086062250952466429790080068524713569021626452171949849509907804479628648008398585072345603147487038170160e-1", false, false); + yield return new("Exp2M1", 1.0, 0.0, "1e+0", false, false); + yield return new("Exp2M1", 3.25, 0.0, "8.51365692002176853373999976448380732234377673971053930415201779775573334581533727896625075630510138972829915819815370549e+0", false, false); + // Exp2M1: dense-and-sparse-low + yield return new("Exp2M1", -1.0, 5.551115123125783e-17, "-4.99999999999999980761301017208449927395839145715461584550039659526711955921405213806161306947046712263190661419361655482e-1", false, false); + yield return new("Exp2M1", -1.0, -5.551115123125783e-17, "-5.00000000000000019238698982791549332349083753355319945239272896296038521360689797886742152734429064724448187428141856202e-1", false, false); + yield return new("Exp2M1", -1.0, 5e-324, "-5.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-1", false, false); + yield return new("Exp2M1", -1.0, -5e-324, "-5.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-1", false, false); + yield return new("Exp2M1", 1.0, 5.551115123125783e-17, "1.00000000000000007695479593116620029041664341713815366179984136189315217631437914477535477221181315094723735432255337807e+0", false, false); + yield return new("Exp2M1", 1.0, -5.551115123125783e-17, "9.99999999999999923045204068833802670603664986578720219042908414815845914557240808453031389062283741102207250287432575192e-1", false, false); + yield return new("Exp2M1", 1.0, 5e-324, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp2M1", 1.0, -5e-324, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + // Exp2M1: tiny + yield return new("Exp2M1", -1e-20, 0.0, "-6.93147180559945271397974601937280894524371589944357697762013153681278242681734018363320249079658219031934314131211295962e-21", false, false); + yield return new("Exp2M1", 1e-20, 0.0, "6.93147180559945271402779132076462908244017094540055080820895837340445945043791570869517416287708473950733796857874762083e-21", false, false); + yield return new("Exp2M1", -1e-100, 0.0, "-6.93147180559945323274561103669292529709311557561290569992180233374923094914456585984237940738350279587910585083533591459e-101", false, false); + yield return new("Exp2M1", 1e-100, 0.0, "6.93147180559945323274561103669292529709311557561290569992180233374923094914456585984237940738350279635955886475353735846e-101", false, false); + yield return new("Exp2M1", -1e-300, 0.0, "-6.93147180559945326786870974425873986894274700919374605593673636620279660269393102485807676278657875062914284858935724438e-301", false, false); + yield return new("Exp2M1", 1e-300, 0.0, "6.93147180559945326786870974425873986894274700919374605593673636620279660269393102485807676278657875062914284858935724438e-301", false, false); + yield return new("Exp2M1", -2.2250738585072014e-308, 0.0, "-1.54230367156190531855910393712847096456890294742459890370396553388262870291588706079233655504875900365429627823154008399e-308", false, false); + yield return new("Exp2M1", 2.2250738585072014e-308, 0.0, "1.54230367156190531855910393712847096456890294742459890370396553388262870291588706079233655504875900365429627823154008399e-308", false, false); + yield return new("Exp2M1", -1e-323, 0.0, "-6.84920418852777021427470999061380292158364356780264127745978767817293643334458446499873465218026337654986814054831649295e-324", false, false); + yield return new("Exp2M1", 1e-323, 0.0, "6.84920418852777021427470999061380292158364356780264127745978767817293643334458446499873465218026337654986814054831649295e-324", false, false); + yield return new("Exp2M1", -5e-324, 0.0, "-3.42460209426388510713735499530690146079182178390132063872989383908646821667229223249936732609013168827493407027415824647e-324", false, false); + yield return new("Exp2M1", 5e-324, 0.0, "3.42460209426388510713735499530690146079182178390132063872989383908646821667229223249936732609013168827493407027415824647e-324", false, false); + // Exp2M1: absolute-log-half-switch + yield return new("Exp2M1", -0.7213475204444817, -1.0177636870465518e-17, "-3.93469340287366576396200465008820417438099792541202652192428640510984866996467126609246899789618731498509990689191240306e-1", false, false); + yield return new("Exp2M1", -0.7213475204444817, -1.0177636870465517e-17, "-3.93469340287366576396200465008819769686204001884905656907964381721986981887448485558283931634670370961885669141173839528e-1", false, false); + yield return new("Exp2M1", -0.7213475204444817, -1.0177636870465515e-17, "-3.93469340287366576396200465008819121934308211228608661623500122932297322155364322708244878742394153333622619986830147957e-1", false, false); + yield return new("Exp2M1", 0.7213475204444817, 1.0177636870465515e-17, "6.48721270700128146848650787814162417406688247281002602415805950728600029067384069417210842542030614979026168642397668251e-1", false, false); + yield return new("Exp2M1", 0.7213475204444817, 1.0177636870465517e-17, "6.48721270700128146848650787814164178178895924919128376762790092850572717154172663805292837189744464569525528256254871274e-1", false, false); + yield return new("Exp2M1", 0.7213475204444817, 1.0177636870465518e-17, "6.48721270700128146848650787814165938951103602557254151109774234974425843628229371597802606772669504937810418198966839579e-1", false, false); + // Exp2M1: input-high-series-switch + yield return new("Exp2M1", -0.5000000000000001, 0.0, "-2.92893218813452530014413685649703955167858165405676925733374390585360532699651656717067373024344502855553936702073773310e-1", false, false); + yield return new("Exp2M1", -0.5, 0.0, "-2.92893218813452475599155637895150960715164062311525963411660131004633760768946480574806232836179213632493076884543851488e-1", false, false); + yield return new("Exp2M1", -0.49999999999999994, 0.0, "-2.92893218813452448391526614017872893170662533251599103883499092807605318109978129696795891336476533197262084559728907784e-1", false, false); + yield return new("Exp2M1", -0.5, -2.7755575615628914e-17, "-2.92893218813452489202970149833789601907876207463268995864344371593967907826653024658969061732863275170576563234829479896e-1", false, false); + yield return new("Exp2M1", -0.5, 2.7755575615628914e-17, "-2.92893218813452461995341125956512057802759504240971011133505105055334106579225243377019120609626314024272048921794865070e-1", false, false); + // Exp2M1: quadratic-series-switch + yield return new("Exp2M1", -8.008566259537294e-17, -1.129944678986474e-33, "-5.55111512312578254804376278947525210997941054353275279722773648933255116925354923236018969777309126680622446520049198092e-17", false, false); + yield return new("Exp2M1", -8.008566259537294e-17, -1.1299446789864738e-33, "-5.55111512312578254804376278947524025321574270017331675503562752039604580100758672073788499063969008912857194783169179131e-17", false, false); + yield return new("Exp2M1", -8.008566259537294e-17, -1.1299446789864737e-33, "-5.55111512312578254804376278947522839645207485681388071284351855145954043276162420770975183675538569412491789275246025718e-17", false, false); + // Exp2M1: input-high-series-switch + yield return new("Exp2M1", 0.49999999999999994, 0.0, "4.14213562373094994386430676455144037238208120607589712724123944476593389531373838698287076714161408538277123570716809775e-1", false, false); + yield return new("Exp2M1", 0.5, 0.0, "4.14213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641572735013846230912297025e-1", false, false); + yield return new("Exp2M1", 0.5000000000000001, 0.0, "4.14213562373095157632204819718812442505217294967231432235173388674957037169706704669406056199264812251701573141419622976e-1", false, false); + yield return new("Exp2M1", 0.5, -2.7755575615628914e-17, "4.14213562373095021594059700332420796184247585073462008271311256812064184346693950682061876534273449658846873530341040209e-1", false, false); + yield return new("Exp2M1", 0.5, 2.7755575615628914e-17, "4.14213562373095076009317748086975884394480991518057977732989789889331786841549513245961758780747371951455902156410269860e-1", false, false); + // Exp2M1: quadratic-series-switch + yield return new("Exp2M1", 8.008566259537294e-17, 1.1299446789864737e-33, "5.55111512312578285619255389143296488540854567040146788574745400787982691625510656198196538826488518799624108409547153570e-17", false, false); + yield return new("Exp2M1", 8.008566259537294e-17, 1.1299446789864738e-33, "5.55111512312578285619255389143297674217221351376222029314172084883206635846687879210449178738034357857131792648832026482e-17", false, false); + yield return new("Exp2M1", 8.008566259537294e-17, 1.129944678986474e-33, "5.55111512312578285619255389143298859893588135712297270053598768978430580067865102363284663324670534255070733217896828964e-17", false, false); + // Exp2M1: random-moderate + yield return new("Exp2M1", 8.78452800194328, -4.440892098500626e-16, "4.39967341476995196084861795619177036986833446346865068997727335172090686266840323292806805468322290875886544488139632805e+2", false, false); + yield return new("Exp2M1", 2.7202982684631807, 1.1102230246251565e-16, "5.59009045860470968459866968889737028605350979742874890837108301124392819835145507472654892784605100090392179416858372658e+0", false, false); + yield return new("Exp2M1", -3.4144589655515833, 1.1102230246251565e-16, "-9.06212397865302249501452386547149636442561439125894793225413740615446561503933959009917468149743846696349004220925890981e-1", false, false); + yield return new("Exp2M1", 3.67643978981615, 1.1102230246251565e-16, "1.17855275975436677510222614327660438420381423323350572823327975657260641462810932407195792018571153800381071193261125222e+1", false, false); + yield return new("Exp2M1", 7.559099749783131, -2.220446049250313e-16, "1.87588742600978100098560391503204728439903232936850014376858946537478147694222110588671754338725081209043377353414161445e+2", false, false); + yield return new("Exp2M1", 1.1272538511782404, 5.551115123125783e-17, "1.18442542093619461179505632322739881400577647657120007022725830276626252983657181904591357904539459020575218758691298325e+0", false, false); + // Exp2M1: random-range + yield return new("Exp2M1", -816.2214079897565, -2.842170943040401e-14, "-1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp2M1", -678.2241907878085, -2.842170943040401e-14, "-1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp2M1", -370.71699147959754, -1.4210854715202004e-14, "-9.99999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999747031956e-1", false, false); + yield return new("Exp2M1", 325.3049427158436, 1.4210854715202004e-14, "8.44394693210093550736509533332028359812000552799434020131118049620518568591929165240041457560841455860244163260586337592e+97", false, false); + yield return new("Exp2M1", 956.9782395104039, 2.842170943040401e-14, "1.19992828387600858106093631065960712137870691174075367927710136390487281510229839467771780838917585057717098231260902141e+288", false, false); + yield return new("Exp2M1", 364.3581585496251, -1.4210854715202004e-14, "4.81653686139180239709139190167158242071136853423964187582462482485565537885773294376979018541476826443148951887942014775e+109", false, false); + // Exp2M1: integer-power + yield return new("Exp2M1", -1076.0, 0.0, "-1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp2M1", -1075.0, 0.0, "-1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp2M1", -1074.0, 0.0, "-1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp2M1", -1022.0, 0.0, "-1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp2M1", -100.0, 0.0, "-9.999999999999999999999999999992111390947789881945882714347172137703267935648909769952297210693359375e-1", false, false); + yield return new("Exp2M1", 10.0, 0.0, "1.023e+3", false, false); + yield return new("Exp2M1", 53.0, 0.0, "9.007199254740991e+15", false, false); + yield return new("Exp2M1", 100.0, 0.0, "1.267650600228229401496703205375e+30", false, false); + yield return new("Exp2M1", 1023.0, 0.0, "8.98846567431157953864652595394512366808988489471153286367150405788663379027504815663542386612037680105600569399356966788e+307", false, false); + yield return new("Exp2M1", 1024.0, 0.0, "1.79769313486231590772930519078902473361797697894230657273430081157732675805500963132708477322407536021120113879871393358e+308", true, false); + // Exp2M1: negative-saturation + yield return new("Exp2M1", -54.0, 0.0, "-9.99999999999999944488848768742172978818416595458984375e-1", false, false); + // Exp2M1: negative-saturation-sparse-low + yield return new("Exp2M1", -1075.0, -5e-324, "-1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp2M1", -1075.0, 5e-324, "-1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + // Exp2M1: overflow-adjacent-low + yield return new("Exp2M1", 1024.0, -8.008566259537295e-17, "1.79769313486231580793728971405302296962136452782486920563609826833849467898852832283391964137375521264887063753644165960e+308", false, false); + yield return new("Exp2M1", 1024.0, -8.008566259537294e-17, "1.79769313486231580793728971405303832856973767181831156142887561931500193969224300858151352073314340505642260157691456128e+308", true, false); + yield return new("Exp2M1", 1024.0, -8.008566259537293e-17, "1.79769313486231580793728971405305368751811081581175391722165297042273142656692749614371314702169582876170167447259768366e+308", true, false); + // Exp10M1: ordinary + yield return new("Exp10M1", -10.25, 0.0, "-9.99999999943765867480965091960504896022351876853174895690130833591831057626411643135693715109414201547377969407132389268e-1", false, false); + yield return new("Exp10M1", -2.0, 0.0, "-9.9e-1", false, false); + yield return new("Exp10M1", -1.0, 0.0, "-9e-1", false, false); + yield return new("Exp10M1", -0.125, 0.0, "-2.50105790667544172697815724384863561558132081835028985379580994570172474832839372019326304016855443753407915992275941455e-1", false, false); + yield return new("Exp10M1", 0.125, 0.0, "3.33521432163324025675931715295331092415667964764370993329549987162758943180195818649013498004732558877445661357678378086e-1", false, false); + yield return new("Exp10M1", 0.75, 0.0, "4.62341325190349080394951039776481231468251043098691664081689423735883568643062848905857984526220305928676107320100325218e+0", false, false); + yield return new("Exp10M1", 1.0, 0.0, "9e+0", false, false); + yield return new("Exp10M1", 3.25, 0.0, "1.77727941003892280122542119519268484473579052640225535801183072277630188153949380490030039927870215508827904815953580779e+3", false, false); + // Exp10M1: dense-and-sparse-low + yield return new("Exp10M1", -1.0, 5.551115123125783e-17, "-8.99999999999999987218085067996765462754319617590916607329559738303097631479573890845779454801469327857031874845490804547e-1", false, false); + yield return new("Exp10M1", -1.0, -5.551115123125783e-17, "-9.00000000000000012781914932003232903472187092736774381834984894590193923691417843226363938571113470345052880638874691579e-1", false, false); + yield return new("Exp10M1", -1.0, 5e-324, "-9.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-1", false, false); + yield return new("Exp10M1", -1.0, -5e-324, "-9.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-1", false, false); + yield return new("Exp10M1", 1.0, 5.551115123125783e-17, "9.00000000000000127819149320032345372456803824090833926704402616969023685204261091542205451985306721429681251545091954530e+0", false, false); + yield return new("Exp10M1", 1.0, -5.551115123125783e-17, "8.99999999999999872180850679967670965278129072632256181650151054098060763085821567736360614288865296549471193611253084211e+0", false, false); + yield return new("Exp10M1", 1.0, 5e-324, "9.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp10M1", 1.0, -5e-324, "9.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + // Exp10M1: tiny + yield return new("Exp10M1", -1e-20, 0.0, "-2.30258509299404555770222241510393875598524143224381069825640182231193232027551245827252136404378259029970201038081656945e-20", false, false); + yield return new("Exp10M1", 1e-20, 0.0, "2.30258509299404555775524139620872273027502276727279330790399614334475670162688907681242268603731733932824173589898447578e-20", false, false); + yield return new("Exp10M1", -1e-100, 0.0, "-2.30258509299404573005104192078836854051056555890772672956315725895733581426274915955308976822915813062782353584671664263e-100", false, false); + yield return new("Exp10M1", 1e-100, 0.0, "2.30258509299404573005104192078836854051056555890772672956315725895733581426274915955308976822915813115801334689455646489e-100", false, false); + yield return new("Exp10M1", -1e-300, 0.0, "-2.30258509299404574171868275840485705189541108963464621458993298511568529915114775677970158135301609867111333610627501526e-300", false, false); + yield return new("Exp10M1", 1e-300, 0.0, "2.30258509299404574171868275840485705189541108963464621458993298511568529915114775677970158135301609867111333610627501526e-300", false, false); + yield return new("Exp10M1", -2.2250738585072014e-308, 0.0, "-5.12342189740942434501569327901461575239505617669057681130375993183183004105261803650543347093083043244937774062452497400e-308", false, false); + yield return new("Exp10M1", 2.2250738585072014e-308, 0.0, "5.12342189740942434501569327901461575239505617669057681130375993183183004105261803650543347093083043244937774062452497400e-308", false, false); + yield return new("Exp10M1", -1e-323, 0.0, "-2.27525638214905982837507870853537983552919151750274587233835884873293322219317237853594514215541907751178107438528862527e-323", false, false); + yield return new("Exp10M1", 1e-323, 0.0, "2.27525638214905982837507870853537983552919151750274587233835884873293322219317237853594514215541907751178107438528862527e-323", false, false); + yield return new("Exp10M1", -5e-324, 0.0, "-1.13762819107452991418753935426768991776459575875137293616917942436646661109658618926797257107770953875589053719264431264e-323", false, false); + yield return new("Exp10M1", 5e-324, 0.0, "1.13762819107452991418753935426768991776459575875137293616917942436646661109658618926797257107770953875589053719264431264e-323", false, false); + // Exp10M1: absolute-log-half-switch + yield return new("Exp10M1", -0.2171472409516259, -5.491598251083826e-18, "-3.93469340287366576396200465008820362882072578421007457645958090117166064388468963572181401431536248330076974888414600998e-1", false, false); + yield return new("Exp10M1", -0.2171472409516259, -5.491598251083825e-18, "-3.93469340287366576396200465008819286989462006654938071200708034981872180311176975058969796694858078679858192326373506407e-1", false, false); + yield return new("Exp10M1", -0.2171472409516259, -5.4915982510838246e-18, "-3.93469340287366576396200465008818211096851434888868684755457979844669827319829825582055598674633324019908386971626983556e-1", false, false); + yield return new("Exp10M1", 0.2171472409516259, 5.4915982510838246e-18, "6.48721270700128146848650787814159941493780812305632348189724146484756143871422484272097497594536517379427525706613712148e-1", false, false); + yield return new("Exp10M1", 0.2171472409516259, 5.491598251083825e-18, "6.48721270700128146848650787814162866073113502901405270661424684716022537920234692762763311713810662847016175115314168844e-1", false, false); + yield return new("Exp10M1", 0.2171472409516259, 5.491598251083826e-18, "6.48721270700128146848650787814165790652446193497178193133125222952476688338302012652507096358970116413935923203943027614e-1", false, false); + // Exp10M1: input-high-series-switch + yield return new("Exp10M1", -0.12500000000000003, 0.0, "-2.50105790667544220623235682831241227091366353848874384997791525723347737745979450739687740450604812196950274788445358612e-1", false, false); + yield return new("Exp10M1", -0.12499999999999999, 0.0, "-2.50105790667544148735105745161673580206563090718922547131556925027730093155363334658244471038575458141695010455686722771e-1", false, false); + yield return new("Exp10M1", -0.125, -6.938893903907228e-18, "-2.50105790667544184679170713996458265087678613615781682045379956151380255522669580317417278415513574852965788877157596059e-1", false, false); + yield return new("Exp10M1", -0.125, 6.938893903907228e-18, "-2.50105790667544160716460734773268666597760240869410293263525178017542720692692645153465620057550733979128107066362824524e-1", false, false); + // Exp10M1: quadratic-series-switch + yield return new("Exp10M1", -2.4108186663832177e-17, -6.096898820344505e-34, "-5.55111512312578254804376278947525111135952916952488118057090928416625710539857409417218238146094275651639899718827039890e-17", false, false); + yield return new("Exp10M1", -2.4108186663832177e-17, -6.0968988203445044e-34, "-5.55111512312578254804376278947523141770135784523200699580201696017754023733808975824261769382695015710450861125941702211e-17", false, false); + yield return new("Exp10M1", -2.4108186663832177e-17, -6.096898820344504e-34, "-5.55111512312578254804376278947521172404318652093913281103312463618882336927760541843465128450327642975730308116022193819e-17", false, false); + // Exp10M1: input-high-series-switch + yield return new("Exp10M1", 0.12499999999999999, 0.0, "3.33521432163323983063537950508510118386148151192037635826279917554245850263297432726335583118584038111090819342091499132e-1", false, false); + yield return new("Exp10M1", 0.12500000000000003, 0.0, "3.33521432163324110900719244868977125484648720885341763163649432278606530942581178394286995865756427602775067340009630110e-1", false, false); + yield return new("Exp10M1", 0.125, -6.938893903907228e-18, "3.33521432163324004369734832901920435192160510937529511447041291879436130189385376993868097672364644153343456779246847864e-1", false, false); + yield return new("Exp10M1", 0.125, 6.938893903907228e-18, "3.33521432163324046982128597688742090056670512672567520458135576124210663897996160077897850952474308729078332610355280926e-1", false, false); + // Exp10M1: quadratic-series-switch + yield return new("Exp10M1", 2.4108186663832177e-17, 6.096898820344504e-34, "5.55111512312578285619255389143294821299965733452486897471468053264121209511768939103939182765135013753615105042745685340e-17", false, false); + yield return new("Exp10M1", 2.4108186663832177e-17, 6.0968988203445044e-34, "5.55111512312578285619255389143296790665782865881992959475766301532431147684229694003121733250823665487940329649836608834e-17", false, false); + yield return new("Exp10M1", 2.4108186663832177e-17, 6.096898820344505e-34, "5.55111512312578285619255389143298760031599998311499021480064549800741085856690449290144455905480473074705970331287834858e-17", false, false); + // Exp10M1: random-moderate + yield return new("Exp10M1", -2.9677357165710085, -1.1102230246251565e-16, "-9.98922879520165034714208159852150757557324583422419498005448293567172081929540244921649768419178178105615911287043460489e-1", false, false); + yield return new("Exp10M1", 0.5995899657426765, -2.7755575615628914e-17, "2.97731479519890352879887600772233924298609686175885392875139750363632580171965156242522363738652933349600506618119299140e+0", false, false); + yield return new("Exp10M1", 0.8836851964549527, -2.7755575615628914e-17, "6.65041856144971858829962647061748349456294701097773926744045197071139270643977619590118903694461196808465346045233138263e+0", false, false); + yield return new("Exp10M1", -0.89395009436352, -2.7755575615628914e-17, "-8.72341450465680217809052941648826847479618841723966114669093459382270793548409479815284018448824046387929644197921673353e-1", false, false); + yield return new("Exp10M1", 2.4580147142578235, -1.1102230246251565e-16, "2.86087784810333896489677828661765156448804286447564220961070487382944745109858177125943722517534602403829765495652338124e+2", false, false); + yield return new("Exp10M1", 1.3118516264245281, 5.551115123125783e-17, "1.95046153376223160193906750616367118569485139159211538092010333804564633619337531743249262303443651659625057726925490890e+1", false, false); + // Exp10M1: random-range + yield return new("Exp10M1", 120.58969150830279, 3.552713678800501e-15, "3.88768893311729066468034240903516540867755441511370722931132770507875572390376800025883066808459259868657675470129694955e+120", false, false); + yield return new("Exp10M1", 128.81118142762747, -7.105427357601002e-15, "6.47413017703446171983916778641885806651744797264667468232608054450205762476305583291383873034783101822647556402661401770e+128", false, false); + yield return new("Exp10M1", 42.94992135127048, 1.7763568394002505e-15, "8.91089551309232685378712887186962248614334667700386757268709525041370952477415725129562497023716824045224543686275322633e+42", false, false); + yield return new("Exp10M1", 116.81408221381518, -3.552713678800501e-15, "6.51751761797928744146122574992738336073782813915351110967844224366907044327161884077729178445278842875234452038465776427e+116", false, false); + yield return new("Exp10M1", 228.96449953361966, -7.105427357601002e-15, "9.21508899144045996730372464171725339874292277767159122312745668742022577561608598769932754185213295353703683516203829072e+228", false, false); + yield return new("Exp10M1", 7.524609271856807, 2.220446049250313e-16, "3.34664201369735369097716632117090216287459757705675042241703792408934408612576386525573391905086740939515500212837287795e+7", false, false); + // Exp10M1: integer-power + yield return new("Exp10M1", -324.0, 0.0, "-1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp10M1", -323.0, 0.0, "-1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp10M1", -308.0, 0.0, "-1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+0", false, false); + yield return new("Exp10M1", 22.0, 0.0, "9.999999999999999999999e+21", false, false); + yield return new("Exp10M1", 308.0, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+308", false, false); + yield return new("Exp10M1", 309.0, 0.0, "1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+309", true, false); + // Exp10M1: negative-saturation + yield return new("Exp10M1", -20.0, 0.0, "-9.9999999999999999999e-1", false, false); + yield return new("Exp10M1", -100.0, 0.0, "-9.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999e-1", false, false); + // Exp10M1: overflow-adjacent-low + yield return new("Exp10M1", 308.25471555991675, -2.895125789515847e-15, "1.79769313486231580793728971405135339792924520420229515147549215774814264965213204794170642362509327672779802464860968196e+308", false, false); + yield return new("Exp10M1", 308.25471555991675, -2.8951257895158466e-15, "1.79769313486231580793728971405298608023672269474350236922097481966929437654805569985823689453925527448616186188699371312e+308", false, false); + yield return new("Exp10M1", 308.25471555991675, -2.8951257895158462e-15, "1.79769313486231580793728971405461876254420018528470958696645896440832277858828718822634729443205898882595653371847623943e+308", true, false); + } +} diff --git a/1-tests/Just.PreciseMath.Tests/ReferenceData/generate_exponential_functions.py b/1-tests/Just.PreciseMath.Tests/ReferenceData/generate_exponential_functions.py new file mode 100644 index 0000000..909f6ab --- /dev/null +++ b/1-tests/Just.PreciseMath.Tests/ReferenceData/generate_exponential_functions.py @@ -0,0 +1,270 @@ +"""Independent stdlib references for the five new exponential-family functions. + +Run: python3 path/to/generate_exponential_functions.py [--check] +Exact binary64 high+low sums are formed at 2200 digits and checked with Fraction. +Decimal ln/exp (or a cancellation-safe expm1 series) runs at 450 and 650 digits; +all inputs, flags, and 120-significant-digit references must agree. Integer base-2 +and base-10 powers use exact rational arithmetic, also checked against ln/exp. +No production arithmetic, binary64 transcendental functions, or rounded DD +constants are used. Zeros, nonfinite inputs, and enormous inputs belong in the +main suite. These rows specify an error bound, not correctly rounded DD results. +""" +from collections import Counter +from decimal import Decimal, localcontext +from fractions import Fraction +from pathlib import Path +import argparse +import math +import random + + +OUTPUT = Path(__file__).with_name('ExponentialFunctionsReferenceData.cs') +OPERATIONS = {'Exp2': 2, 'Exp10': 10, 'ExpM1': None, 'Exp2M1': 2, 'Exp10M1': 10} +EPSILON = math.ulp(0.0) + + +def exact_sum(high, low): + with localcontext() as context: + context.prec = 2200 + value = Decimal.from_float(high) + Decimal.from_float(low) + assert Fraction(value) == Fraction(high) + Fraction(low), (high, low) + return value + + +def exact_decimal(value): + with localcontext() as context: + context.prec = 2200 + result = Decimal(value.numerator) / Decimal(value.denominator) + assert Fraction(result) == value + return result + + +OVERFLOW = exact_decimal(Fraction(2) ** 1024 - Fraction(2) ** 970) +HALF_EPSILON = exact_decimal(Fraction(2) ** -1075) + + +def expm1(value): + # exp(value)-1 would erase tiny inputs even with hundreds of digits. For + # |value| <= 0.5, sum x + x^2/2! + ... until the context stops changing. + if abs(value) > Decimal('0.5'): + return value.exp() - 1 + total = term = value + for denominator in range(2, 10000): + term = term * value / denominator + updated = total + term + if updated == total: + return updated + total = updated + raise AssertionError('expm1 series did not converge') + + +def neighbors(boundary): + high = float(boundary) + low = float(boundary - Decimal.from_float(high)) + adjacent = [math.nextafter(low, -math.inf), low, math.nextafter(low, math.inf)] + assert exact_sum(high, adjacent[0]) < boundary < exact_sum(high, adjacent[-1]) + assert math.nextafter(adjacent[0], math.inf) == adjacent[1] + assert math.nextafter(adjacent[1], math.inf) == adjacent[2] + return [(high, residual) for residual in adjacent] + + +def inputs(operation, log_base): + base = OPERATIONS[operation] + minus_one = operation.endswith('M1') + for high in [-10.25, -2.0, -1.0, -0.125, 0.125, 0.75, 1.0, 3.25]: + yield 'ordinary', high, 0.0 + for high in [-1.0, 1.0]: + for low in [math.ulp(high) / 4, -math.ulp(high) / 4, EPSILON, -EPSILON]: + yield 'dense-and-sparse-low', high, low + for magnitude in [1e-20, 1e-100, 1e-300, math.ldexp(1.0, -1022), 2 * EPSILON, EPSILON]: + for sign in [-1, 1]: + yield 'tiny', sign * magnitude, 0.0 + for sign in [-1, 1]: + for high, low in neighbors(Decimal(sign) / (2 * log_base)): + yield 'absolute-log-half-switch', high, low + if minus_one: + # Concrete dispatch uses the input high (0.5 for e/2, 0.125 for 10). + # Include both neighboring highs and lows on either side of the anchor. + switch = 0.125 if base == 10 else 0.5 + for sign in [-1, 1]: + high = sign * switch + for adjacent in [math.nextafter(high, -math.inf), high, math.nextafter(high, math.inf)]: + yield 'input-high-series-switch', adjacent, 0.0 + for low in [-math.ulp(high) / 4, math.ulp(high) / 4]: + yield 'input-high-series-switch', high, low + for high, low in neighbors(Decimal(sign) * Decimal(2) ** -54 / log_base): + yield 'quadratic-series-switch', high, low + if operation == 'Exp2': + for sign in [-1, 1]: + low = sign * math.ldexp(1.0, -500) + for adjacent in [math.nextafter(low, -math.inf), low, math.nextafter(low, math.inf)]: + yield 'sparse-correction-switch', 1000.0, adjacent + rng = random.Random(20260916 + (base or 1)) + # Both local/moderate arguments and the complete useful exponential output + # range, expressed in natural-log coordinates, with reproducible dense lows. + for label, lower, upper in [('random-moderate', -8, 8), ('random-range', -748, 712)]: + for _ in range(6): + high = float(Decimal.from_float(rng.uniform(lower, upper)) / log_base) + yield label, high, rng.choice([-1, 1]) * math.ulp(high) / 4 + if base == 2: + for exponent in [-1076, -1075, -1074, -1022, -100, -1, 1, 10, 53, 100, 1023, 1024]: + yield 'integer-power', float(exponent), 0.0 + elif base == 10: + for exponent in [-324, -323, -308, -1, 1, 22, 308, 309]: + yield 'integer-power', float(exponent), 0.0 + if minus_one: + for high in ([-40.0, -100.0, -745.0, -746.0] if base is None else + [-54.0, -100.0, -1075.0] if base == 2 else [-20.0, -100.0, -324.0]): + yield 'negative-saturation', high, 0.0 + if base == 2: + # The exponential correction is exactly half epsilon at -1075, + # but exp2m1 is near -1, NOT an underflowing result. + for low in [-EPSILON, EPSILON]: + yield 'negative-saturation-sparse-low', -1075.0, low + # For M1, exp(x ln b) must reach OVERFLOW + 1, not OVERFLOW. + overflow_log = (OVERFLOW + int(minus_one)).ln() / log_base + for high, low in neighbors(overflow_log): + yield 'overflow-adjacent-low', high, low + if not minus_one: + # Base 2 has the exactly attainable threshold x=-1075; do not obtain + # it by an inexact quotient of two rounded logarithms. + underflow_log = Decimal(-1075) if base == 2 else HALF_EPSILON.ln() / log_base + for high, low in neighbors(underflow_log): + yield 'underflow-adjacent-low', high, low + normal_log = Decimal(-1022) if base == 2 else Decimal(2).ln() * -1022 / log_base + for high, low in neighbors(normal_log): + yield 'min-normal-adjacent-low', high, low + + +def evaluate(operation, high, low, log_base): + value = exact_sum(high, low) + base = OPERATIONS[operation] + minus_one = operation.endswith('M1') + argument = value * log_base + result = expm1(argument) if minus_one else argument.exp() + if base is not None and value == value.to_integral_value(): + rational = Fraction(base) ** int(value) - int(minus_one) + exact = exact_decimal(rational) + # Cross-check integer fixtures by an independent mathematical identity. + assert abs(result - exact) <= abs(exact) * Decimal('1e-400') + result = exact + # Use exact input versus the independently computed thresholds: comparing a + # rounded exponential with half epsilon can misclassify an exact midpoint. + overflow = value >= (OVERFLOW + int(minus_one)).ln() / log_base + if minus_one: + underflow = abs(result) <= HALF_EPSILON + assert (result < 0) == (value < 0) + else: + threshold = Decimal(-1075) if base == 2 else HALF_EPSILON.ln() / log_base + underflow = value <= threshold + with localcontext() as output_context: + output_context.prec = 120 + rounded = +result + assert abs(result - rounded) <= abs(result) * Decimal(2) ** -350 + return format(rounded, 'e'), overflow, underflow + + +def generate_at_precision(precision): + rows = [] + counts = Counter() + with localcontext() as context: + context.prec = precision + for operation, base in OPERATIONS.items(): + log_base = Decimal(base).ln() if base else Decimal(1) + seen = set() + boundary_flags = {} + for label, high, low in inputs(operation, log_base): + assert math.isfinite(high) and math.isfinite(low) + assert high != 0 and abs(low) <= math.ulp(high) / 2 + assert float(Fraction(high) + Fraction(low)) == high, (high, low) + result = evaluate(operation, high, low, log_base) + if label in ['overflow-adjacent-low', 'underflow-adjacent-low']: + flag_index = 1 if label.startswith('overflow') else 2 + boundary_flags.setdefault(label, []).append(result[flag_index]) + key = (high, low) + if key in seen: + continue + seen.add(key) + rows.append((operation, label, high, low, *result)) + counts[operation] += 1 + assert boundary_flags['overflow-adjacent-low'][0] is False + assert boundary_flags['overflow-adjacent-low'][-1] is True + if not operation.endswith('M1'): + assert boundary_flags['underflow-adjacent-low'][0] is True + assert boundary_flags['underflow-adjacent-low'][-1] is False + return rows, counts + + +def verify_sparse_splits(precision): + # Component-retention sentinels in DoubleDoubleExponentialFunctionsTests: + # unlike the 120-digit rows, these retain the tiny correction beside 2^n. + with localcontext() as context: + context.prec = precision + ln2 = Decimal(2).ln() + for exponent, expected_low in [(1, EPSILON), (500, 1.1210060331144859e-173), + (1000, 3.6694906201918696e-23)]: + for sign in [-1, 1]: + value = exact_sum(float(exponent), sign * EPSILON) + result = (value * ln2).exp() + high = math.ldexp(1.0, exponent) + assert float(result) == high + assert float(result - Decimal.from_float(high)) == sign * expected_low + + +def generate(): + first, counts = generate_at_precision(450) + second, second_counts = generate_at_precision(650) + assert first == second and counts == second_counts, 'Precision stability check failed' + verify_sparse_splits(450) + verify_sparse_splits(650) + assert 200 <= len(first) <= 400 + lines = [] + previous_group = None + for operation, label, high, low, reference, overflow, underflow in first: + group = operation, label + if group != previous_group: + lines.append(f' // {operation}: {label}') + previous_group = group + lines.append(f' yield return new("{operation}", {high!r}, {low!r}, ' + f'"{reference}", {str(overflow).lower()}, {str(underflow).lower()});') + content = '''// Generated by generate_exponential_functions.py; do not hand-edit. +// Exact binary64 sums verified with Fraction at 2200 digits; Decimal ln/exp and +// cancellation-safe expm1 at 450/650 digits; 120-digit references agree. +// Integer base-2/base-10 cases additionally use exact rational powers. +// Row: operation, high, low, reference, overflow, underflow. +// Overflow: exact result >= 2^1024 - 2^970 (binary64 overflow midpoint). +// Underflow: |exact result| <= 2^-1075, NOT exp(x) <= 2^-1075 for M1. +// Compare exact component sums with relative tolerance 2^-100 + double.Epsilon +// absolute; reference rounding uncertainty is bounded by 2^-350 relative. +// Tiny corrections to 1/-1 below 120 digits are not component-retention oracles. +// Signed zeros, infinities, NaNs, and enormous arguments are tested separately. +using Xunit; + +namespace Just.PreciseMath.Tests.ReferenceData; + +internal static class ExponentialFunctionsReferenceData +{ + internal static IEnumerable> Cases() + { +''' + '\n'.join(lines) + '\n }\n}\n' + return content, counts + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument('--check', action='store_true', help='Verify fixtures without writing') + args = parser.parse_args() + content, counts = generate() + if args.check: + assert OUTPUT.read_text() == content, 'Fixture is stale; regenerate it' + action = 'Verified' + else: + OUTPUT.write_text(content) + action = 'Generated' + print(f'{action} {sum(counts.values())} references at 450/650 digits: {dict(counts)}') + print('Exact 2200-digit sums, rational integer powers, adjacent-low boundary flags, ' + 'normalized pairs, and reference rounding bounds verified.') + + +if __name__ == '__main__': + main() diff --git a/README.md b/README.md index a154440..2e34be6 100644 --- a/README.md +++ b/README.md @@ -93,7 +93,11 @@ constants. The natural logarithm function is provided separately by `DDMath.Log` Root functions are implemented on `DoubleDouble` in `DoubleDouble.RootFunctions.cs`: `Sqrt`, `Cbrt`, `Hypot`, and `RootN` complete `IRootFunctions`. The additional `InvSqrt` helper and three-argument `Hypot` overload live alongside -them. `DDMath` exposes thin forwarding wrappers; the other mathematical kernels +them. Exponential functions live in `DoubleDouble.ExponentialFunctions.cs`: +`Exp`, `Exp2`, `Exp10`, `ExpM1`, `Exp2M1`, and `Exp10M1` complete +`IExponentialFunctions`, including cancellation-safe overrides of +the interface's default minus-one methods. `DDMath` exposes thin forwarding +wrappers for roots and exponentials; logarithm, power, and reciprocal kernels remain in `DDMath`. The `DDMath` static class provides: @@ -163,6 +167,20 @@ The `DDMath` static class provides: components affect range boundaries; representable subnormals are retained. Either zero maps to one, negative infinity to positive zero, positive infinity to positive infinity, and NaN to canonical NaN. +- `Exp2(DoubleDouble)`: reduces the argument in base two before evaluating a + bounded natural exponential. Integer powers from -1074 through 1023 are exact. + An exponent-tracked correction preserves sparse lows that would otherwise round + prematurely before final scaling. The exact half-minimum-subnormal threshold + uses both input components. Special-value behavior is the same as `Exp`. +- `Exp10(DoubleDouble)`: subtracts three split `log10(2)` products before converting + the bounded remainder to a natural exponent, avoiding large-argument amplification + of a rounded `ln(10)` product. Special-value behavior is the same as `Exp`. +- `ExpM1(DoubleDouble)`, `Exp2M1(DoubleDouble)`, and `Exp10M1(DoubleDouble)`: + compute the corresponding exponential minus one, using a direct series near + zero to avoid cancellation. Signed zeros are preserved, negative infinity maps + to negative one, positive infinity to positive infinity, and NaN to canonical NaN. + Their tested relative error is measured against the minus-one result itself, + not against the exponential before subtraction. - `Log(DoubleDouble)`: natural logarithm using binary range reduction and a centered atanh series. The bounded mantissa avoids denominator overflow for large inputs; a separate near-one path retains even minimum-subnormal low @@ -181,7 +199,10 @@ DoubleDouble fifthRoot = DoubleDouble.RootN(new DoubleDouble(2.0), 5); DoubleDouble reciprocal = DDMath.Reciprocal(new DoubleDouble(3.0)); DoubleDouble magnitude = DDMath.Abs(-root); DoubleDouble smallPower = DDMath.Pow(new DoubleDouble(2.0), -1024); -DoubleDouble exponential = DDMath.Exp(new DoubleDouble(1.0)); +DoubleDouble exponential = DoubleDouble.Exp(new DoubleDouble(1.0)); +DoubleDouble binaryPower = DoubleDouble.Exp2(new DoubleDouble(-1000.0)); +DoubleDouble decimalPower = DoubleDouble.Exp10(new DoubleDouble(0.5)); +DoubleDouble tinyChange = DoubleDouble.ExpM1(new DoubleDouble(1e-30)); DoubleDouble logarithm = DDMath.Log(new DoubleDouble(10.0)); DoubleDouble fractionalPower = DDMath.Pow(new DoubleDouble(2.0), 0.5); DoubleDouble preciseExponent = DoubleDouble.FromComponents(0.5, 1e-30); @@ -225,13 +246,17 @@ additional tests cover extreme magnitudes, near-one cancellation, sparse lows of either sign, and range-reduction transitions. This is sampled approximate accuracy, not a universal error proof or a correct-rounding guarantee. -Exponential tests use independent high-precision decimal references evaluated at +Exponential-family tests use independent high-precision decimal references evaluated at two precisions, with exact integer comparisons of the stored component sum. They check `2^-100` relative error plus one minimum binary64 subnormal, with a separately -bounded reference-rounding allowance. Integer-power tests use exact rational +bounded reference-rounding allowance. The base-two/base-ten and minus-one fixtures +are checked at 450/650 digits, including adjacent lows at true range boundaries, +series transitions, and tiny inputs. Additional component checks pin sparse +base-two corrections; constrained generic calls and facade calls match direct +calls bit-for-bit. Integer-power tests use exact rational references and high-precision fixtures for large exponents; their tested absolute error bound is `|exact result| * (|exponent| + 1) * 2^-100 + 2^-1074`, not uniform -relative accuracy independent of the exponent. Both functions are approximate; +relative accuracy independent of the exponent. These functions are approximate; precision decreases near underflow and approximation can affect results extremely close to a rounding boundary. Final scaling uses the existing allocating exact boundary machinery where necessary. No performance measurements are claimed. @@ -342,9 +367,10 @@ bool success = DoubleDouble.TryParse("1.25e-2".AsSpan(), CultureInfo.InvariantCu ## Deferred scope -Natural `DDMath.Log`, `Exp`, and all three `Pow` overloads are implemented. +Natural `DDMath.Log`, the complete exponential family on `DoubleDouble`, and all +three `DDMath.Pow` overloads are implemented. Logarithms in other bases, generic-math interfaces beyond `ISignedNumber`, -`IFloatingPointConstants`, and `IRootFunctions`, additional text formats/general +`IFloatingPointConstants`, `IRootFunctions`, and `IExponentialFunctions`, additional text formats/general round-trip formatting, and non-arithmetic performance benchmarks remain deferred. Replacing allocating arithmetic boundary fallbacks is also deferred; the current `BigInteger` paths remain in place. That optimization does not require removing