241 lines
10 KiB
C#
241 lines
10 KiB
C#
using System.Globalization;
|
|
using Just.PreciseMath.Tests.ReferenceData;
|
|
|
|
namespace Just.PreciseMath.Tests;
|
|
|
|
public class DoubleDoubleRootNTests
|
|
{
|
|
public static IEnumerable<TheoryDataRow<double, double, int, string>> ReferenceCases => RootNReferenceData.Cases();
|
|
|
|
[Theory]
|
|
[MemberData(nameof(ReferenceCases))]
|
|
public void LargeDegreesMeetIndependentHighPrecisionReferences(double high, double low, int degree, string reference)
|
|
{
|
|
DoubleDouble input = DoubleDouble.FromComponents(high, low);
|
|
(Units(input.High) + Units(input.Low)).ShouldBe(Units(high) + Units(low));
|
|
DoubleDouble actual = DoubleDouble.RootN(input, degree);
|
|
DoubleDouble.IsFinite(actual).ShouldBeTrue();
|
|
DoubleDouble.IsCanonical(actual).ShouldBeTrue();
|
|
Math.Sign(actual.High).ShouldBe(Math.Sign(high));
|
|
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));
|
|
// 2^-100 relative plus the separately bounded 2^-350 relative reference
|
|
// rounding allowance. For |degree|>=2 every nonzero root is normal.
|
|
BigInteger magnitude = BigInteger.Abs(numerator);
|
|
BigInteger bound = (magnitude << 1324) + (magnitude << 1074);
|
|
((error << 350) <= bound).ShouldBeTrue(
|
|
$"RootN reference failed for ({high:R}, {low:R}), {degree}: ({actual.High:R}, {actual.Low:R})");
|
|
}
|
|
|
|
[Fact]
|
|
public void SpecialValuesMatchRuntimeRootNSemantics()
|
|
{
|
|
double[] values = [0.0, -0.0, 1.0, -1.0, double.PositiveInfinity, double.NegativeInfinity, double.NaN];
|
|
int[] degrees = [int.MinValue, int.MinValue + 1, -1000, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 1000, int.MaxValue];
|
|
foreach (double value in values)
|
|
{
|
|
foreach (int degree in degrees)
|
|
{
|
|
AssertBits(new DoubleDouble(double.RootN(value, degree)), DoubleDouble.RootN(new DoubleDouble(value), degree));
|
|
}
|
|
}
|
|
foreach (int degree in new[] { int.MinValue, -4, -2, 0, 2, 4, int.MaxValue - 1 })
|
|
{
|
|
DoubleDouble.IsNaN(DoubleDouble.RootN(new DoubleDouble(-2.0), degree)).ShouldBeTrue();
|
|
}
|
|
}
|
|
|
|
[Fact]
|
|
public void IdentitiesAndExistingKernelsKeepTheirComponentBits()
|
|
{
|
|
DoubleDouble[] values = [new(double.Epsilon), new(double.MaxValue), new(2.0),
|
|
DoubleDouble.FromComponents(1.0, double.Epsilon), DoubleDouble.FromComponents(1.0, -Math.ScaleB(1.0, -54))];
|
|
foreach (DoubleDouble value in values)
|
|
{
|
|
AssertBits(value, DoubleDouble.RootN(value, 1));
|
|
AssertBits(-value, DoubleDouble.RootN(-value, 1));
|
|
AssertBits(DDMath.Reciprocal(value), DoubleDouble.RootN(value, -1));
|
|
AssertBits(DoubleDouble.Sqrt(value), DoubleDouble.RootN(value, 2));
|
|
AssertBits(DoubleDouble.InvSqrt(value), DoubleDouble.RootN(value, -2));
|
|
AssertBits(DoubleDouble.Cbrt(value), DoubleDouble.RootN(value, 3));
|
|
AssertBits(1.0 / DoubleDouble.Cbrt(value), DoubleDouble.RootN(value, -3));
|
|
}
|
|
}
|
|
|
|
[Fact]
|
|
public void SmallDegreesMeetIndependentExactPowerInequalities()
|
|
{
|
|
Random random = new(8675309);
|
|
for (int exponent = -1074; exponent <= 1023; ++exponent)
|
|
{
|
|
double high = Math.ScaleB(1.0 + random.NextDouble(), exponent);
|
|
foreach (int degree in new[] { -5, -4, -3, -2, 2, 3, 4, 5 })
|
|
{
|
|
DoubleDouble input = DoubleDouble.FromComponents(high, Math.ScaleB((degree < 0 ? -1.0 : 1.0), exponent - 54));
|
|
AssertBound(input, DoubleDouble.RootN(input, degree), degree);
|
|
if ((degree & 1) != 0)
|
|
{
|
|
AssertBits(-DoubleDouble.RootN(input, degree), DoubleDouble.RootN(-input, degree));
|
|
}
|
|
}
|
|
}
|
|
foreach (double high in new[] { double.Epsilon, 1e-308, 0.5, 1.0, 2.0, 81.0, 1e308, double.MaxValue })
|
|
{
|
|
foreach (int degree in new[] { -31, -17, -7, -6, 6, 7, 17, 31 })
|
|
{
|
|
DoubleDouble input = new(high);
|
|
AssertBound(input, DoubleDouble.RootN(input, degree), degree);
|
|
}
|
|
}
|
|
}
|
|
|
|
[Fact]
|
|
public void ExactPowerOfTwoRootsAndSparseCorrectionsSurvive()
|
|
{
|
|
foreach (int degree in new[] { -31, -7, -5, -4, 4, 5, 7, 31 })
|
|
{
|
|
for (int exponent = -1074; exponent <= 1023; ++exponent)
|
|
{
|
|
if (exponent % degree == 0)
|
|
{
|
|
AssertBits(new DoubleDouble(Math.ScaleB(1.0, exponent / degree)),
|
|
DoubleDouble.RootN(new DoubleDouble(Math.ScaleB(1.0, exponent)), degree));
|
|
}
|
|
}
|
|
}
|
|
foreach (int degree in new[] { int.MinValue, int.MinValue + 1, -7, -4, 4, 7, int.MaxValue })
|
|
{
|
|
foreach (int exponent in new[] { -100, -500, -1000, -1070 })
|
|
{
|
|
foreach (double sign in new[] { -1.0, 1.0 })
|
|
{
|
|
double low = Math.ScaleB(sign, exponent);
|
|
// (1+d)^(1/n) = 1+d/n+O(d²); the quadratic term is
|
|
// below half an ulp of this rounded correction.
|
|
AssertBits(DoubleDouble.FromComponents(1.0, low / degree),
|
|
DoubleDouble.RootN(DoubleDouble.FromComponents(1.0, low), degree));
|
|
}
|
|
}
|
|
}
|
|
foreach (double sign in new[] { -1.0, 1.0 })
|
|
{
|
|
DoubleDouble input = DoubleDouble.FromComponents(Math.ScaleB(1.0, 900), Math.ScaleB(sign, -174));
|
|
// d*r/(n*high), with exact r=2^225 and n=4, is sign*2^-851.
|
|
AssertBits(DoubleDouble.FromComponents(Math.ScaleB(1.0, 225), Math.ScaleB(sign, -851)), DoubleDouble.RootN(input, 4));
|
|
}
|
|
}
|
|
|
|
[Fact]
|
|
public void NegativeDegreeRetainsSparseCorrectionAfterLargeRescaling()
|
|
{
|
|
// (2^-900+d)^(-1/4) = 2^225 - d*2^225/(4*2^-900) + O(d²).
|
|
// For d=±2^-1074, the correction is exactly ∓2^49 at binary64
|
|
// residual precision; the quadratic contribution is below half an ulp.
|
|
foreach (double sign in new[] { -1.0, 1.0 })
|
|
{
|
|
DoubleDouble input = DoubleDouble.FromComponents(Math.ScaleB(1.0, -900), sign * double.Epsilon);
|
|
AssertBits(DoubleDouble.FromComponents(Math.ScaleB(1.0, 225), Math.ScaleB(-sign, 49)), DoubleDouble.RootN(input, -4));
|
|
}
|
|
}
|
|
|
|
[Fact]
|
|
public void GenericContractAndFacadeDispatchToTypeMembers()
|
|
{
|
|
DoubleDouble[] values = [DoubleDouble.Zero, DoubleDouble.NegativeZero, DoubleDouble.NaN,
|
|
DoubleDouble.PositiveInfinity, DoubleDouble.NegativeInfinity, new(-2.0), new(double.Epsilon),
|
|
new(double.MaxValue), DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -54))];
|
|
foreach (DoubleDouble value in values)
|
|
{
|
|
AssertBits(DoubleDouble.Sqrt(value), GenericSqrt(value));
|
|
AssertBits(DoubleDouble.Cbrt(value), GenericCbrt(value));
|
|
AssertBits(DoubleDouble.Hypot(value, DoubleDouble.One), GenericHypot(value, DoubleDouble.One));
|
|
foreach (int degree in new[] { int.MinValue, -7, -3, -2, -1, 0, 1, 2, 3, 4, 7, int.MaxValue })
|
|
{
|
|
DoubleDouble result = DoubleDouble.RootN(value, degree);
|
|
AssertBits(result, GenericRootN(value, degree));
|
|
AssertBits(result, DDMath.RootN(value, degree));
|
|
}
|
|
}
|
|
}
|
|
|
|
private static T GenericSqrt<T>(T value) where T : IRootFunctions<T>
|
|
{
|
|
return T.Sqrt(value);
|
|
}
|
|
|
|
private static T GenericCbrt<T>(T value) where T : IRootFunctions<T>
|
|
{
|
|
return T.Cbrt(value);
|
|
}
|
|
|
|
private static T GenericHypot<T>(T x, T y) where T : IRootFunctions<T>
|
|
{
|
|
return T.Hypot(x, y);
|
|
}
|
|
|
|
private static T GenericRootN<T>(T value, int n) where T : IRootFunctions<T>
|
|
{
|
|
return T.RootN(value, n);
|
|
}
|
|
|
|
private static void AssertBound(DoubleDouble input, DoubleDouble actual, int degree)
|
|
{
|
|
DoubleDouble.IsFinite(actual).ShouldBeTrue();
|
|
DoubleDouble.IsCanonical(actual).ShouldBeTrue();
|
|
(actual.High > 0.0).ShouldBeTrue();
|
|
int n = Math.Abs(degree);
|
|
BigInteger x = Units(input.High) + Units(input.Low);
|
|
BigInteger y = Units(actual.High) + Units(actual.Low);
|
|
BigInteger scale = BigInteger.One << 100;
|
|
BigInteger power = BigInteger.Pow(y, n) << (100 * n);
|
|
BigInteger target;
|
|
if (degree > 0)
|
|
{
|
|
target = x << (1074 * (n - 1));
|
|
}
|
|
else
|
|
{
|
|
power *= x;
|
|
target = BigInteger.One << (1074 * (n + 1));
|
|
}
|
|
// Positive degree: x*(1-t)^n <= y^n <= x*(1+t)^n.
|
|
// Negative degree: (1-t)^n <= x*y^n <= (1+t)^n.
|
|
(power >= target * BigInteger.Pow(scale - 1, n)
|
|
&& power <= target * BigInteger.Pow(scale + 1, n)).ShouldBeTrue(
|
|
$"RootN bound failed for ({input.High:R}, {input.Low:R}), {degree}: ({actual.High:R}, {actual.Low:R})");
|
|
}
|
|
|
|
private static BigInteger Units(double value)
|
|
{
|
|
long bits = BitConverter.DoubleToInt64Bits(value);
|
|
int exponent = (int)((bits >> 52) & 0x7ff);
|
|
BigInteger significand = bits & 0xfffffffffffffL;
|
|
if (exponent != 0)
|
|
{
|
|
significand += BigInteger.One << 52;
|
|
significand <<= exponent - 1;
|
|
}
|
|
return bits < 0 ? -significand : significand;
|
|
}
|
|
|
|
private static void AssertBits(DoubleDouble expected, DoubleDouble actual)
|
|
{
|
|
BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(BitConverter.DoubleToInt64Bits(expected.High));
|
|
BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(BitConverter.DoubleToInt64Bits(expected.Low));
|
|
}
|
|
}
|