225 lines
9.7 KiB
C#
225 lines
9.7 KiB
C#
using System.Globalization;
|
|
using Just.PreciseMath.Tests.ReferenceData;
|
|
|
|
namespace Just.PreciseMath.Tests;
|
|
|
|
public class DoubleDoubleExponentialFunctionsTests
|
|
{
|
|
public static IEnumerable<TheoryDataRow<string, double, double, string, bool, bool>> 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<T>(string operation, T value) where T : IExponentialFunctions<T>
|
|
{
|
|
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));
|
|
}
|
|
}
|