operators refactoring
.NET Test / .NET tests (push) Successful in 2m11s

This commit is contained in:
2026-09-24 00:39:59 +04:00
parent 3de9457fdc
commit ac93744579
26 changed files with 651 additions and 334 deletions
-46
View File
@@ -1,46 +0,0 @@
namespace Just.PreciseMath;
public static partial class DDMath
{
/// <inheritdoc cref="DoubleDouble.Exp"/>
[Pure]
public static DoubleDouble Exp(DoubleDouble value)
{
return DoubleDouble.Exp(value);
}
/// <inheritdoc cref="DoubleDouble.Exp2"/>
[Pure]
public static DoubleDouble Exp2(DoubleDouble value)
{
return DoubleDouble.Exp2(value);
}
/// <inheritdoc cref="DoubleDouble.Exp10"/>
[Pure]
public static DoubleDouble Exp10(DoubleDouble value)
{
return DoubleDouble.Exp10(value);
}
/// <inheritdoc cref="DoubleDouble.ExpM1"/>
[Pure]
public static DoubleDouble ExpM1(DoubleDouble value)
{
return DoubleDouble.ExpM1(value);
}
/// <inheritdoc cref="DoubleDouble.Exp2M1"/>
[Pure]
public static DoubleDouble Exp2M1(DoubleDouble value)
{
return DoubleDouble.Exp2M1(value);
}
/// <inheritdoc cref="DoubleDouble.Exp10M1"/>
[Pure]
public static DoubleDouble Exp10M1(DoubleDouble value)
{
return DoubleDouble.Exp10M1(value);
}
}
-70
View File
@@ -1,70 +0,0 @@
namespace Just.PreciseMath;
/// <summary>Provides mathematical functions for normalized double-double values.</summary>
/// <remarks>Results use fixed-size double-double precision, not arbitrary precision.</remarks>
public static partial class DDMath
{
/// <summary>Returns the absolute value without discarding the low component.</summary>
/// <remarks>Both signs of zero become positive zero; NaN is canonicalized and either infinity becomes positive infinity.</remarks>
[Pure]
public static DoubleDouble Abs(DoubleDouble value)
{
return DoubleDouble.Abs(value);
}
/// <inheritdoc cref="DoubleDouble.Square"/>
[Pure]
public static DoubleDouble Square(DoubleDouble value)
{
return DoubleDouble.Square(value);
}
/// <inheritdoc cref="DoubleDouble.Reciprocal"/>
[Pure]
public static DoubleDouble Reciprocal(DoubleDouble value)
{
return DoubleDouble.Reciprocal(value);
}
/// <inheritdoc cref="DoubleDouble.Cbrt"/>
[Pure]
public static DoubleDouble Cbrt(DoubleDouble value)
{
return DoubleDouble.Cbrt(value);
}
/// <inheritdoc cref="DoubleDouble.Hypot(DoubleDouble, DoubleDouble)"/>
[Pure]
public static DoubleDouble Hypot(DoubleDouble x, DoubleDouble y)
{
return DoubleDouble.Hypot(x, y);
}
/// <inheritdoc cref="DoubleDouble.Hypot(DoubleDouble, DoubleDouble, DoubleDouble)"/>
[Pure]
public static DoubleDouble Hypot(DoubleDouble x, DoubleDouble y, DoubleDouble z)
{
return DoubleDouble.Hypot(x, y, z);
}
/// <inheritdoc cref="DoubleDouble.RootN"/>
[Pure]
public static DoubleDouble RootN(DoubleDouble value, int n)
{
return DoubleDouble.RootN(value, n);
}
/// <inheritdoc cref="DoubleDouble.Sqrt"/>
[Pure]
public static DoubleDouble Sqrt(DoubleDouble value)
{
return DoubleDouble.Sqrt(value);
}
/// <inheritdoc cref="DoubleDouble.InvSqrt"/>
[Pure]
public static DoubleDouble InvSqrt(DoubleDouble value)
{
return DoubleDouble.InvSqrt(value);
}
}
@@ -44,13 +44,13 @@ public readonly partial struct DoubleDouble :
} }
/// <summary>Returns the operand unchanged.</summary> /// <summary>Returns the operand unchanged.</summary>
public static DoubleDouble operator +(DoubleDouble value) public static DoubleDouble Plus(DoubleDouble value)
{ {
return value; return value;
} }
/// <summary>Negates the value, including the high zero's sign, retaining canonical NaN and zero residuals.</summary> /// <summary>Negates the value, including the high zero's sign, retaining canonical NaN and zero residuals.</summary>
public static DoubleDouble operator -(DoubleDouble value) public static DoubleDouble Negate(DoubleDouble value)
{ {
// Negation preserves normalization; only NaN and zero residuals need canonicalization. // Negation preserves normalization; only NaN and zero residuals need canonicalization.
return new DoubleDouble(double.IsNaN(value._high) ? double.NaN : -value._high, return new DoubleDouble(double.IsNaN(value._high) ? double.NaN : -value._high,
@@ -58,7 +58,7 @@ public readonly partial struct DoubleDouble :
} }
/// <summary>Adds normalized expansions, retaining low-sum residuals under cancellation.</summary> /// <summary>Adds normalized expansions, retaining low-sum residuals under cancellation.</summary>
public static DoubleDouble operator +(DoubleDouble left, DoubleDouble right) public static DoubleDouble Add(DoubleDouble left, DoubleDouble right)
{ {
if (!IsFinite(left) || !IsFinite(right) || (left._high == 0.0 && right._high == 0.0)) if (!IsFinite(left) || !IsFinite(right) || (left._high == 0.0 && right._high == 0.0))
{ {
@@ -73,7 +73,7 @@ public readonly partial struct DoubleDouble :
} }
/// <summary>Subtracts normalized expansions.</summary> /// <summary>Subtracts normalized expansions.</summary>
public static DoubleDouble operator -(DoubleDouble left, DoubleDouble right) public static DoubleDouble Subtract(DoubleDouble left, DoubleDouble right)
{ {
if (!IsFinite(left) || !IsFinite(right) || (left._high == 0.0 && right._high == 0.0)) if (!IsFinite(left) || !IsFinite(right) || (left._high == 0.0 && right._high == 0.0))
{ {
@@ -95,7 +95,7 @@ public readonly partial struct DoubleDouble :
/// Swapping operands preserves both result components bit-for-bit. /// Swapping operands preserves both result components bit-for-bit.
/// Results are approximate double-double values, not universally correctly rounded. /// Results are approximate double-double values, not universally correctly rounded.
/// </remarks> /// </remarks>
public static DoubleDouble operator *(DoubleDouble left, DoubleDouble right) public static DoubleDouble Multiply(DoubleDouble left, DoubleDouble right)
{ {
if (!IsFinite(left) || !IsFinite(right) || left._high == 0.0 || right._high == 0.0) if (!IsFinite(left) || !IsFinite(right) || left._high == 0.0 || right._high == 0.0)
{ {
@@ -127,7 +127,7 @@ public readonly partial struct DoubleDouble :
/// <summary>Divides expansions using a quotient estimate and two residual corrections.</summary> /// <summary>Divides expansions using a quotient estimate and two residual corrections.</summary>
/// <remarks>Zero and nonfinite operands follow binary64 rules; precision decreases near underflow.</remarks> /// <remarks>Zero and nonfinite operands follow binary64 rules; precision decreases near underflow.</remarks>
public static DoubleDouble operator /(DoubleDouble left, DoubleDouble right) public static DoubleDouble Divide(DoubleDouble left, DoubleDouble right)
{ {
if (!IsFinite(left) || !IsFinite(right) || left._high == 0.0 || right._high == 0.0) if (!IsFinite(left) || !IsFinite(right) || left._high == 0.0 || right._high == 0.0)
{ {
@@ -159,35 +159,35 @@ public readonly partial struct DoubleDouble :
/// <summary>Applies the expansion operation without discarding the low component.</summary> /// <summary>Applies the expansion operation without discarding the low component.</summary>
/// <remarks>Preserves the component bits of addition with the equivalent zero-low scalar expansion.</remarks> /// <remarks>Preserves the component bits of addition with the equivalent zero-low scalar expansion.</remarks>
public static DoubleDouble operator +(DoubleDouble left, double right) public static DoubleDouble Add(DoubleDouble left, double right)
{ {
return PreciseMathHelper.AddScalar(left._high, left._low, right); return PreciseMathHelper.AddScalar(left._high, left._low, right);
} }
/// <summary>Applies the expansion operation without discarding the low component.</summary> /// <summary>Applies the expansion operation without discarding the low component.</summary>
/// <remarks>Preserves the component bits of addition with the equivalent zero-low scalar expansion.</remarks> /// <remarks>Preserves the component bits of addition with the equivalent zero-low scalar expansion.</remarks>
public static DoubleDouble operator +(double left, DoubleDouble right) public static DoubleDouble Add(double left, DoubleDouble right)
{ {
return right + left; return Add(right, left);
} }
/// <summary>Applies the expansion operation without discarding the low component.</summary> /// <summary>Applies the expansion operation without discarding the low component.</summary>
/// <remarks>Preserves the component bits of subtraction with the equivalent zero-low scalar expansion.</remarks> /// <remarks>Preserves the component bits of subtraction with the equivalent zero-low scalar expansion.</remarks>
public static DoubleDouble operator -(DoubleDouble left, double right) public static DoubleDouble Subtract(DoubleDouble left, double right)
{ {
return PreciseMathHelper.AddScalar(left._high, left._low, -right); return PreciseMathHelper.AddScalar(left._high, left._low, -right);
} }
/// <summary>Applies the expansion operation without discarding the low component.</summary> /// <summary>Applies the expansion operation without discarding the low component.</summary>
/// <remarks>Preserves the component bits of subtraction with the equivalent zero-low scalar expansion.</remarks> /// <remarks>Preserves the component bits of subtraction with the equivalent zero-low scalar expansion.</remarks>
public static DoubleDouble operator -(double left, DoubleDouble right) public static DoubleDouble Subtract(double left, DoubleDouble right)
{ {
// Negate the components, not the result: exact cancellation must yield +0. // Negate the components, not the result: exact cancellation must yield +0.
return PreciseMathHelper.AddScalar(-right._high, -right._low, left); return PreciseMathHelper.AddScalar(-right._high, -right._low, left);
} }
/// <summary>Applies the expansion operation without discarding the low component.</summary> /// <summary>Applies the expansion operation without discarding the low component.</summary>
public static DoubleDouble operator *(DoubleDouble left, double right) public static DoubleDouble Multiply(DoubleDouble left, double right)
{ {
if (!IsFinite(left) || !double.IsFinite(right) || left._high == 0.0 || right == 0.0) if (!IsFinite(left) || !double.IsFinite(right) || left._high == 0.0 || right == 0.0)
{ {
@@ -207,14 +207,14 @@ public readonly partial struct DoubleDouble :
} }
/// <summary>Applies the expansion operation without discarding the low component.</summary> /// <summary>Applies the expansion operation without discarding the low component.</summary>
public static DoubleDouble operator *(double left, DoubleDouble right) public static DoubleDouble Multiply(double left, DoubleDouble right)
{ {
return right * left; return Multiply(right, left);
} }
/// <summary>Divides by a scalar using a quotient estimate and two residual corrections.</summary> /// <summary>Divides by a scalar using a quotient estimate and two residual corrections.</summary>
/// <remarks>Preserves the component bits of division by the equivalent zero-low expansion.</remarks> /// <remarks>Preserves the component bits of division by the equivalent zero-low expansion.</remarks>
public static DoubleDouble operator /(DoubleDouble left, double right) public static DoubleDouble Divide(DoubleDouble left, double right)
{ {
if (!IsFinite(left) || !double.IsFinite(right) || left._high == 0.0 || right == 0.0) if (!IsFinite(left) || !double.IsFinite(right) || left._high == 0.0 || right == 0.0)
{ {
@@ -246,7 +246,7 @@ public readonly partial struct DoubleDouble :
/// <summary>Divides a scalar using a quotient estimate and two residual corrections.</summary> /// <summary>Divides a scalar using a quotient estimate and two residual corrections.</summary>
/// <remarks>Preserves the component bits of division with the equivalent zero-low numerator.</remarks> /// <remarks>Preserves the component bits of division with the equivalent zero-low numerator.</remarks>
public static DoubleDouble operator /(double left, DoubleDouble right) public static DoubleDouble Divide(double left, DoubleDouble right)
{ {
if (!double.IsFinite(left) || !IsFinite(right) || left == 0.0 || right._high == 0.0) if (!double.IsFinite(left) || !IsFinite(right) || left == 0.0 || right._high == 0.0)
{ {
@@ -1,6 +1,6 @@
namespace Just.PreciseMath; namespace Just.PreciseMath;
public static partial class DDMath public readonly partial struct DoubleDouble
{ {
/// <summary>Returns the natural logarithm of a double-double value.</summary> /// <summary>Returns the natural logarithm of a double-double value.</summary>
/// <remarks> /// <remarks>
@@ -17,7 +17,7 @@ public static partial class DDMath
{ {
if (double.IsNaN(value.High) || value.High < 0.0) if (double.IsNaN(value.High) || value.High < 0.0)
{ {
return DoubleDouble.NaN; return NaN;
} }
if (value.High == 0.0) if (value.High == 0.0)
{ {
@@ -27,9 +27,9 @@ public static partial class DDMath
{ {
return value; return value;
} }
if (value == DoubleDouble.One) if (value == One)
{ {
return DoubleDouble.Zero; return Zero;
} }
return LogPositiveFinite(value); return LogPositiveFinite(value);
} }
@@ -49,9 +49,9 @@ public static partial class DDMath
} }
int exponent = Math.ILogB(value.High); int exponent = Math.ILogB(value.High);
DoubleDouble mantissa = DoubleDouble.FromComponents(Math.ScaleB(value.High, -exponent), DoubleDouble mantissa = FromComponents(Math.ScaleB(value.High, -exponent),
Math.ScaleB(value.Low, -exponent)); Math.ScaleB(value.Low, -exponent));
if (mantissa.High > DoubleDouble.Sqrt2.High) if (mantissa.High > Sqrt2.High)
{ {
mantissa *= 0.5; mantissa *= 0.5;
++exponent; ++exponent;
@@ -72,6 +72,6 @@ public static partial class DDMath
} }
// log(m) = 2*atanh(t). The omitted series tail is bounded by // log(m) = 2*atanh(t). The omitted series tail is bounded by
// 2*|t|^51/(51*(1-t^2)); rounding, not truncation, dominates. // 2*|t|^51/(51*(1-t^2)); rounding, not truncation, dominates.
return (sum * 2.0) + (DoubleDouble.Ln2 * exponent); return (sum * 2.0) + (Ln2 * exponent);
} }
} }
@@ -86,18 +86,6 @@ public readonly partial struct DoubleDouble
return IsNaN(y) ? x : MinMagnitude(x, y); return IsNaN(y) ? x : MinMagnitude(x, y);
} }
/// <summary>Adds one using double-double arithmetic, including its overflow and nonfinite behavior.</summary>
public static DoubleDouble operator ++(DoubleDouble value)
{
return value + 1.0;
}
/// <summary>Subtracts one using double-double arithmetic, including its overflow and nonfinite behavior.</summary>
public static DoubleDouble operator --(DoubleDouble value)
{
return value - 1.0;
}
/// <summary>Tests normalization and canonical NaN and zero-residual bit patterns.</summary> /// <summary>Tests normalization and canonical NaN and zero-residual bit patterns.</summary>
public static bool IsCanonical(DoubleDouble value) public static bool IsCanonical(DoubleDouble value)
{ {
@@ -0,0 +1,100 @@
namespace Just.PreciseMath;
public readonly partial struct DoubleDouble
{
/// <inheritdoc cref="Plus"/>
public static DoubleDouble operator +(DoubleDouble value)
{
return Plus(value);
}
/// <inheritdoc cref="Negate"/>
public static DoubleDouble operator -(DoubleDouble value)
{
return Negate(value);
}
/// <inheritdoc cref="Add(DoubleDouble, DoubleDouble)"/>
public static DoubleDouble operator +(DoubleDouble left, DoubleDouble right)
{
return Add(left, right);
}
/// <inheritdoc cref="Subtract(DoubleDouble, DoubleDouble)"/>
public static DoubleDouble operator -(DoubleDouble left, DoubleDouble right)
{
return Subtract(left, right);
}
/// <inheritdoc cref="Multiply(DoubleDouble, DoubleDouble)"/>
public static DoubleDouble operator *(DoubleDouble left, DoubleDouble right)
{
return Multiply(left, right);
}
/// <inheritdoc cref="Divide(DoubleDouble, DoubleDouble)"/>
public static DoubleDouble operator /(DoubleDouble left, DoubleDouble right)
{
return Divide(left, right);
}
/// <inheritdoc cref="Add(DoubleDouble, double)"/>
public static DoubleDouble operator +(DoubleDouble left, double right)
{
return Add(left, right);
}
/// <inheritdoc cref="Add(double, DoubleDouble)"/>
public static DoubleDouble operator +(double left, DoubleDouble right)
{
return Add(left, right);
}
/// <inheritdoc cref="Subtract(DoubleDouble, double)"/>
public static DoubleDouble operator -(DoubleDouble left, double right)
{
return Subtract(left, right);
}
/// <inheritdoc cref="Subtract(double, DoubleDouble)"/>
public static DoubleDouble operator -(double left, DoubleDouble right)
{
return Subtract(left, right);
}
/// <inheritdoc cref="Multiply(DoubleDouble, double)"/>
public static DoubleDouble operator *(DoubleDouble left, double right)
{
return Multiply(left, right);
}
/// <inheritdoc cref="Multiply(double, DoubleDouble)"/>
public static DoubleDouble operator *(double left, DoubleDouble right)
{
return Multiply(left, right);
}
/// <inheritdoc cref="Divide(DoubleDouble, double)"/>
public static DoubleDouble operator /(DoubleDouble left, double right)
{
return Divide(left, right);
}
/// <inheritdoc cref="Divide(double, DoubleDouble)"/>
public static DoubleDouble operator /(double left, DoubleDouble right)
{
return Divide(left, right);
}
/// <summary>Adds one using double-double arithmetic, including its overflow and nonfinite behavior.</summary>
public static DoubleDouble operator ++(DoubleDouble value)
{
return Add(value, 1.0);
}
/// <summary>Subtracts one using double-double arithmetic, including its overflow and nonfinite behavior.</summary>
public static DoubleDouble operator --(DoubleDouble value)
{
return Subtract(value, 1.0);
}
}
@@ -1,6 +1,6 @@
namespace Just.PreciseMath; namespace Just.PreciseMath;
public static partial class DDMath public readonly partial struct DoubleDouble
{ {
/// <summary>Raises a double-double value to an integer power.</summary> /// <summary>Raises a double-double value to an integer power.</summary>
/// <remarks> /// <remarks>
@@ -21,11 +21,11 @@ public static partial class DDMath
{ {
if (exponent == 0) if (exponent == 0)
{ {
return DoubleDouble.One; return One;
} }
if (double.IsNaN(value.High)) if (double.IsNaN(value.High))
{ {
return DoubleDouble.NaN; return NaN;
} }
if (exponent == 1) if (exponent == 1)
{ {
@@ -43,10 +43,10 @@ public static partial class DDMath
// Widen BEFORE negation: abs(int.MinValue) is not an int. // Widen BEFORE negation: abs(int.MinValue) is not an int.
long remaining = exponent < 0 ? -(long)exponent : exponent; long remaining = exponent < 0 ? -(long)exponent : exponent;
long factorExponent = Math.ILogB(value.High); long factorExponent = Math.ILogB(value.High);
DoubleDouble factor = DoubleDouble.FromComponents( DoubleDouble factor = FromComponents(
Math.ScaleB(Math.Abs(value.High), -(int)factorExponent), Math.ScaleB(Math.Abs(value.High), -(int)factorExponent),
Math.ScaleB(value.High < 0.0 ? -value.Low : value.Low, -(int)factorExponent)); Math.ScaleB(value.High < 0.0 ? -value.Low : value.Low, -(int)factorExponent));
DoubleDouble result = DoubleDouble.One; DoubleDouble result = One;
long resultExponent = 0; long resultExponent = 0;
while (remaining != 0) while (remaining != 0)
{ {
@@ -105,35 +105,35 @@ public static partial class DDMath
[Pure] [Pure]
public static DoubleDouble Pow(DoubleDouble value, DoubleDouble exponent) public static DoubleDouble Pow(DoubleDouble value, DoubleDouble exponent)
{ {
if (DoubleDouble.IsZero(exponent) || value == DoubleDouble.One) if (IsZero(exponent) || value == One)
{ {
return DoubleDouble.One; return One;
} }
if (DoubleDouble.IsNaN(value) || DoubleDouble.IsNaN(exponent)) if (IsNaN(value) || IsNaN(exponent))
{ {
return DoubleDouble.NaN; return NaN;
} }
DoubleDouble magnitude = DoubleDouble.Abs(value); DoubleDouble magnitude = Abs(value);
if (DoubleDouble.IsInfinity(exponent)) if (IsInfinity(exponent))
{ {
if (magnitude == DoubleDouble.One) if (magnitude == One)
{ {
return DoubleDouble.One; return One;
} }
bool grows = (magnitude > DoubleDouble.One) == (exponent.High > 0.0); bool grows = (magnitude > One) == (exponent.High > 0.0);
return grows ? new DoubleDouble(double.PositiveInfinity) : DoubleDouble.Zero; return grows ? new DoubleDouble(double.PositiveInfinity) : Zero;
} }
bool integer = DoubleDouble.IsInteger(exponent); bool integer = IsInteger(exponent);
if (integer && exponent.High >= int.MinValue && exponent.High <= int.MaxValue) if (integer && exponent.High >= int.MinValue && exponent.High <= int.MaxValue)
{ {
// A normalized integer in this range has a zero low component. // A normalized integer in this range has a zero low component.
// Reuse the integer algorithm, including exact identities and ties. // Reuse the integer algorithm, including exact identities and ties.
return Pow(value, (int)exponent.High); return Pow(value, (int)exponent.High);
} }
bool negative = double.IsNegative(value.High) && DoubleDouble.IsOddInteger(exponent); bool negative = double.IsNegative(value.High) && IsOddInteger(exponent);
if (value.High == 0.0 || DoubleDouble.IsInfinity(value)) if (value.High == 0.0 || IsInfinity(value))
{ {
bool infinite = (value.High == 0.0) == (exponent.High < 0.0); bool infinite = (value.High == 0.0) == (exponent.High < 0.0);
double result = infinite ? double.PositiveInfinity : 0.0; double result = infinite ? double.PositiveInfinity : 0.0;
@@ -141,11 +141,11 @@ public static partial class DDMath
} }
if (value.High < 0.0 && !integer) if (value.High < 0.0 && !integer)
{ {
return DoubleDouble.NaN; return NaN;
} }
if (magnitude == DoubleDouble.One) if (magnitude == One)
{ {
return negative ? DoubleDouble.NegativeOne : DoubleDouble.One; return negative ? NegativeOne : One;
} }
DoubleDouble resultMagnitude = Exp(exponent * LogPositiveFinite(magnitude)); DoubleDouble resultMagnitude = Exp(exponent * LogPositiveFinite(magnitude));
return negative ? -resultMagnitude : resultMagnitude; return negative ? -resultMagnitude : resultMagnitude;
@@ -158,7 +158,7 @@ public static partial class DDMath
// int.MinValue times any finite input's binary64 exponent. // int.MinValue times any finite input's binary64 exponent.
int adjustment = Math.ILogB(value.High); int adjustment = Math.ILogB(value.High);
exponent += adjustment; exponent += adjustment;
return DoubleDouble.FromComponents(Math.ScaleB(value.High, -adjustment), return FromComponents(Math.ScaleB(value.High, -adjustment),
Math.ScaleB(value.Low, -adjustment)); Math.ScaleB(value.Low, -adjustment));
} }
} }
@@ -200,7 +200,7 @@ public readonly partial struct DoubleDouble : IRootFunctions<DoubleDouble>
case 1: case 1:
return value; return value;
case -1: case -1:
return DDMath.Reciprocal(value); return Reciprocal(value);
case 2: case 2:
return Sqrt(value); return Sqrt(value);
case -2: case -2:
@@ -0,0 +1,128 @@
namespace Just.PreciseMath;
/// <summary>Provides instance-style calls to the static DoubleDouble mathematical functions.</summary>
/// <remarks>
/// Every extension returns a new result without modifying its receiver. Extensions
/// forward directly to DoubleDouble, retaining its numerical and special-value contracts.
/// </remarks>
public static class DoubleDoubleExtensions
{
/// <inheritdoc cref="DoubleDouble.Abs(DoubleDouble)"/>
[Pure]
public static DoubleDouble Abs(this DoubleDouble value)
{
return DoubleDouble.Abs(value);
}
/// <inheritdoc cref="DoubleDouble.Square(DoubleDouble)"/>
[Pure]
public static DoubleDouble Square(this DoubleDouble value)
{
return DoubleDouble.Square(value);
}
/// <inheritdoc cref="DoubleDouble.Reciprocal(DoubleDouble)"/>
[Pure]
public static DoubleDouble Reciprocal(this DoubleDouble value)
{
return DoubleDouble.Reciprocal(value);
}
/// <inheritdoc cref="DoubleDouble.Sqrt(DoubleDouble)"/>
[Pure]
public static DoubleDouble Sqrt(this DoubleDouble value)
{
return DoubleDouble.Sqrt(value);
}
/// <inheritdoc cref="DoubleDouble.InvSqrt(DoubleDouble)"/>
[Pure]
public static DoubleDouble InvSqrt(this DoubleDouble value)
{
return DoubleDouble.InvSqrt(value);
}
/// <inheritdoc cref="DoubleDouble.Cbrt(DoubleDouble)"/>
[Pure]
public static DoubleDouble Cbrt(this DoubleDouble value)
{
return DoubleDouble.Cbrt(value);
}
/// <inheritdoc cref="DoubleDouble.Exp(DoubleDouble)"/>
[Pure]
public static DoubleDouble Exp(this DoubleDouble value)
{
return DoubleDouble.Exp(value);
}
/// <inheritdoc cref="DoubleDouble.Exp2(DoubleDouble)"/>
[Pure]
public static DoubleDouble Exp2(this DoubleDouble value)
{
return DoubleDouble.Exp2(value);
}
/// <inheritdoc cref="DoubleDouble.Exp10(DoubleDouble)"/>
[Pure]
public static DoubleDouble Exp10(this DoubleDouble value)
{
return DoubleDouble.Exp10(value);
}
/// <inheritdoc cref="DoubleDouble.ExpM1(DoubleDouble)"/>
[Pure]
public static DoubleDouble ExpM1(this DoubleDouble value)
{
return DoubleDouble.ExpM1(value);
}
/// <inheritdoc cref="DoubleDouble.Exp2M1(DoubleDouble)"/>
[Pure]
public static DoubleDouble Exp2M1(this DoubleDouble value)
{
return DoubleDouble.Exp2M1(value);
}
/// <inheritdoc cref="DoubleDouble.Exp10M1(DoubleDouble)"/>
[Pure]
public static DoubleDouble Exp10M1(this DoubleDouble value)
{
return DoubleDouble.Exp10M1(value);
}
/// <inheritdoc cref="DoubleDouble.Log(DoubleDouble)"/>
[Pure]
public static DoubleDouble Log(this DoubleDouble value)
{
return DoubleDouble.Log(value);
}
/// <inheritdoc cref="DoubleDouble.RootN(DoubleDouble, int)"/>
[Pure]
public static DoubleDouble RootN(this DoubleDouble value, int n)
{
return DoubleDouble.RootN(value, n);
}
/// <inheritdoc cref="DoubleDouble.Pow(DoubleDouble, int)"/>
[Pure]
public static DoubleDouble Pow(this DoubleDouble value, int exponent)
{
return DoubleDouble.Pow(value, exponent);
}
/// <inheritdoc cref="DoubleDouble.Pow(DoubleDouble, double)"/>
[Pure]
public static DoubleDouble Pow(this DoubleDouble value, double exponent)
{
return DoubleDouble.Pow(value, exponent);
}
/// <inheritdoc cref="DoubleDouble.Pow(DoubleDouble, DoubleDouble)"/>
[Pure]
public static DoubleDouble Pow(this DoubleDouble value, DoubleDouble exponent)
{
return DoubleDouble.Pow(value, exponent);
}
}
@@ -2,6 +2,55 @@ namespace Just.PreciseMath.Tests;
public class DoubleDoubleArithmeticTests public class DoubleDoubleArithmeticTests
{ {
[Theory]
[InlineData("+")]
[InlineData("-")]
[InlineData("*")]
[InlineData("/")]
public void NamedArithmeticMatchesOperatorComponentBits(string operation)
{
// Compatibility, not an accuracy oracle. The existing rational tests
// independently validate the kernels reached through these operators.
DoubleDouble[] values = MultiplicationCommutativityCases();
foreach (DoubleDouble left in values)
{
foreach (DoubleDouble right in values)
{
(DoubleDouble actual, DoubleDouble expected)[] pairs = operation switch
{
"+" => [(DoubleDouble.Add(left, right), left + right),
(DoubleDouble.Add(left, right.High), left + right.High),
(DoubleDouble.Add(left.High, right), left.High + right)],
"-" => [(DoubleDouble.Subtract(left, right), left - right),
(DoubleDouble.Subtract(left, right.High), left - right.High),
(DoubleDouble.Subtract(left.High, right), left.High - right)],
"*" => [(DoubleDouble.Multiply(left, right), left * right),
(DoubleDouble.Multiply(left, right.High), left * right.High),
(DoubleDouble.Multiply(left.High, right), left.High * right)],
"/" => [(DoubleDouble.Divide(left, right), left / right),
(DoubleDouble.Divide(left, right.High), left / right.High),
(DoubleDouble.Divide(left.High, right), left.High / right)],
_ => throw new ArgumentOutOfRangeException(nameof(operation))
};
foreach ((DoubleDouble actual, DoubleDouble expected) in pairs)
{
CheckBoundary(actual, expected.High, expected.Low);
}
}
}
}
[Fact]
public void NamedUnaryArithmeticMatchesOperatorComponentBits()
{
foreach (DoubleDouble value in MultiplicationCommutativityCases())
{
CheckBoundary(DoubleDouble.Plus(value), value.High, value.Low);
DoubleDouble negated = -value;
CheckBoundary(DoubleDouble.Negate(value), negated.High, negated.Low);
}
}
[Theory] [Theory]
[InlineData(6.0, 2.0, 8.0, 4.0, 12.0, 3.0)] [InlineData(6.0, 2.0, 8.0, 4.0, 12.0, 3.0)]
[InlineData(-6.0, 2.0, -4.0, -8.0, -12.0, -3.0)] [InlineData(-6.0, 2.0, -4.0, -8.0, -12.0, -3.0)]
@@ -1050,7 +1099,7 @@ public class DoubleDoubleArithmeticTests
foreach (double sign in new[] { -1.0, 1.0 }) foreach (double sign in new[] { -1.0, 1.0 })
{ {
CheckBoundary(DoubleDouble.Square(sign * below), Math.BitDecrement(highSquared), Math.ScaleB(1.0, (2 * exponent) - 108)); CheckBoundary(DoubleDouble.Square(sign * below), Math.BitDecrement(highSquared), Math.ScaleB(1.0, (2 * exponent) - 108));
CheckBoundary(DDMath.Square(sign * above), Math.BitIncrement(highSquared), Math.ScaleB(1.0, (2 * exponent) - 106)); CheckBoundary(DoubleDouble.Square(sign * above), Math.BitIncrement(highSquared), Math.ScaleB(1.0, (2 * exponent) - 106));
} }
} }
@@ -1060,7 +1109,6 @@ public class DoubleDoubleArithmeticTests
foreach (double value in new[] { 0.0, -0.0, double.NaN, double.PositiveInfinity, double.NegativeInfinity }) foreach (double value in new[] { 0.0, -0.0, double.NaN, double.PositiveInfinity, double.NegativeInfinity })
{ {
CheckBits(DoubleDouble.Square(new DoubleDouble(value)), value * value); CheckBits(DoubleDouble.Square(new DoubleDouble(value)), value * value);
CheckBits(DDMath.Square(new DoubleDouble(value)), value * value);
} }
for (int exponent = -1074; exponent <= 1023; exponent++) for (int exponent = -1074; exponent <= 1023; exponent++)
{ {
@@ -1139,7 +1187,6 @@ public class DoubleDoubleArithmeticTests
DoubleDouble.IsZero(actual).ShouldBe((numerator << 1) <= denominator); DoubleDouble.IsZero(actual).ShouldBe((numerator << 1) <= denominator);
} }
CheckBoundary(DoubleDouble.Square(-value), actual.High, actual.Low); CheckBoundary(DoubleDouble.Square(-value), actual.High, actual.Low);
CheckBoundary(DDMath.Square(value), actual.High, actual.Low);
} }
private static DoubleDouble[] MultiplicationCommutativityCases() private static DoubleDouble[] MultiplicationCommutativityCases()
@@ -36,7 +36,6 @@ public class DoubleDoubleCbrtTests
{ {
AssertBound(input, DoubleDouble.Cbrt(input)); AssertBound(input, DoubleDouble.Cbrt(input));
AssertBound(-input, DoubleDouble.Cbrt(-input)); AssertBound(-input, DoubleDouble.Cbrt(-input));
AssertBits(DoubleDouble.Cbrt(input), DDMath.Cbrt(input));
} }
} }
} }
@@ -20,7 +20,6 @@ public class DoubleDoubleExponentialFunctionsTests
DoubleDouble actual = Evaluate(operation, input); DoubleDouble actual = Evaluate(operation, input);
DoubleDouble.IsCanonical(actual).ShouldBeTrue(); DoubleDouble.IsCanonical(actual).ShouldBeTrue();
AssertBits(actual, EvaluateGeneric(operation, input)); AssertBits(actual, EvaluateGeneric(operation, input));
AssertBits(actual, EvaluateFacade(operation, input));
if (overflow) if (overflow)
{ {
AssertBits(DoubleDouble.PositiveInfinity, actual); AssertBits(DoubleDouble.PositiveInfinity, actual);
@@ -57,7 +56,7 @@ public class DoubleDoubleExponentialFunctionsTests
} }
[Fact] [Fact]
public void NaturalExponentialIsAvailableOnTheTypeAndMatchesFacadeBits() public void NaturalExponentialMatchesGenericDispatchBits()
{ {
DoubleDouble[] values = [DoubleDouble.Zero, DoubleDouble.NegativeZero, DoubleDouble[] values = [DoubleDouble.Zero, DoubleDouble.NegativeZero,
DoubleDouble.NaN, DoubleDouble.PositiveInfinity, DoubleDouble.NegativeInfinity, DoubleDouble.NaN, DoubleDouble.PositiveInfinity, DoubleDouble.NegativeInfinity,
@@ -65,7 +64,7 @@ public class DoubleDoubleExponentialFunctionsTests
new(double.MaxValue), DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -1000))]; new(double.MaxValue), DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -1000))];
foreach (DoubleDouble value in values) foreach (DoubleDouble value in values)
{ {
AssertBits(DDMath.Exp(value), DoubleDouble.Exp(value)); AssertBits(DoubleDouble.Exp(value), EvaluateGeneric("Exp", value));
} }
} }
@@ -117,7 +116,6 @@ public class DoubleDoubleExponentialFunctionsTests
DoubleDouble input = DoubleDouble.FromComponents(exponent, sign * double.Epsilon); DoubleDouble input = DoubleDouble.FromComponents(exponent, sign * double.Epsilon);
DoubleDouble expected = DoubleDouble.FromComponents(Math.ScaleB(1.0, exponent), sign * expectedLow); DoubleDouble expected = DoubleDouble.FromComponents(Math.ScaleB(1.0, exponent), sign * expectedLow);
AssertBits(expected, DoubleDouble.Exp2(input)); AssertBits(expected, DoubleDouble.Exp2(input));
AssertBits(expected, DDMath.Exp2(input));
} }
} }
@@ -128,7 +126,7 @@ public class DoubleDoubleExponentialFunctionsTests
[InlineData("ExpM1")] [InlineData("ExpM1")]
[InlineData("Exp2M1")] [InlineData("Exp2M1")]
[InlineData("Exp10M1")] [InlineData("Exp10M1")]
public void CompleteContractHasGenericDispatchAndThinFacade(string operation) public void CompleteContractHasGenericDispatch(string operation)
{ {
bool minusOne = operation.EndsWith("M1", StringComparison.Ordinal); bool minusOne = operation.EndsWith("M1", StringComparison.Ordinal);
DoubleDouble[] values = [DoubleDouble.Zero, DoubleDouble.NegativeZero, DoubleDouble[] values = [DoubleDouble.Zero, DoubleDouble.NegativeZero,
@@ -141,7 +139,6 @@ public class DoubleDoubleExponentialFunctionsTests
DoubleDouble actual = Evaluate(operation, value); DoubleDouble actual = Evaluate(operation, value);
DoubleDouble.IsCanonical(actual).ShouldBeTrue(); DoubleDouble.IsCanonical(actual).ShouldBeTrue();
AssertBits(actual, EvaluateGeneric(operation, value)); AssertBits(actual, EvaluateGeneric(operation, value));
AssertBits(actual, EvaluateFacade(operation, value));
if (DoubleDouble.IsZero(value)) if (DoubleDouble.IsZero(value))
{ {
AssertBits(minusOne ? value : DoubleDouble.One, actual); AssertBits(minusOne ? value : DoubleDouble.One, actual);
@@ -189,20 +186,6 @@ public class DoubleDoubleExponentialFunctionsTests
}; };
} }
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) private static BigInteger Units(double value)
{ {
long bits = BitConverter.DoubleToInt64Bits(value); long bits = BitConverter.DoubleToInt64Bits(value);
@@ -0,0 +1,165 @@
namespace Just.PreciseMath.Tests;
public class DoubleDoubleExtensionsTests
{
[Theory]
[InlineData("Abs")]
[InlineData("Square")]
[InlineData("Reciprocal")]
[InlineData("Sqrt")]
[InlineData("InvSqrt")]
[InlineData("Cbrt")]
[InlineData("Exp")]
[InlineData("Exp2")]
[InlineData("Exp10")]
[InlineData("ExpM1")]
[InlineData("Exp2M1")]
[InlineData("Exp10M1")]
[InlineData("Log")]
public void UnaryExtensionsPreserveStaticFunctionComponentBits(string operation)
{
// These are forwarding contracts, not accuracy oracles. The static
// kernels have separate independent numerical-reference tests.
foreach (DoubleDouble value in Values())
{
(DoubleDouble expected, DoubleDouble actual) = operation switch
{
"Abs" => (DoubleDouble.Abs(value), value.Abs()),
"Square" => (DoubleDouble.Square(value), value.Square()),
"Reciprocal" => (DoubleDouble.Reciprocal(value), value.Reciprocal()),
"Sqrt" => (DoubleDouble.Sqrt(value), value.Sqrt()),
"InvSqrt" => (DoubleDouble.InvSqrt(value), value.InvSqrt()),
"Cbrt" => (DoubleDouble.Cbrt(value), value.Cbrt()),
"Exp" => (DoubleDouble.Exp(value), value.Exp()),
"Exp2" => (DoubleDouble.Exp2(value), value.Exp2()),
"Exp10" => (DoubleDouble.Exp10(value), value.Exp10()),
"ExpM1" => (DoubleDouble.ExpM1(value), value.ExpM1()),
"Exp2M1" => (DoubleDouble.Exp2M1(value), value.Exp2M1()),
"Exp10M1" => (DoubleDouble.Exp10M1(value), value.Exp10M1()),
"Log" => (DoubleDouble.Log(value), value.Log()),
_ => throw new ArgumentOutOfRangeException(nameof(operation))
};
AssertBits(expected, actual);
}
}
[Fact]
public void ExtensionsCanBeChainedWithoutChangingTheReceiver()
{
DoubleDouble value = new(-4.0);
DoubleDouble actual = value.Abs().Square().Sqrt().Reciprocal();
// All intermediate results are exactly representable binary powers.
AssertBits(new DoubleDouble(0.25), actual);
AssertBits(new DoubleDouble(-4.0), value);
}
[Theory]
[InlineData(int.MinValue)]
[InlineData(-7)]
[InlineData(-2)]
[InlineData(-1)]
[InlineData(0)]
[InlineData(1)]
[InlineData(2)]
[InlineData(7)]
[InlineData(int.MaxValue)]
public void IntegerPowerExtensionPreservesStaticFunctionComponentBits(int exponent)
{
foreach (DoubleDouble value in Values())
{
AssertBits(DoubleDouble.Pow(value, exponent), value.Pow(exponent));
}
}
[Fact]
public void ScalarPowerExtensionPreservesStaticFunctionComponentBits()
{
double[] exponents = [0.0, -0.0, 1.0, -1.0, 2.0, -2.0, 0.5, -0.5,
1.25, -1.25, double.Epsilon, double.MaxValue, double.MinValue,
double.NaN, double.PositiveInfinity, double.NegativeInfinity];
foreach (DoubleDouble value in Values())
{
foreach (double exponent in exponents)
{
AssertBits(DoubleDouble.Pow(value, exponent), value.Pow(exponent));
}
}
}
[Fact]
public void ExpansionPowerExtensionPreservesBothExponentComponents()
{
DoubleDouble[] exponents = [DoubleDouble.Zero, DoubleDouble.NegativeZero,
DoubleDouble.One, DoubleDouble.NegativeOne, DoubleDouble.NaN,
DoubleDouble.PositiveInfinity, DoubleDouble.NegativeInfinity,
DoubleDouble.FromComponents(0.5, Math.ScaleB(1.0, -54)),
DoubleDouble.FromComponents(0.5, -Math.ScaleB(1.0, -55)),
DoubleDouble.FromComponents(1.0, double.Epsilon),
DoubleDouble.FromComponents(Math.ScaleB(1.0, 54), 1.0)];
foreach (DoubleDouble value in Values())
{
foreach (DoubleDouble exponent in exponents)
{
AssertBits(DoubleDouble.Pow(value, exponent), value.Pow(exponent));
}
}
}
[Theory]
[InlineData(int.MinValue)]
[InlineData(-7)]
[InlineData(-3)]
[InlineData(-2)]
[InlineData(-1)]
[InlineData(0)]
[InlineData(1)]
[InlineData(2)]
[InlineData(3)]
[InlineData(7)]
[InlineData(int.MaxValue)]
public void RootExtensionPreservesStaticFunctionComponentBits(int degree)
{
foreach (DoubleDouble value in Values())
{
AssertBits(DoubleDouble.RootN(value, degree), value.RootN(degree));
}
}
private static IEnumerable<DoubleDouble> Values()
{
yield return DoubleDouble.Zero;
yield return DoubleDouble.NegativeZero;
yield return DoubleDouble.NaN;
yield return DoubleDouble.PositiveInfinity;
yield return DoubleDouble.NegativeInfinity;
yield return new DoubleDouble(double.MaxValue);
yield return new DoubleDouble(double.MinValue);
yield return DoubleDouble.One;
yield return DoubleDouble.NegativeOne;
yield return DoubleDouble.Pi;
yield return -DoubleDouble.Pi;
// Square's dedicated kernel can differ from multiplication; near-one
// and sparse pairs also exercise low-component retention in roots/log.
yield return DoubleDouble.FromComponents(1.0, -Math.ScaleB(1.0, -54));
yield return DoubleDouble.FromComponents(1.0, double.Epsilon);
yield return DoubleDouble.FromComponents(1.0, -double.Epsilon);
foreach (int exponent in new[] { -1074, -1022, -451, -450, -1, 0, 1, 450, 451, 1020, 1021, 1023 })
{
foreach (double sign in new[] { -1.0, 1.0 })
{
double high = Math.ScaleB(sign * 1.25, exponent);
double denseLow = Math.ScaleB(1.0, exponent - 54);
foreach (double low in new[] { 0.0, denseLow, -denseLow, double.Epsilon, -double.Epsilon })
{
yield return DoubleDouble.FromComponents(high, low);
}
}
}
}
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));
}
}
@@ -52,7 +52,6 @@ public class DoubleDoubleHypotTests
DoubleDouble result = DoubleDouble.Hypot(x, y); DoubleDouble result = DoubleDouble.Hypot(x, y);
AssertBound(x, y, result); AssertBound(x, y, result);
AssertBits(result, DoubleDouble.Hypot(-y, x)); AssertBits(result, DoubleDouble.Hypot(-y, x));
AssertBits(result, DDMath.Hypot(x, y));
} }
} }
} }
@@ -134,7 +133,6 @@ public class DoubleDoubleHypotTests
{ {
AssertBound(new DoubleDouble(x), new DoubleDouble(y), new DoubleDouble(z), result); AssertBound(new DoubleDouble(x), new DoubleDouble(y), new DoubleDouble(z), result);
} }
AssertBits(result, DDMath.Hypot(new DoubleDouble(x), new DoubleDouble(y), new DoubleDouble(z)));
} }
} }
} }
@@ -186,7 +184,6 @@ public class DoubleDoubleHypotTests
AssertBound(x, y, z, result); AssertBound(x, y, z, result);
AssertBits(result, DoubleDouble.Hypot(-z, x, -y)); AssertBits(result, DoubleDouble.Hypot(-z, x, -y));
AssertBits(result, DoubleDouble.Hypot(y, -z, x)); AssertBits(result, DoubleDouble.Hypot(y, -z, x));
AssertBits(result, DDMath.Hypot(x, y, z));
} }
} }
} }
@@ -271,7 +268,6 @@ public class DoubleDoubleHypotTests
AssertBits(expected, DoubleDouble.Hypot(b, c, a)); AssertBits(expected, DoubleDouble.Hypot(b, c, a));
AssertBits(expected, DoubleDouble.Hypot(c, a, b)); AssertBits(expected, DoubleDouble.Hypot(c, a, b));
AssertBits(expected, DoubleDouble.Hypot(c, b, a)); AssertBits(expected, DoubleDouble.Hypot(c, b, a));
AssertBits(expected, DDMath.Hypot(a, b, c));
} }
} }
} }
@@ -294,7 +290,6 @@ public class DoubleDoubleHypotTests
AssertBound(first, second, result); AssertBound(first, second, result);
} }
AssertBits(result, DoubleDouble.Hypot(second, first)); AssertBits(result, DoubleDouble.Hypot(second, first));
AssertBits(result, DDMath.Hypot(first, second));
} }
} }
} }
@@ -3,7 +3,7 @@ namespace Just.PreciseMath.Tests;
public class DoubleDoubleRootFunctionsTests public class DoubleDoubleRootFunctionsTests
{ {
[Fact] [Fact]
public void ExistingRootsAreAvailableOnTheTypeAndMatchFacadeBits() public void RootsAreCanonicalAndSqrtMatchesGenericDispatchBits()
{ {
DoubleDouble[] values = [DoubleDouble.Zero, DoubleDouble.NegativeZero, DoubleDouble[] values = [DoubleDouble.Zero, DoubleDouble.NegativeZero,
DoubleDouble.NaN, DoubleDouble.PositiveInfinity, DoubleDouble.NegativeInfinity, DoubleDouble.NaN, DoubleDouble.PositiveInfinity, DoubleDouble.NegativeInfinity,
@@ -11,11 +11,16 @@ public class DoubleDoubleRootFunctionsTests
DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -1000))]; DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -1000))];
foreach (DoubleDouble value in values) foreach (DoubleDouble value in values)
{ {
AssertBits(DDMath.Sqrt(value), DoubleDouble.Sqrt(value)); AssertBits(DoubleDouble.Sqrt(value), GenericSqrt(value));
AssertBits(DDMath.InvSqrt(value), DoubleDouble.InvSqrt(value)); DoubleDouble.IsCanonical(DoubleDouble.InvSqrt(value)).ShouldBeTrue();
} }
} }
private static T GenericSqrt<T>(T value) where T : IRootFunctions<T>
{
return T.Sqrt(value);
}
private static void AssertBits(DoubleDouble expected, DoubleDouble actual) private static void AssertBits(DoubleDouble expected, DoubleDouble actual)
{ {
BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(BitConverter.DoubleToInt64Bits(expected.High)); BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(BitConverter.DoubleToInt64Bits(expected.High));
@@ -68,7 +68,7 @@ public class DoubleDoubleRootNTests
{ {
AssertBits(value, DoubleDouble.RootN(value, 1)); AssertBits(value, DoubleDouble.RootN(value, 1));
AssertBits(-value, DoubleDouble.RootN(-value, 1)); AssertBits(-value, DoubleDouble.RootN(-value, 1));
AssertBits(DDMath.Reciprocal(value), DoubleDouble.RootN(value, -1)); AssertBits(DoubleDouble.Reciprocal(value), DoubleDouble.RootN(value, -1));
AssertBits(DoubleDouble.Sqrt(value), DoubleDouble.RootN(value, 2)); AssertBits(DoubleDouble.Sqrt(value), DoubleDouble.RootN(value, 2));
AssertBits(DoubleDouble.InvSqrt(value), DoubleDouble.RootN(value, -2)); AssertBits(DoubleDouble.InvSqrt(value), DoubleDouble.RootN(value, -2));
AssertBits(DoubleDouble.Cbrt(value), DoubleDouble.RootN(value, 3)); AssertBits(DoubleDouble.Cbrt(value), DoubleDouble.RootN(value, 3));
@@ -153,7 +153,7 @@ public class DoubleDoubleRootNTests
} }
[Fact] [Fact]
public void GenericContractAndFacadeDispatchToTypeMembers() public void GenericContractDispatchesToTypeMembers()
{ {
DoubleDouble[] values = [DoubleDouble.Zero, DoubleDouble.NegativeZero, DoubleDouble.NaN, DoubleDouble[] values = [DoubleDouble.Zero, DoubleDouble.NegativeZero, DoubleDouble.NaN,
DoubleDouble.PositiveInfinity, DoubleDouble.NegativeInfinity, new(-2.0), new(double.Epsilon), DoubleDouble.PositiveInfinity, DoubleDouble.NegativeInfinity, new(-2.0), new(double.Epsilon),
@@ -167,7 +167,6 @@ public class DoubleDoubleRootNTests
{ {
DoubleDouble result = DoubleDouble.RootN(value, degree); DoubleDouble result = DoubleDouble.RootN(value, degree);
AssertBits(result, GenericRootN(value, degree)); AssertBits(result, GenericRootN(value, degree));
AssertBits(result, DDMath.RootN(value, degree));
} }
} }
} }
@@ -16,7 +16,7 @@ public class PreciseMathExpTests
// Inputs are exact stored binary64 sums, not rounded irrational constants. // Inputs are exact stored binary64 sums, not rounded irrational constants.
DoubleDouble input = DoubleDouble.FromComponents(high, low); DoubleDouble input = DoubleDouble.FromComponents(high, low);
(Units(input.High) + Units(input.Low)).ShouldBe(Units(high) + Units(low)); (Units(input.High) + Units(input.Low)).ShouldBe(Units(high) + Units(low));
DoubleDouble actual = DDMath.Exp(input); DoubleDouble actual = DoubleDouble.Exp(input);
DoubleDouble.IsCanonical(actual).ShouldBeTrue(); DoubleDouble.IsCanonical(actual).ShouldBeTrue();
if (overflow) if (overflow)
{ {
@@ -60,7 +60,7 @@ public class PreciseMathExpTests
{ {
// Independently round Decimal.exp(exact binary64 input)/2^-1074 to an // Independently round Decimal.exp(exact binary64 input)/2^-1074 to an
// integer at 180 digits. Pin more than the general one-subnormal bound. // integer at 180 digits. Pin more than the general one-subnormal bound.
DoubleDouble actual = DDMath.Exp(new DoubleDouble(input)); DoubleDouble actual = DoubleDouble.Exp(new DoubleDouble(input));
BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(expectedBits); BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(expectedBits);
BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L); BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L);
} }
@@ -70,7 +70,7 @@ public class PreciseMathExpTests
[InlineData(709.781)] [InlineData(709.781)]
public void LegacyFiniteResultsDoNotOverflowPrematurely(double input) public void LegacyFiniteResultsDoNotOverflowPrematurely(double input)
{ {
DoubleDouble actual = DDMath.Exp(new DoubleDouble(input)); DoubleDouble actual = DoubleDouble.Exp(new DoubleDouble(input));
DoubleDouble.IsFinite(actual).ShouldBeTrue(); DoubleDouble.IsFinite(actual).ShouldBeTrue();
actual.Low.ShouldNotBe(0.0); actual.Low.ShouldNotBe(0.0);
// Component-aware accuracy for these exact inputs is pinned by ReferenceCases. // Component-aware accuracy for these exact inputs is pinned by ReferenceCases.
@@ -83,7 +83,7 @@ public class PreciseMathExpTests
double.MaxValue, double.MinValue]; double.MaxValue, double.MinValue];
foreach (double input in inputs) foreach (double input in inputs)
{ {
DoubleDouble actual = DDMath.Exp(new DoubleDouble(input)); DoubleDouble actual = DoubleDouble.Exp(new DoubleDouble(input));
double expected = Math.Exp(input); double expected = Math.Exp(input);
BitConverter.DoubleToInt64Bits(actual.High).ShouldBe( BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(
BitConverter.DoubleToInt64Bits(double.IsNaN(expected) ? double.NaN : expected)); BitConverter.DoubleToInt64Bits(double.IsNaN(expected) ? double.NaN : expected));
@@ -11,7 +11,7 @@ public class PreciseMathInvSqrtTests
public void IrrationalInverseRootsRetainMoreThanBinary64Precision(double value) public void IrrationalInverseRootsRetainMoreThanBinary64Precision(double value)
{ {
DoubleDouble input = new(value); DoubleDouble input = new(value);
DoubleDouble actual = DDMath.InvSqrt(input); DoubleDouble actual = DoubleDouble.InvSqrt(input);
actual.Low.ShouldNotBe(0.0); actual.Low.ShouldNotBe(0.0);
AssertInvSqrtBound(input, actual); AssertInvSqrtBound(input, actual);
} }
@@ -23,7 +23,7 @@ public class PreciseMathInvSqrtTests
double.NaN, -1.0, -double.Epsilon, double.MinValue]; double.NaN, -1.0, -double.Epsilon, double.MinValue];
foreach (double value in values) foreach (double value in values)
{ {
DoubleDouble actual = DDMath.InvSqrt(new DoubleDouble(value)); DoubleDouble actual = DoubleDouble.InvSqrt(new DoubleDouble(value));
// Binary64 is an independent oracle only for these special values. // Binary64 is an independent oracle only for these special values.
double expected = 1.0 / Math.Sqrt(value); double expected = 1.0 / Math.Sqrt(value);
BitConverter.DoubleToInt64Bits(actual.High).ShouldBe( BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(
@@ -33,7 +33,7 @@ public class PreciseMathInvSqrtTests
} }
DoubleDouble negative = DoubleDouble.FromComponents(-1.0, Math.ScaleB(1.0, -54)); DoubleDouble negative = DoubleDouble.FromComponents(-1.0, Math.ScaleB(1.0, -54));
DoubleDouble.IsNaN(DDMath.InvSqrt(negative)).ShouldBeTrue(); DoubleDouble.IsNaN(DoubleDouble.InvSqrt(negative)).ShouldBeTrue();
} }
[Fact] [Fact]
@@ -42,7 +42,7 @@ public class PreciseMathInvSqrtTests
// 1/sqrt(2^(2k)) = 2^-k exactly, including the minimum subnormal input. // 1/sqrt(2^(2k)) = 2^-k exactly, including the minimum subnormal input.
for (int exponent = -1074; exponent <= 1022; exponent += 2) for (int exponent = -1074; exponent <= 1022; exponent += 2)
{ {
DoubleDouble actual = DDMath.InvSqrt(new DoubleDouble(Math.ScaleB(1.0, exponent))); DoubleDouble actual = DoubleDouble.InvSqrt(new DoubleDouble(Math.ScaleB(1.0, exponent)));
actual.High.ShouldBe(Math.ScaleB(1.0, -(exponent / 2))); actual.High.ShouldBe(Math.ScaleB(1.0, -(exponent / 2)));
BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L); BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L);
} }
@@ -54,7 +54,7 @@ public class PreciseMathInvSqrtTests
foreach (double low in new[] { Math.ScaleB(1.0, -54), -Math.ScaleB(1.0, -54) }) foreach (double low in new[] { Math.ScaleB(1.0, -54), -Math.ScaleB(1.0, -54) })
{ {
DoubleDouble input = DoubleDouble.FromComponents(1.0, low); DoubleDouble input = DoubleDouble.FromComponents(1.0, low);
DoubleDouble actual = DDMath.InvSqrt(input); DoubleDouble actual = DoubleDouble.InvSqrt(input);
actual.High.ShouldBe(1.0); actual.High.ShouldBe(1.0);
Math.Sign(actual.Low).ShouldBe(-Math.Sign(low)); Math.Sign(actual.Low).ShouldBe(-Math.Sign(low));
AssertInvSqrtBound(input, actual); AssertInvSqrtBound(input, actual);
@@ -78,7 +78,7 @@ public class PreciseMathInvSqrtTests
if (boundary > 0.0) if (boundary > 0.0)
{ {
DoubleDouble input = new(boundary); DoubleDouble input = new(boundary);
AssertInvSqrtBound(input, DDMath.InvSqrt(input)); AssertInvSqrtBound(input, DoubleDouble.InvSqrt(input));
} }
} }
foreach (double residual in new[] { 0.0, low, -low, double.Epsilon, -double.Epsilon }) foreach (double residual in new[] { 0.0, low, -low, double.Epsilon, -double.Epsilon })
@@ -86,7 +86,7 @@ public class PreciseMathInvSqrtTests
DoubleDouble input = DoubleDouble.FromComponents(high, residual); DoubleDouble input = DoubleDouble.FromComponents(high, residual);
if (input.High > 0.0) if (input.High > 0.0)
{ {
AssertInvSqrtBound(input, DDMath.InvSqrt(input)); AssertInvSqrtBound(input, DoubleDouble.InvSqrt(input));
} }
} }
} }
@@ -103,7 +103,7 @@ public class PreciseMathInvSqrtTests
DoubleDouble.FromComponents(0.0, double.Epsilon)]; DoubleDouble.FromComponents(0.0, double.Epsilon)];
foreach (DoubleDouble input in boundaries) foreach (DoubleDouble input in boundaries)
{ {
AssertInvSqrtBound(input, DDMath.InvSqrt(input)); AssertInvSqrtBound(input, DoubleDouble.InvSqrt(input));
} }
} }
@@ -123,7 +123,7 @@ public class PreciseMathInvSqrtTests
DoubleDouble input = DoubleDouble.FromComponents(high, low); DoubleDouble input = DoubleDouble.FromComponents(high, low);
if (DoubleDouble.IsFinite(input) && input.High > 0.0) if (DoubleDouble.IsFinite(input) && input.High > 0.0)
{ {
AssertInvSqrtBound(input, DDMath.InvSqrt(input)); AssertInvSqrtBound(input, DoubleDouble.InvSqrt(input));
} }
} }
} }
@@ -140,7 +140,7 @@ public class PreciseMathInvSqrtTests
foreach (double sign in new[] { -1.0, 1.0 }) foreach (double sign in new[] { -1.0, 1.0 })
{ {
double low = Math.ScaleB(sign, exponent); double low = Math.ScaleB(sign, exponent);
DoubleDouble actual = DDMath.InvSqrt(DoubleDouble.FromComponents(1.0, low)); DoubleDouble actual = DoubleDouble.InvSqrt(DoubleDouble.FromComponents(1.0, low));
actual.High.ShouldBe(1.0); actual.High.ShouldBe(1.0);
actual.Low.ShouldBe(Math.ScaleB(-sign, exponent - 1)); actual.Low.ShouldBe(Math.ScaleB(-sign, exponent - 1));
DoubleDouble.IsCanonical(actual).ShouldBeTrue(); DoubleDouble.IsCanonical(actual).ShouldBeTrue();
@@ -17,10 +17,10 @@ public class PreciseMathLogTests
foreach (double value in values) foreach (double value in values)
{ {
DoubleDouble expected = new(Math.Log(value)); DoubleDouble expected = new(Math.Log(value));
AssertBits(DDMath.Log(new DoubleDouble(value)), expected.High, 0.0); AssertBits(DoubleDouble.Log(new DoubleDouble(value)), expected.High, 0.0);
} }
DoubleDouble negative = DoubleDouble.FromComponents(-1.0, Math.ScaleB(1.0, -54)); DoubleDouble negative = DoubleDouble.FromComponents(-1.0, Math.ScaleB(1.0, -54));
AssertBits(DDMath.Log(negative), double.NaN, 0.0); AssertBits(DoubleDouble.Log(negative), double.NaN, 0.0);
} }
[Theory] [Theory]
@@ -31,7 +31,7 @@ public class PreciseMathLogTests
// 2200 digits and checked with Fraction. No binary64 log supplies the oracle. // 2200 digits and checked with Fraction. No binary64 log supplies the oracle.
DoubleDouble input = DoubleDouble.FromComponents(high, low); DoubleDouble input = DoubleDouble.FromComponents(high, low);
(Units(input.High) + Units(input.Low)).ShouldBe(Units(high) + Units(low)); (Units(input.High) + Units(input.Low)).ShouldBe(Units(high) + Units(low));
DoubleDouble actual = DDMath.Log(input); DoubleDouble actual = DoubleDouble.Log(input);
AssertReferenceBound(actual, reference, 1, $"{label}: ({high:R}, {low:R})"); AssertReferenceBound(actual, reference, 1, $"{label}: ({high:R}, {low:R})");
if (input == DoubleDouble.One) if (input == DoubleDouble.One)
{ {
@@ -51,7 +51,7 @@ public class PreciseMathLogTests
for (int exponent = -1074; exponent <= 1023; ++exponent) for (int exponent = -1074; exponent <= 1023; ++exponent)
{ {
DoubleDouble input = new(Math.ScaleB(1.0, exponent)); DoubleDouble input = new(Math.ScaleB(1.0, exponent));
AssertReferenceBound(DDMath.Log(input), LogReferenceData.Ln2, exponent, $"2^{exponent}"); AssertReferenceBound(DoubleDouble.Log(input), LogReferenceData.Ln2, exponent, $"2^{exponent}");
} }
} }
@@ -63,7 +63,7 @@ public class PreciseMathLogTests
{ {
// The legacy Log(1e308) formed u+value and overflowed. The finite // The legacy Log(1e308) formed u+value and overflowed. The finite
// kernel must reduce the mantissa before constructing its denominator. // kernel must reduce the mantissa before constructing its denominator.
DoubleDouble actual = DDMath.Log(new DoubleDouble(input)); DoubleDouble actual = DoubleDouble.Log(new DoubleDouble(input));
DoubleDouble.IsFinite(actual).ShouldBeTrue(); DoubleDouble.IsFinite(actual).ShouldBeTrue();
actual.Low.ShouldNotBe(0.0); actual.Low.ShouldNotBe(0.0);
DoubleDouble.IsCanonical(actual).ShouldBeTrue(); DoubleDouble.IsCanonical(actual).ShouldBeTrue();
@@ -75,12 +75,12 @@ public class PreciseMathLogTests
{ {
foreach (double sign in new[] { -1.0, 1.0 }) foreach (double sign in new[] { -1.0, 1.0 })
{ {
AssertBits(DDMath.Log(DoubleDouble.FromComponents(1.0, sign * double.Epsilon)), AssertBits(DoubleDouble.Log(DoubleDouble.FromComponents(1.0, sign * double.Epsilon)),
sign * double.Epsilon, 0.0); sign * double.Epsilon, 0.0);
double delta = sign * Math.ScaleB(1.0, -500); double delta = sign * Math.ScaleB(1.0, -500);
// ln(1+d) = d-d^2/2+O(d^3). Here the cubic tail is below half // ln(1+d) = d-d^2/2+O(d^3). Here the cubic tail is below half
// the minimum subnormal, so these two binary components are exact. // the minimum subnormal, so these two binary components are exact.
AssertBits(DDMath.Log(DoubleDouble.FromComponents(1.0, delta)), AssertBits(DoubleDouble.Log(DoubleDouble.FromComponents(1.0, delta)),
delta, -Math.ScaleB(1.0, -1001)); delta, -Math.ScaleB(1.0, -1001));
} }
} }
@@ -89,8 +89,8 @@ public class PreciseMathLogTests
public void ArbitraryPairsAreNormalizedByThePublicFactoryBeforeLog() public void ArbitraryPairsAreNormalizedByThePublicFactoryBeforeLog()
{ {
DoubleDouble positive = DoubleDouble.FromComponents(0.0, 2.0); DoubleDouble positive = DoubleDouble.FromComponents(0.0, 2.0);
AssertReferenceBound(DDMath.Log(positive), LogReferenceData.Ln2, 1, "zero-high public construction"); AssertReferenceBound(DoubleDouble.Log(positive), LogReferenceData.Ln2, 1, "zero-high public construction");
AssertBits(DDMath.Log(DoubleDouble.FromComponents(0.0, -2.0)), double.NaN, 0.0); AssertBits(DoubleDouble.Log(DoubleDouble.FromComponents(0.0, -2.0)), double.NaN, 0.0);
} }
private static void AssertReferenceBound(DoubleDouble actual, string reference, int factor, string context) private static void AssertReferenceBound(DoubleDouble actual, string reference, int factor, string context)
@@ -6,17 +6,17 @@ public class PreciseMathPowTests
public void NegativePowerDoesNotOverflowBeforeReciprocation() public void NegativePowerDoesNotOverflowBeforeReciprocation()
{ {
// Legacy review: Pow(2, -1024) overflowed before taking the reciprocal. // Legacy review: Pow(2, -1024) overflowed before taking the reciprocal.
DoubleDouble actual = DDMath.Pow(new DoubleDouble(2.0), -1024); DoubleDouble actual = DoubleDouble.Pow(new DoubleDouble(2.0), -1024);
AssertBits(actual, Math.ScaleB(1.0, -1024), 0.0); AssertBits(actual, Math.ScaleB(1.0, -1024), 0.0);
} }
[Fact] [Fact]
public void MinimumIntegerExponentDoesNotOverflowItsMagnitude() public void MinimumIntegerExponentDoesNotOverflowItsMagnitude()
{ {
AssertBits(DDMath.Pow(DoubleDouble.One, int.MinValue), 1.0, 0.0); AssertBits(DoubleDouble.Pow(DoubleDouble.One, int.MinValue), 1.0, 0.0);
AssertBits(DDMath.Pow(DoubleDouble.NegativeOne, int.MinValue), 1.0, 0.0); AssertBits(DoubleDouble.Pow(DoubleDouble.NegativeOne, int.MinValue), 1.0, 0.0);
AssertBits(DDMath.Pow(new DoubleDouble(2.0), int.MinValue), 0.0, 0.0); AssertBits(DoubleDouble.Pow(new DoubleDouble(2.0), int.MinValue), 0.0, 0.0);
AssertBits(DDMath.Pow(new DoubleDouble(0.5), int.MinValue), double.PositiveInfinity, 0.0); AssertBits(DoubleDouble.Pow(new DoubleDouble(0.5), int.MinValue), double.PositiveInfinity, 0.0);
} }
[Fact] [Fact]
@@ -25,7 +25,7 @@ public class PreciseMathPowTests
foreach (double sign in new[] { -1.0, 1.0 }) foreach (double sign in new[] { -1.0, 1.0 })
{ {
DoubleDouble input = DoubleDouble.FromComponents(sign * double.MaxValue, sign * double.Epsilon); DoubleDouble input = DoubleDouble.FromComponents(sign * double.MaxValue, sign * double.Epsilon);
AssertBits(DDMath.Pow(input, 1), input.High, input.Low); AssertBits(DoubleDouble.Pow(input, 1), input.High, input.Low);
} }
} }
@@ -42,7 +42,7 @@ public class PreciseMathPowTests
// Every int converts exactly to double. BCL is an independent // Every int converts exactly to double. BCL is an independent
// special-value/sign oracle, not the finite precision oracle. // special-value/sign oracle, not the finite precision oracle.
double expected = Math.Pow(input, exponent); double expected = Math.Pow(input, exponent);
AssertBits(DDMath.Pow(new DoubleDouble(input), exponent), AssertBits(DoubleDouble.Pow(new DoubleDouble(input), exponent),
double.IsNaN(expected) ? double.NaN : expected, 0.0); double.IsNaN(expected) ? double.NaN : expected, 0.0);
} }
} }
@@ -63,7 +63,7 @@ public class PreciseMathPowTests
foreach (double sign in new[] { -1.0, 1.0 }) foreach (double sign in new[] { -1.0, 1.0 })
{ {
double expected = sign < 0.0 && (power & 1) != 0 ? -magnitude : magnitude; double expected = sign < 0.0 && (power & 1) != 0 ? -magnitude : magnitude;
AssertBits(DDMath.Pow(new DoubleDouble(sign * Math.ScaleB(1.0, inputExponent)), power), AssertBits(DoubleDouble.Pow(new DoubleDouble(sign * Math.ScaleB(1.0, inputExponent)), power),
expected, 0.0); expected, 0.0);
} }
} }
@@ -118,13 +118,13 @@ public class PreciseMathPowTests
AssertExactPowerBound(input, -1); AssertExactPowerBound(input, -1);
if (scale == -500 || scale == 500) if (scale == -500 || scale == 500)
{ {
DoubleDouble actual = DDMath.Pow(input, -2); DoubleDouble actual = DoubleDouble.Pow(input, -2);
actual.Low.ShouldNotBe(0.0); actual.Low.ShouldNotBe(0.0);
AssertExactPowerBound(input, -2); AssertExactPowerBound(input, -2);
} }
if (scale < 0) if (scale < 0)
{ {
DDMath.Pow(input, -1).Low.ShouldNotBe(0.0); DoubleDouble.Pow(input, -1).Low.ShouldNotBe(0.0);
} }
} }
} }
@@ -141,11 +141,11 @@ public class PreciseMathPowTests
foreach (double sign in new[] { -1.0, 1.0 }) foreach (double sign in new[] { -1.0, 1.0 })
{ {
DoubleDouble tiny = DoubleDouble.FromComponents(Math.ScaleB(1.0, -500), sign * double.Epsilon); DoubleDouble tiny = DoubleDouble.FromComponents(Math.ScaleB(1.0, -500), sign * double.Epsilon);
AssertBits(DDMath.Pow(tiny, -1), Math.ScaleB(1.0, 500), -sign * Math.ScaleB(1.0, -74)); AssertBits(DoubleDouble.Pow(tiny, -1), Math.ScaleB(1.0, 500), -sign * Math.ScaleB(1.0, -74));
DoubleDouble huge = DoubleDouble.FromComponents(Math.ScaleB(1.0, 500), sign * Math.ScaleB(1.0, -74)); DoubleDouble huge = DoubleDouble.FromComponents(Math.ScaleB(1.0, 500), sign * Math.ScaleB(1.0, -74));
AssertBits(DDMath.Pow(huge, -1), Math.ScaleB(1.0, -500), -sign * double.Epsilon); AssertBits(DoubleDouble.Pow(huge, -1), Math.ScaleB(1.0, -500), -sign * double.Epsilon);
DoubleDouble tinySquare = DoubleDouble.FromComponents(Math.ScaleB(1.0, -250), sign * double.Epsilon); DoubleDouble tinySquare = DoubleDouble.FromComponents(Math.ScaleB(1.0, -250), sign * double.Epsilon);
AssertBits(DDMath.Pow(tinySquare, -2), Math.ScaleB(1.0, 500), -sign * Math.ScaleB(1.0, -323)); AssertBits(DoubleDouble.Pow(tinySquare, -2), Math.ScaleB(1.0, 500), -sign * Math.ScaleB(1.0, -323));
} }
} }
@@ -165,7 +165,7 @@ public class PreciseMathPowTests
{ {
DoubleDouble input = DoubleDouble.FromComponents(sign * Math.ScaleB(1.0, scale), DoubleDouble input = DoubleDouble.FromComponents(sign * Math.ScaleB(1.0, scale),
sign * lowSign * Math.ScaleB(1.0, scale - 54)); sign * lowSign * Math.ScaleB(1.0, scale - 54));
AssertBits(DDMath.Pow(input, exponent), sign * (aboveTie ? double.Epsilon : 0.0), 0.0); AssertBits(DoubleDouble.Pow(input, exponent), sign * (aboveTie ? double.Epsilon : 0.0), 0.0);
} }
} }
} }
@@ -183,7 +183,7 @@ public class PreciseMathPowTests
AssertExactPowerBound(finite, exponent); AssertExactPowerBound(finite, exponent);
DoubleDouble overflowing = DoubleDouble.FromComponents(Math.ScaleB(1.0, scale), DoubleDouble overflowing = DoubleDouble.FromComponents(Math.ScaleB(1.0, scale),
lowSign * Math.ScaleB(1.0, scale - 57)); lowSign * Math.ScaleB(1.0, scale - 57));
AssertBits(DDMath.Pow(overflowing, exponent), double.PositiveInfinity, 0.0); AssertBits(DoubleDouble.Pow(overflowing, exponent), double.PositiveInfinity, 0.0);
} }
} }
@@ -208,7 +208,7 @@ public class PreciseMathPowTests
foreach (double sign in new[] { -1.0, 1.0 }) foreach (double sign in new[] { -1.0, 1.0 })
{ {
DoubleDouble input = DoubleDouble.FromComponents(sign * high, sign * low); DoubleDouble input = DoubleDouble.FromComponents(sign * high, sign * low);
DoubleDouble actual = DDMath.Pow(input, exponent); DoubleDouble actual = DoubleDouble.Pow(input, exponent);
actual.Low.ShouldNotBe(0.0); actual.Low.ShouldNotBe(0.0);
int resultSign = sign < 0.0 && (exponent & 1) != 0 ? -1 : 1; int resultSign = sign < 0.0 && (exponent & 1) != 0 ? -1 : 1;
// Check both interval endpoints so the reference uncertainty is not // Check both interval endpoints so the reference uncertainty is not
@@ -239,7 +239,7 @@ public class PreciseMathPowTests
{ {
(numerator, denominator) = (denominator * numerator.Sign, BigInteger.Abs(numerator)); (numerator, denominator) = (denominator * numerator.Sign, BigInteger.Abs(numerator));
} }
AssertRationalBound(DDMath.Pow(input, exponent), numerator, denominator, exponent, AssertRationalBound(DoubleDouble.Pow(input, exponent), numerator, denominator, exponent,
$"({input.High:R}, {input.Low:R})^{exponent}"); $"({input.High:R}, {input.Low:R})^{exponent}");
} }
@@ -10,10 +10,10 @@ public class PreciseMathRealPowSpecialTests
DoubleDouble.FromComponents(double.MaxValue, double.Epsilon)]; DoubleDouble.FromComponents(double.MaxValue, double.Epsilon)];
foreach (DoubleDouble value in values) foreach (DoubleDouble value in values)
{ {
AssertBits(DDMath.Pow(value, 0.0), DoubleDouble.One); AssertBits(DoubleDouble.Pow(value, 0.0), DoubleDouble.One);
AssertBits(DDMath.Pow(value, new DoubleDouble(-0.0)), DoubleDouble.One); AssertBits(DoubleDouble.Pow(value, new DoubleDouble(-0.0)), DoubleDouble.One);
AssertBits(DDMath.Pow(value, 1.0), value); AssertBits(DoubleDouble.Pow(value, 1.0), value);
AssertBits(DDMath.Pow(value, DoubleDouble.One), value); AssertBits(DoubleDouble.Pow(value, DoubleDouble.One), value);
} }
} }
@@ -37,8 +37,8 @@ public class PreciseMathRealPowSpecialTests
continue; continue;
} }
DoubleDouble input = new(value); DoubleDouble input = new(value);
AssertBits(DDMath.Pow(input, exponent), new DoubleDouble(expected)); AssertBits(DoubleDouble.Pow(input, exponent), new DoubleDouble(expected));
AssertBits(DDMath.Pow(input, new DoubleDouble(exponent)), new DoubleDouble(expected)); AssertBits(DoubleDouble.Pow(input, new DoubleDouble(exponent)), new DoubleDouble(expected));
} }
} }
} }
@@ -53,27 +53,27 @@ public class PreciseMathRealPowSpecialTests
DoubleDouble odd = DoubleDouble.FromComponents(direction * high, 1.0); DoubleDouble odd = DoubleDouble.FromComponents(direction * high, 1.0);
DoubleDouble even = DoubleDouble.FromComponents(direction * high, 2.0); DoubleDouble even = DoubleDouble.FromComponents(direction * high, 2.0);
DoubleDouble fractional = DoubleDouble.FromComponents(direction * high, 0.5); DoubleDouble fractional = DoubleDouble.FromComponents(direction * high, 0.5);
AssertBits(DDMath.Pow(DoubleDouble.NegativeOne, odd), DoubleDouble.NegativeOne); AssertBits(DoubleDouble.Pow(DoubleDouble.NegativeOne, odd), DoubleDouble.NegativeOne);
AssertBits(DDMath.Pow(DoubleDouble.NegativeOne, even), DoubleDouble.One); AssertBits(DoubleDouble.Pow(DoubleDouble.NegativeOne, even), DoubleDouble.One);
AssertBits(DDMath.Pow(DoubleDouble.NegativeOne, fractional), DoubleDouble.NaN); AssertBits(DoubleDouble.Pow(DoubleDouble.NegativeOne, fractional), DoubleDouble.NaN);
AssertBits(DDMath.Pow(new DoubleDouble(-2.0), odd), AssertBits(DoubleDouble.Pow(new DoubleDouble(-2.0), odd),
new DoubleDouble(direction > 0.0 ? double.NegativeInfinity : -0.0)); new DoubleDouble(direction > 0.0 ? double.NegativeInfinity : -0.0));
AssertBits(DDMath.Pow(new DoubleDouble(-0.5), odd), AssertBits(DoubleDouble.Pow(new DoubleDouble(-0.5), odd),
new DoubleDouble(direction > 0.0 ? -0.0 : double.NegativeInfinity)); new DoubleDouble(direction > 0.0 ? -0.0 : double.NegativeInfinity));
AssertBits(DDMath.Pow(new DoubleDouble(-2.0), even), AssertBits(DoubleDouble.Pow(new DoubleDouble(-2.0), even),
new DoubleDouble(direction > 0.0 ? double.PositiveInfinity : 0.0)); new DoubleDouble(direction > 0.0 ? double.PositiveInfinity : 0.0));
AssertBits(DDMath.Pow(new DoubleDouble(-0.0), odd), AssertBits(DoubleDouble.Pow(new DoubleDouble(-0.0), odd),
new DoubleDouble(direction > 0.0 ? -0.0 : double.NegativeInfinity)); new DoubleDouble(direction > 0.0 ? -0.0 : double.NegativeInfinity));
AssertBits(DDMath.Pow(new DoubleDouble(double.NegativeInfinity), odd), AssertBits(DoubleDouble.Pow(new DoubleDouble(double.NegativeInfinity), odd),
new DoubleDouble(direction > 0.0 ? double.NegativeInfinity : -0.0)); new DoubleDouble(direction > 0.0 ? double.NegativeInfinity : -0.0));
AssertBits(DDMath.Pow(new DoubleDouble(-0.0), fractional), AssertBits(DoubleDouble.Pow(new DoubleDouble(-0.0), fractional),
new DoubleDouble(direction > 0.0 ? 0.0 : double.PositiveInfinity)); new DoubleDouble(direction > 0.0 ? 0.0 : double.PositiveInfinity));
AssertBits(DDMath.Pow(new DoubleDouble(double.NegativeInfinity), fractional), AssertBits(DoubleDouble.Pow(new DoubleDouble(double.NegativeInfinity), fractional),
new DoubleDouble(direction > 0.0 ? double.PositiveInfinity : 0.0)); new DoubleDouble(direction > 0.0 ? double.PositiveInfinity : 0.0));
} }
} }
DoubleDouble belowInteger = DoubleDouble.FromComponents(3.0, -double.Epsilon); DoubleDouble belowInteger = DoubleDouble.FromComponents(3.0, -double.Epsilon);
AssertBits(DDMath.Pow(new DoubleDouble(-2.0), belowInteger), DoubleDouble.NaN); AssertBits(DoubleDouble.Pow(new DoubleDouble(-2.0), belowInteger), DoubleDouble.NaN);
} }
[Fact] [Fact]
@@ -84,9 +84,9 @@ public class PreciseMathRealPowSpecialTests
foreach (double sign in new[] { -1.0, 1.0 }) foreach (double sign in new[] { -1.0, 1.0 })
{ {
DoubleDouble value = DoubleDouble.FromComponents(sign, sign * low); DoubleDouble value = DoubleDouble.FromComponents(sign, sign * low);
AssertBits(DDMath.Pow(value, double.PositiveInfinity), AssertBits(DoubleDouble.Pow(value, double.PositiveInfinity),
new DoubleDouble(low > 0.0 ? double.PositiveInfinity : 0.0)); new DoubleDouble(low > 0.0 ? double.PositiveInfinity : 0.0));
AssertBits(DDMath.Pow(value, new DoubleDouble(double.NegativeInfinity)), AssertBits(DoubleDouble.Pow(value, new DoubleDouble(double.NegativeInfinity)),
new DoubleDouble(low < 0.0 ? double.PositiveInfinity : 0.0)); new DoubleDouble(low < 0.0 ? double.PositiveInfinity : 0.0));
} }
} }
@@ -102,11 +102,11 @@ public class PreciseMathRealPowSpecialTests
{ {
foreach (int exponent in exponents) foreach (int exponent in exponents)
{ {
DoubleDouble expected = DDMath.Pow(value, exponent); DoubleDouble expected = DoubleDouble.Pow(value, exponent);
// Dispatch consistency, not an accuracy oracle: the integer suite // Dispatch consistency, not an accuracy oracle: the integer suite
// separately checks these algorithms against independent references. // separately checks these algorithms against independent references.
AssertBits(DDMath.Pow(value, (double)exponent), expected); AssertBits(DoubleDouble.Pow(value, (double)exponent), expected);
AssertBits(DDMath.Pow(value, new DoubleDouble(exponent)), expected); AssertBits(DoubleDouble.Pow(value, new DoubleDouble(exponent)), expected);
} }
} }
} }
@@ -116,8 +116,8 @@ public class PreciseMathRealPowSpecialTests
{ {
// Independently rounded sqrt(2) binary64 components; also used by the // Independently rounded sqrt(2) binary64 components; also used by the
// constant-reference generator. Compare components, not high+low in double. // constant-reference generator. Compare components, not high+low in double.
foreach (DoubleDouble actual in new[] { DDMath.Pow(new DoubleDouble(2.0), 0.5), foreach (DoubleDouble actual in new[] { DoubleDouble.Pow(new DoubleDouble(2.0), 0.5),
DDMath.Pow(new DoubleDouble(2.0), new DoubleDouble(0.5)) }) DoubleDouble.Pow(new DoubleDouble(2.0), new DoubleDouble(0.5)) })
{ {
actual.High.ShouldBe(1.4142135623730951); actual.High.ShouldBe(1.4142135623730951);
Math.Abs(actual.Low - -9.667293313452913e-17).ShouldBeLessThan(Math.ScaleB(1.0, -90)); Math.Abs(actual.Low - -9.667293313452913e-17).ShouldBeLessThan(Math.ScaleB(1.0, -90));
@@ -132,8 +132,8 @@ public class PreciseMathRealPowSpecialTests
{ {
// 4^-537.25 = 2^-1074.5 lies above the half-epsilon threshold; // 4^-537.25 = 2^-1074.5 lies above the half-epsilon threshold;
// 4^-537.75 = 2^-1075.5 lies below it. These are well separated from ties. // 4^-537.75 = 2^-1075.5 lies below it. These are well separated from ties.
foreach (DoubleDouble actual in new[] { DDMath.Pow(new DoubleDouble(4.0), exponent), foreach (DoubleDouble actual in new[] { DoubleDouble.Pow(new DoubleDouble(4.0), exponent),
DDMath.Pow(new DoubleDouble(4.0), new DoubleDouble(exponent)) }) DoubleDouble.Pow(new DoubleDouble(4.0), new DoubleDouble(exponent)) })
{ {
BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(expectedBits); BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(expectedBits);
BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L); BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L);
@@ -149,8 +149,8 @@ public class PreciseMathRealPowSpecialTests
{ {
bool grows = (value > 1.0) == (exponent > 0.0); bool grows = (value > 1.0) == (exponent > 0.0);
DoubleDouble expected = new(grows ? double.PositiveInfinity : 0.0); DoubleDouble expected = new(grows ? double.PositiveInfinity : 0.0);
AssertBits(DDMath.Pow(new DoubleDouble(value), exponent), expected); AssertBits(DoubleDouble.Pow(new DoubleDouble(value), exponent), expected);
AssertBits(DDMath.Pow(new DoubleDouble(value), new DoubleDouble(exponent)), expected); AssertBits(DoubleDouble.Pow(new DoubleDouble(value), new DoubleDouble(exponent)), expected);
} }
} }
} }
@@ -21,7 +21,7 @@ public class PreciseMathRealPowTests
(Units(input.High) + Units(input.Low)).ShouldBe(Units(high) + Units(low)); (Units(input.High) + Units(input.Low)).ShouldBe(Units(high) + Units(low));
BigInteger exponentUnits = Units(exponentHigh) + Units(scalar ? 0.0 : exponentLow); BigInteger exponentUnits = Units(exponentHigh) + Units(scalar ? 0.0 : exponentLow);
(Units(exponent.High) + Units(exponent.Low)).ShouldBe(exponentUnits); (Units(exponent.High) + Units(exponent.Low)).ShouldBe(exponentUnits);
DoubleDouble actual = scalar ? DDMath.Pow(input, exponentHigh) : DDMath.Pow(input, exponent); DoubleDouble actual = scalar ? DoubleDouble.Pow(input, exponentHigh) : DoubleDouble.Pow(input, exponent);
string context = $"{label}, scalar={scalar}: ({high:R}, {low:R})^({exponentHigh:R}, {exponentLow:R})"; string context = $"{label}, scalar={scalar}: ({high:R}, {low:R})^({exponentHigh:R}, {exponentLow:R})";
AssertReferenceBound(actual, reference, exponentUnits, context); AssertReferenceBound(actual, reference, exponentUnits, context);
} }
@@ -36,9 +36,9 @@ public class PreciseMathRealPowTests
// These references are exact rational identities. The general ln/exp // These references are exact rational identities. The general ln/exp
// implementation is approximate, so do not require universal bit equality. // implementation is approximate, so do not require universal bit equality.
string reference = expected.ToString("e17", CultureInfo.InvariantCulture); string reference = expected.ToString("e17", CultureInfo.InvariantCulture);
AssertReferenceBound(DDMath.Pow(new DoubleDouble(input), exponent), reference, AssertReferenceBound(DoubleDouble.Pow(new DoubleDouble(input), exponent), reference,
Units(exponent), "scalar exact binary reference"); Units(exponent), "scalar exact binary reference");
AssertReferenceBound(DDMath.Pow(new DoubleDouble(input), new DoubleDouble(exponent)), reference, AssertReferenceBound(DoubleDouble.Pow(new DoubleDouble(input), new DoubleDouble(exponent)), reference,
Units(exponent), "DD exact binary reference"); Units(exponent), "DD exact binary reference");
} }
@@ -55,17 +55,17 @@ public class PreciseMathRealPowTests
{ {
foreach (int exponent in new[] { -63, -7, -1, 0, 1, 7, 63 }) foreach (int exponent in new[] { -63, -7, -1, 0, 1, 7, 63 })
{ {
DoubleDouble expected = DDMath.Pow(input, exponent); DoubleDouble expected = DoubleDouble.Pow(input, exponent);
AssertBits(DDMath.Pow(input, (double)exponent), expected); AssertBits(DoubleDouble.Pow(input, (double)exponent), expected);
AssertBits(DDMath.Pow(input, new DoubleDouble(exponent)), expected); AssertBits(DoubleDouble.Pow(input, new DoubleDouble(exponent)), expected);
} }
} }
foreach (int exponent in new[] { int.MinValue, int.MaxValue }) foreach (int exponent in new[] { int.MinValue, int.MaxValue })
{ {
DoubleDouble input = DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -54)); DoubleDouble input = DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -54));
DoubleDouble expected = DDMath.Pow(input, exponent); DoubleDouble expected = DoubleDouble.Pow(input, exponent);
AssertBits(DDMath.Pow(input, (double)exponent), expected); AssertBits(DoubleDouble.Pow(input, (double)exponent), expected);
AssertBits(DDMath.Pow(input, new DoubleDouble(exponent)), expected); AssertBits(DoubleDouble.Pow(input, new DoubleDouble(exponent)), expected);
} }
} }
@@ -77,8 +77,8 @@ public class PreciseMathRealPowTests
foreach (double low in new[] { double.Epsilon, -double.Epsilon }) foreach (double low in new[] { double.Epsilon, -double.Epsilon })
{ {
DoubleDouble input = DoubleDouble.FromComponents(high, low); DoubleDouble input = DoubleDouble.FromComponents(high, low);
AssertBits(DDMath.Pow(input, 1.0), input); AssertBits(DoubleDouble.Pow(input, 1.0), input);
AssertBits(DDMath.Pow(input, DoubleDouble.One), input); AssertBits(DoubleDouble.Pow(input, DoubleDouble.One), input);
} }
} }
} }
@@ -12,7 +12,7 @@ public class PreciseMathReciprocalTests
foreach (double low in new[] { 0.0, Math.ScaleB(1.0, -54), -Math.ScaleB(1.0, -54) }) foreach (double low in new[] { 0.0, Math.ScaleB(1.0, -54), -Math.ScaleB(1.0, -54) })
{ {
DoubleDouble input = DoubleDouble.FromComponents(high, low); DoubleDouble input = DoubleDouble.FromComponents(high, low);
DoubleDouble actual = DDMath.Reciprocal(input); DoubleDouble actual = DoubleDouble.Reciprocal(input);
actual.Low.ShouldNotBe(0.0); actual.Low.ShouldNotBe(0.0);
AssertReciprocalBound(input, actual); AssertReciprocalBound(input, actual);
} }
@@ -25,7 +25,7 @@ public class PreciseMathReciprocalTests
double.NaN, BitConverter.Int64BitsToDouble(0x7ff0000000000001L)]; double.NaN, BitConverter.Int64BitsToDouble(0x7ff0000000000001L)];
foreach (double value in values) foreach (double value in values)
{ {
DoubleDouble actual = DDMath.Reciprocal(new DoubleDouble(value)); DoubleDouble actual = DoubleDouble.Reciprocal(new DoubleDouble(value));
// Binary64 is an independent oracle for special values only. // Binary64 is an independent oracle for special values only.
double expected = 1.0 / value; double expected = 1.0 / value;
BitConverter.DoubleToInt64Bits(actual.High).ShouldBe( BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(
@@ -44,7 +44,7 @@ public class PreciseMathReciprocalTests
{ {
foreach (double sign in new[] { -1.0, 1.0 }) foreach (double sign in new[] { -1.0, 1.0 })
{ {
DoubleDouble actual = DDMath.Reciprocal(new DoubleDouble(Math.ScaleB(sign, exponent))); DoubleDouble actual = DoubleDouble.Reciprocal(new DoubleDouble(Math.ScaleB(sign, exponent)));
actual.High.ShouldBe(Math.ScaleB(sign, -exponent)); actual.High.ShouldBe(Math.ScaleB(sign, -exponent));
BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L); BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L);
DoubleDouble.IsCanonical(actual).ShouldBeTrue(); DoubleDouble.IsCanonical(actual).ShouldBeTrue();
@@ -104,7 +104,7 @@ public class PreciseMathReciprocalTests
{ {
double low = Math.ScaleB(lowSign, exponent); double low = Math.ScaleB(lowSign, exponent);
DoubleDouble input = DoubleDouble.FromComponents(sign, low); DoubleDouble input = DoubleDouble.FromComponents(sign, low);
DoubleDouble actual = DDMath.Reciprocal(input); DoubleDouble actual = DoubleDouble.Reciprocal(input);
actual.High.ShouldBe(sign); actual.High.ShouldBe(sign);
actual.Low.ShouldBe(-low); actual.Low.ShouldBe(-low);
AssertReciprocalBound(input, actual); AssertReciprocalBound(input, actual);
@@ -116,7 +116,7 @@ public class PreciseMathReciprocalTests
private static void AssertMatchesDivision(DoubleDouble input) private static void AssertMatchesDivision(DoubleDouble input)
{ {
DoubleDouble expected = 1.0 / input; DoubleDouble expected = 1.0 / input;
DoubleDouble actual = DDMath.Reciprocal(input); DoubleDouble actual = DoubleDouble.Reciprocal(input);
string context = $"Input: ({input.High:R}, {input.Low:R})"; string context = $"Input: ({input.High:R}, {input.Low:R})";
BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(BitConverter.DoubleToInt64Bits(expected.High), context); BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(BitConverter.DoubleToInt64Bits(expected.High), context);
BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(BitConverter.DoubleToInt64Bits(expected.Low), context); BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(BitConverter.DoubleToInt64Bits(expected.Low), context);
@@ -10,7 +10,7 @@ public class PreciseMathSqrtTests
[InlineData(2.25, 1.5)] [InlineData(2.25, 1.5)]
public void ExactBinarySquaresHaveExactRoots(double square, double root) public void ExactBinarySquaresHaveExactRoots(double square, double root)
{ {
DoubleDouble actual = DDMath.Sqrt(new DoubleDouble(square)); DoubleDouble actual = DoubleDouble.Sqrt(new DoubleDouble(square));
actual.High.ShouldBe(root); actual.High.ShouldBe(root);
BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L); BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L);
} }
@@ -21,7 +21,7 @@ public class PreciseMathSqrtTests
// sqrt(2^(2k)) = 2^k exactly, including the minimum subnormal input. // sqrt(2^(2k)) = 2^k exactly, including the minimum subnormal input.
for (int exponent = -1074; exponent <= 1022; exponent += 2) for (int exponent = -1074; exponent <= 1022; exponent += 2)
{ {
DoubleDouble actual = DDMath.Sqrt(new DoubleDouble(Math.ScaleB(1.0, exponent))); DoubleDouble actual = DoubleDouble.Sqrt(new DoubleDouble(Math.ScaleB(1.0, exponent)));
actual.High.ShouldBe(Math.ScaleB(1.0, exponent / 2)); actual.High.ShouldBe(Math.ScaleB(1.0, exponent / 2));
BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L); BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L);
} }
@@ -36,7 +36,7 @@ public class PreciseMathSqrtTests
public void IrrationalRootsRetainMoreThanBinary64Precision(double value) public void IrrationalRootsRetainMoreThanBinary64Precision(double value)
{ {
DoubleDouble input = new(value); DoubleDouble input = new(value);
DoubleDouble actual = DDMath.Sqrt(input); DoubleDouble actual = DoubleDouble.Sqrt(input);
actual.Low.ShouldNotBe(0.0); actual.Low.ShouldNotBe(0.0);
AssertSqrtBound(input, actual); AssertSqrtBound(input, actual);
} }
@@ -47,7 +47,7 @@ public class PreciseMathSqrtTests
foreach (double low in new[] { Math.ScaleB(1.0, -54), -Math.ScaleB(1.0, -54) }) foreach (double low in new[] { Math.ScaleB(1.0, -54), -Math.ScaleB(1.0, -54) })
{ {
DoubleDouble input = DoubleDouble.FromComponents(1.0, low); DoubleDouble input = DoubleDouble.FromComponents(1.0, low);
DoubleDouble actual = DDMath.Sqrt(input); DoubleDouble actual = DoubleDouble.Sqrt(input);
actual.High.ShouldBe(1.0); actual.High.ShouldBe(1.0);
(Math.Sign(actual.Low) == Math.Sign(low)).ShouldBeTrue(); (Math.Sign(actual.Low) == Math.Sign(low)).ShouldBeTrue();
AssertSqrtBound(input, actual); AssertSqrtBound(input, actual);
@@ -61,7 +61,7 @@ public class PreciseMathSqrtTests
double.NaN, -1.0, -double.Epsilon, double.MinValue]; double.NaN, -1.0, -double.Epsilon, double.MinValue];
foreach (double value in values) foreach (double value in values)
{ {
DoubleDouble actual = DDMath.Sqrt(new DoubleDouble(value)); DoubleDouble actual = DoubleDouble.Sqrt(new DoubleDouble(value));
double expected = Math.Sqrt(value); double expected = Math.Sqrt(value);
// Construction canonicalizes all NaNs to double.NaN, including // Construction canonicalizes all NaNs to double.NaN, including
// domain errors; zeros must instead be compared by their sign bits. // domain errors; zeros must instead be compared by their sign bits.
@@ -72,7 +72,7 @@ public class PreciseMathSqrtTests
} }
DoubleDouble negative = DoubleDouble.FromComponents(-1.0, Math.ScaleB(1.0, -54)); DoubleDouble negative = DoubleDouble.FromComponents(-1.0, Math.ScaleB(1.0, -54));
DoubleDouble.IsNaN(DDMath.Sqrt(negative)).ShouldBeTrue(); DoubleDouble.IsNaN(DoubleDouble.Sqrt(negative)).ShouldBeTrue();
} }
[Fact] [Fact]
@@ -92,7 +92,7 @@ public class PreciseMathSqrtTests
if (boundary > 0.0) if (boundary > 0.0)
{ {
DoubleDouble input = new(boundary); DoubleDouble input = new(boundary);
AssertSqrtBound(input, DDMath.Sqrt(input)); AssertSqrtBound(input, DoubleDouble.Sqrt(input));
} }
} }
foreach (double residual in new[] { 0.0, low, -low, double.Epsilon, -double.Epsilon }) foreach (double residual in new[] { 0.0, low, -low, double.Epsilon, -double.Epsilon })
@@ -100,7 +100,7 @@ public class PreciseMathSqrtTests
DoubleDouble input = DoubleDouble.FromComponents(high, residual); DoubleDouble input = DoubleDouble.FromComponents(high, residual);
if (input.High > 0.0) if (input.High > 0.0)
{ {
AssertSqrtBound(input, DDMath.Sqrt(input)); AssertSqrtBound(input, DoubleDouble.Sqrt(input));
} }
} }
} }
@@ -120,7 +120,7 @@ public class PreciseMathSqrtTests
DoubleDouble.FromComponents(4.0, -Math.ScaleB(1.0, -52))]; DoubleDouble.FromComponents(4.0, -Math.ScaleB(1.0, -52))];
foreach (DoubleDouble input in boundaries) foreach (DoubleDouble input in boundaries)
{ {
AssertSqrtBound(input, DDMath.Sqrt(input)); AssertSqrtBound(input, DoubleDouble.Sqrt(input));
} }
} }
@@ -142,7 +142,7 @@ public class PreciseMathSqrtTests
DoubleDouble input = DoubleDouble.FromComponents(Math.ScaleB(high, exponent), Math.ScaleB(low, exponent)); DoubleDouble input = DoubleDouble.FromComponents(Math.ScaleB(high, exponent), Math.ScaleB(low, exponent));
if (DoubleDouble.IsFinite(input) && input.High > 0.0) if (DoubleDouble.IsFinite(input) && input.High > 0.0)
{ {
AssertSqrtBound(input, DDMath.Sqrt(input)); AssertSqrtBound(input, DoubleDouble.Sqrt(input));
} }
} }
} }
@@ -156,9 +156,9 @@ public class PreciseMathSqrtTests
foreach (double low in new[] { double.Epsilon, 1e-308, 2.0, double.MaxValue }) foreach (double low in new[] { double.Epsilon, 1e-308, 2.0, double.MaxValue })
{ {
DoubleDouble input = DoubleDouble.FromComponents(0.0, low); DoubleDouble input = DoubleDouble.FromComponents(0.0, low);
AssertSqrtBound(input, DDMath.Sqrt(input)); AssertSqrtBound(input, DoubleDouble.Sqrt(input));
} }
DoubleDouble.IsNaN(DDMath.Sqrt(DoubleDouble.FromComponents(0.0, -1.0))).ShouldBeTrue(); DoubleDouble.IsNaN(DoubleDouble.Sqrt(DoubleDouble.FromComponents(0.0, -1.0))).ShouldBeTrue();
} }
[Fact] [Fact]
@@ -177,7 +177,7 @@ public class PreciseMathSqrtTests
DoubleDouble input = DoubleDouble.FromComponents(high, low); DoubleDouble input = DoubleDouble.FromComponents(high, low);
if (DoubleDouble.IsFinite(input) && input.High > 0.0) if (DoubleDouble.IsFinite(input) && input.High > 0.0)
{ {
AssertSqrtBound(input, DDMath.Sqrt(input)); AssertSqrtBound(input, DoubleDouble.Sqrt(input));
} }
} }
} }
@@ -194,7 +194,7 @@ public class PreciseMathSqrtTests
foreach (double sign in new[] { -1.0, 1.0 }) foreach (double sign in new[] { -1.0, 1.0 })
{ {
double low = Math.ScaleB(sign, exponent); double low = Math.ScaleB(sign, exponent);
DoubleDouble actual = DDMath.Sqrt(DoubleDouble.FromComponents(1.0, low)); DoubleDouble actual = DoubleDouble.Sqrt(DoubleDouble.FromComponents(1.0, low));
actual.High.ShouldBe(1.0); actual.High.ShouldBe(1.0);
actual.Low.ShouldBe(Math.ScaleB(sign, exponent - 1)); actual.Low.ShouldBe(Math.ScaleB(sign, exponent - 1));
DoubleDouble.IsCanonical(actual).ShouldBeTrue(); DoubleDouble.IsCanonical(actual).ShouldBeTrue();
@@ -358,9 +358,6 @@ public class PreciseMathSqrtTests
DoubleDouble actual = DoubleDouble.Sqrt(input); DoubleDouble actual = DoubleDouble.Sqrt(input);
actual.High.ShouldBe(high); actual.High.ShouldBe(high);
BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(BitConverter.DoubleToInt64Bits(low)); BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(BitConverter.DoubleToInt64Bits(low));
DoubleDouble facade = DDMath.Sqrt(input);
BitConverter.DoubleToInt64Bits(facade.High).ShouldBe(BitConverter.DoubleToInt64Bits(actual.High));
BitConverter.DoubleToInt64Bits(facade.Low).ShouldBe(BitConverter.DoubleToInt64Bits(actual.Low));
AssertSqrtBound(input, actual); AssertSqrtBound(input, actual);
} }
@@ -10,7 +10,7 @@ public class PreciseMathTests
foreach (double scalar in scalars) foreach (double scalar in scalars)
{ {
// Binary64 is an independent reference when the input has no residual. // Binary64 is an independent reference when the input has no residual.
DoubleDouble actual = DDMath.Abs(new DoubleDouble(scalar)); DoubleDouble actual = DoubleDouble.Abs(new DoubleDouble(scalar));
AssertBits(actual, double.IsNaN(scalar) ? double.NaN : Math.Abs(scalar), 0.0); AssertBits(actual, double.IsNaN(scalar) ? double.NaN : Math.Abs(scalar), 0.0);
DoubleDouble.IsCanonical(actual).ShouldBeTrue(); DoubleDouble.IsCanonical(actual).ShouldBeTrue();
} }
@@ -19,8 +19,8 @@ public class PreciseMathTests
{ {
DoubleDouble positive = DoubleDouble.FromComponents(1.0, low); DoubleDouble positive = DoubleDouble.FromComponents(1.0, low);
DoubleDouble negative = DoubleDouble.FromComponents(-1.0, -low); DoubleDouble negative = DoubleDouble.FromComponents(-1.0, -low);
AssertBits(DDMath.Abs(positive), 1.0, low); AssertBits(DoubleDouble.Abs(positive), 1.0, low);
AssertBits(DDMath.Abs(negative), 1.0, low); AssertBits(DoubleDouble.Abs(negative), 1.0, low);
} }
} }
+48 -21
View File
@@ -22,8 +22,11 @@ already-normalized results. Mathematical constants use precomputed high/low pair
checked against independently computed high-precision values; accessing them does checked against independently computed high-precision values; accessing them does
not perform double-double arithmetic or allocate on the heap. not perform double-double arithmetic or allocate on the heap.
- Arithmetic: unary `+`/`-`, binary `+`, `-`, `*`, `/`, and both operand orders with - Arithmetic: static `Add`, `Subtract`, `Multiply`, and `Divide` methods accept
a `double`. Addition retains residuals under cancellation; multiplication uses `DoubleDouble` pairs and both operand orders with a `double`. `Plus` and `Negate`
provide unary operations. The arithmetic operators (`+`, `-`, `*`, `/`, unary
`+`/`-`, and `++`/`--`) are thin wrappers around these methods.
Addition retains residuals under cancellation; multiplication uses
fused multiply-add; division uses residual corrections. Mixed `double` operators fused multiply-add; division uses residual corrections. Mixed `double` operators
use specialized scalar paths rather than promoting the scalar to `DoubleDouble`. use specialized scalar paths rather than promoting the scalar to `DoubleDouble`.
Division retains the complete first remainder and applies two quotient corrections, Division retains the complete first remainder and applies two quotient corrections,
@@ -92,7 +95,7 @@ DoubleDouble quarterTurn = DoubleDouble.PiOver2;
The factors avoid deriving constants at runtime; the multiplication itself remains The factors avoid deriving constants at runtime; the multiplication itself remains
approximate DD arithmetic, so conversions are not guaranteed exact round trips. approximate DD arithmetic, so conversions are not guaranteed exact round trips.
The list is mathematical and dimensionless, not a table of unit-dependent physical The list is mathematical and dimensionless, not a table of unit-dependent physical
constants. The natural logarithm function is provided separately by `DDMath.Log`. constants. The natural logarithm function is provided separately by `DoubleDouble.Log`.
## Mathematical functions ## Mathematical functions
@@ -102,16 +105,40 @@ The additional `InvSqrt` helper and three-argument `Hypot` overload live alongsi
them. Exponential functions live in `DoubleDouble.ExponentialFunctions.cs`: them. Exponential functions live in `DoubleDouble.ExponentialFunctions.cs`:
`Exp`, `Exp2`, `Exp10`, `ExpM1`, `Exp2M1`, and `Exp10M1` complete `Exp`, `Exp2`, `Exp10`, `ExpM1`, `Exp2M1`, and `Exp10M1` complete
`IExponentialFunctions<DoubleDouble>`, including cancellation-safe overrides of `IExponentialFunctions<DoubleDouble>`, including cancellation-safe overrides of
the interface's default minus-one methods. `DDMath` exposes thin forwarding the interface's default minus-one methods. Natural logarithm and power kernels
wrappers for roots and exponentials; logarithm, power, and reciprocal kernels live in `DoubleDouble.Log.cs` and `DoubleDouble.Pow.cs`; reciprocal lives in
remain in `DDMath`. `DoubleDouble.NumberBase.cs`. Arithmetic implementations live in
`DoubleDouble.Arithmetic.cs`, with forwarding operators in `DoubleDouble.Operators.cs`.
The `DDMath` static class provides: All mathematical functions are static members of `DoubleDouble`. The former
`DDMath` facade has been removed: migrate `DDMath.Function(...)` calls to
`DoubleDouble.Function(...)`. This is a breaking API change; no compatibility
facade is provided.
`DoubleDoubleExtensions`, in the same `Just.PreciseMath` namespace, also provides
instance-style calls for `Abs`, `Square`, `Reciprocal`, `Sqrt`, `InvSqrt`, `Cbrt`,
all six exponential functions, `Log`, `RootN(int)`, and `Pow` with `int`, `double`,
or `DoubleDouble` exponents. These extensions return a result without modifying
the receiver and forward directly to the corresponding static method, preserving
its numerical contract and component bits. Arithmetic remains available through
operators and named static methods; `Hypot` remains a static multi-operand function.
```csharp
using Just.PreciseMath;
DoubleDouble value = new(-4.0);
DoubleDouble magnitude = value.Abs();
DoubleDouble inverse = value.Abs().Square().Sqrt().Reciprocal();
DoubleDouble power = magnitude.Pow(3);
DoubleDouble root = magnitude.RootN(2);
// value is still -4; inverse is 0.25.
```
`DoubleDouble` provides:
- `Abs(DoubleDouble)`: preserves both components, maps either signed zero to - `Abs(DoubleDouble)`: preserves both components, maps either signed zero to
positive zero and either infinity to positive infinity, and returns canonical NaN. positive zero and either infinity to positive infinity, and returns canonical NaN.
It shares the existing `DoubleDouble.Abs` implementation. - `Square(DoubleDouble)`: uses a specialized double-double squaring kernel.
- `Square(DoubleDouble)`: forwards to the specialized `DoubleDouble.Square` kernel.
It combines the equal cross terms with FMA, normalizes, then incorporates the It combines the equal cross terms with FMA, normalizes, then incorporates the
low-component square without losing selected exactly representable tails. low-component square without losing selected exactly representable tails.
Either signed zero maps to positive zero, either infinity to positive infinity, Either signed zero maps to positive zero, either infinity to positive infinity,
@@ -216,17 +243,17 @@ DoubleDouble cubeRoot = DoubleDouble.Cbrt(new DoubleDouble(-8.0));
DoubleDouble distance = DoubleDouble.Hypot(new DoubleDouble(3.0), new DoubleDouble(4.0)); DoubleDouble distance = DoubleDouble.Hypot(new DoubleDouble(3.0), new DoubleDouble(4.0));
DoubleDouble norm3D = DoubleDouble.Hypot(new DoubleDouble(2.0), new DoubleDouble(3.0), new DoubleDouble(6.0)); DoubleDouble norm3D = DoubleDouble.Hypot(new DoubleDouble(2.0), new DoubleDouble(3.0), new DoubleDouble(6.0));
DoubleDouble fifthRoot = DoubleDouble.RootN(new DoubleDouble(2.0), 5); DoubleDouble fifthRoot = DoubleDouble.RootN(new DoubleDouble(2.0), 5);
DoubleDouble reciprocal = DDMath.Reciprocal(new DoubleDouble(3.0)); DoubleDouble reciprocal = DoubleDouble.Reciprocal(new DoubleDouble(3.0));
DoubleDouble magnitude = DDMath.Abs(-root); DoubleDouble magnitude = DoubleDouble.Abs(-root);
DoubleDouble smallPower = DDMath.Pow(new DoubleDouble(2.0), -1024); DoubleDouble smallPower = DoubleDouble.Pow(new DoubleDouble(2.0), -1024);
DoubleDouble exponential = DoubleDouble.Exp(new DoubleDouble(1.0)); DoubleDouble exponential = DoubleDouble.Exp(new DoubleDouble(1.0));
DoubleDouble binaryPower = DoubleDouble.Exp2(new DoubleDouble(-1000.0)); DoubleDouble binaryPower = DoubleDouble.Exp2(new DoubleDouble(-1000.0));
DoubleDouble decimalPower = DoubleDouble.Exp10(new DoubleDouble(0.5)); DoubleDouble decimalPower = DoubleDouble.Exp10(new DoubleDouble(0.5));
DoubleDouble tinyChange = DoubleDouble.ExpM1(new DoubleDouble(1e-30)); DoubleDouble tinyChange = DoubleDouble.ExpM1(new DoubleDouble(1e-30));
DoubleDouble logarithm = DDMath.Log(new DoubleDouble(10.0)); DoubleDouble logarithm = DoubleDouble.Log(new DoubleDouble(10.0));
DoubleDouble fractionalPower = DDMath.Pow(new DoubleDouble(2.0), 0.5); DoubleDouble fractionalPower = DoubleDouble.Pow(new DoubleDouble(2.0), 0.5);
DoubleDouble preciseExponent = DoubleDouble.FromComponents(0.5, 1e-30); DoubleDouble preciseExponent = DoubleDouble.FromComponents(0.5, 1e-30);
DoubleDouble precisePower = DDMath.Pow(new DoubleDouble(2.0), preciseExponent); DoubleDouble precisePower = DoubleDouble.Pow(new DoubleDouble(2.0), preciseExponent);
``` ```
Reciprocal tests check bitwise equivalence with `1.0 / value` and independently Reciprocal tests check bitwise equivalence with `1.0 / value` and independently
@@ -254,7 +281,7 @@ at 450 and 650 digits, including `int.MinValue` and `int.MaxValue`.
`RootN(x, -1)` retains the reciprocal contract, including its subnormal error floor. `RootN(x, -1)` retains the reciprocal contract, including its subnormal error floor.
Hypotenuse tests use exact squared inequalities with `2^-100` relative error plus Hypotenuse tests use exact squared inequalities with `2^-100` relative error plus
one minimum binary64 subnormal. Tests also check exact binary cases, sparse lows, one minimum binary64 subnormal. Tests also check exact binary cases, sparse lows,
special-value signs, facade component-bit agreement, and constrained generic dispatch. special-value signs and constrained generic dispatch.
The 2D and 3D hypotenuse overflow tests bracket the exact threshold, including scaled The 2D and 3D hypotenuse overflow tests bracket the exact threshold, including scaled
55/48/73 and 1/2/2/3 Pythagorean cases at the midpoint. Three-coordinate tests include 55/48/73 and 1/2/2/3 Pythagorean cases at the midpoint. Three-coordinate tests include
all signs and permutations of boundary cases and sparse corrections, including all signs and permutations of boundary cases and sparse corrections, including
@@ -275,9 +302,9 @@ check `2^-100` relative error plus one minimum binary64 subnormal, with a separa
bounded reference-rounding allowance. The base-two/base-ten and minus-one fixtures 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, are checked at 450/650 digits, including adjacent lows at true range boundaries,
series transitions, and tiny inputs. Additional component checks pin sparse series transitions, and tiny inputs. Additional component checks pin sparse
base-two corrections; constrained generic calls and facade calls match direct base-two corrections; constrained generic calls match direct calls bit-for-bit.
calls bit-for-bit. Integer-power tests use exact rational Integer-power tests use exact rational references and high-precision fixtures for
references and high-precision fixtures for large exponents; their tested absolute large exponents; their tested absolute
error bound is `|exact result| * (|exponent| + 1) * 2^-100 + 2^-1074`, not uniform error bound is `|exact result| * (|exponent| + 1) * 2^-100 + 2^-1074`, not uniform
relative accuracy independent of the exponent. These functions are approximate; relative accuracy independent of the exponent. These functions are approximate;
precision decreases near underflow and approximation can affect results extremely precision decreases near underflow and approximation can affect results extremely
@@ -402,8 +429,8 @@ bool success = DoubleDouble.TryParse("1.25e-2".AsSpan(), CultureInfo.InvariantCu
## Deferred scope ## Deferred scope
Natural `DDMath.Log`, the complete exponential family on `DoubleDouble`, and all Natural `DoubleDouble.Log`, the complete exponential family on `DoubleDouble`, and all
three `DDMath.Pow` overloads are implemented. three `DoubleDouble.Pow` overloads are implemented.
Logarithms in other bases, generic-math interfaces beyond `ISignedNumber`, Logarithms in other bases, generic-math interfaces beyond `ISignedNumber`,
`IFloatingPointConstants`, `IRootFunctions`, and `IExponentialFunctions`, additional text formats/shortest- `IFloatingPointConstants`, `IRootFunctions`, and `IExponentialFunctions`, additional text formats/shortest-
round-trip formatting, and non-arithmetic performance benchmarks remain deferred. round-trip formatting, and non-arithmetic performance benchmarks remain deferred.