diff --git a/0-source/Just.PreciseMath/DDMath.Exp.cs b/0-source/Just.PreciseMath/DDMath.Exp.cs deleted file mode 100644 index e06c287..0000000 --- a/0-source/Just.PreciseMath/DDMath.Exp.cs +++ /dev/null @@ -1,46 +0,0 @@ -namespace Just.PreciseMath; - -public static partial class DDMath -{ - /// - [Pure] - public static DoubleDouble Exp(DoubleDouble value) - { - return DoubleDouble.Exp(value); - } - - /// - [Pure] - public static DoubleDouble Exp2(DoubleDouble value) - { - return DoubleDouble.Exp2(value); - } - - /// - [Pure] - public static DoubleDouble Exp10(DoubleDouble value) - { - return DoubleDouble.Exp10(value); - } - - /// - [Pure] - public static DoubleDouble ExpM1(DoubleDouble value) - { - return DoubleDouble.ExpM1(value); - } - - /// - [Pure] - public static DoubleDouble Exp2M1(DoubleDouble value) - { - return DoubleDouble.Exp2M1(value); - } - - /// - [Pure] - public static DoubleDouble Exp10M1(DoubleDouble value) - { - return DoubleDouble.Exp10M1(value); - } -} diff --git a/0-source/Just.PreciseMath/DDMath.cs b/0-source/Just.PreciseMath/DDMath.cs deleted file mode 100644 index ce2a445..0000000 --- a/0-source/Just.PreciseMath/DDMath.cs +++ /dev/null @@ -1,70 +0,0 @@ -namespace Just.PreciseMath; - -/// Provides mathematical functions for normalized double-double values. -/// Results use fixed-size double-double precision, not arbitrary precision. -public static partial class DDMath -{ - /// Returns the absolute value without discarding the low component. - /// Both signs of zero become positive zero; NaN is canonicalized and either infinity becomes positive infinity. - [Pure] - public static DoubleDouble Abs(DoubleDouble value) - { - return DoubleDouble.Abs(value); - } - - /// - [Pure] - public static DoubleDouble Square(DoubleDouble value) - { - return DoubleDouble.Square(value); - } - - /// - [Pure] - public static DoubleDouble Reciprocal(DoubleDouble value) - { - return DoubleDouble.Reciprocal(value); - } - - /// - [Pure] - public static DoubleDouble Cbrt(DoubleDouble value) - { - return DoubleDouble.Cbrt(value); - } - - /// - [Pure] - public static DoubleDouble Hypot(DoubleDouble x, DoubleDouble y) - { - return DoubleDouble.Hypot(x, y); - } - - /// - [Pure] - public static DoubleDouble Hypot(DoubleDouble x, DoubleDouble y, DoubleDouble z) - { - return DoubleDouble.Hypot(x, y, z); - } - - /// - [Pure] - public static DoubleDouble RootN(DoubleDouble value, int n) - { - return DoubleDouble.RootN(value, n); - } - - /// - [Pure] - public static DoubleDouble Sqrt(DoubleDouble value) - { - return DoubleDouble.Sqrt(value); - } - - /// - [Pure] - public static DoubleDouble InvSqrt(DoubleDouble value) - { - return DoubleDouble.InvSqrt(value); - } -} diff --git a/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs b/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs index 0c7ac16..785322b 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs @@ -44,13 +44,13 @@ public readonly partial struct DoubleDouble : } /// Returns the operand unchanged. - public static DoubleDouble operator +(DoubleDouble value) + public static DoubleDouble Plus(DoubleDouble value) { return value; } /// Negates the value, including the high zero's sign, retaining canonical NaN and zero residuals. - public static DoubleDouble operator -(DoubleDouble value) + public static DoubleDouble Negate(DoubleDouble value) { // Negation preserves normalization; only NaN and zero residuals need canonicalization. return new DoubleDouble(double.IsNaN(value._high) ? double.NaN : -value._high, @@ -58,7 +58,7 @@ public readonly partial struct DoubleDouble : } /// Adds normalized expansions, retaining low-sum residuals under cancellation. - 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)) { @@ -73,7 +73,7 @@ public readonly partial struct DoubleDouble : } /// Subtracts normalized expansions. - 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)) { @@ -95,7 +95,7 @@ public readonly partial struct DoubleDouble : /// Swapping operands preserves both result components bit-for-bit. /// Results are approximate double-double values, not universally correctly rounded. /// - 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) { @@ -127,7 +127,7 @@ public readonly partial struct DoubleDouble : /// Divides expansions using a quotient estimate and two residual corrections. /// Zero and nonfinite operands follow binary64 rules; precision decreases near underflow. - 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) { @@ -159,35 +159,35 @@ public readonly partial struct DoubleDouble : /// Applies the expansion operation without discarding the low component. /// Preserves the component bits of addition with the equivalent zero-low scalar expansion. - public static DoubleDouble operator +(DoubleDouble left, double right) + public static DoubleDouble Add(DoubleDouble left, double right) { return PreciseMathHelper.AddScalar(left._high, left._low, right); } /// Applies the expansion operation without discarding the low component. /// Preserves the component bits of addition with the equivalent zero-low scalar expansion. - public static DoubleDouble operator +(double left, DoubleDouble right) + public static DoubleDouble Add(double left, DoubleDouble right) { - return right + left; + return Add(right, left); } /// Applies the expansion operation without discarding the low component. /// Preserves the component bits of subtraction with the equivalent zero-low scalar expansion. - public static DoubleDouble operator -(DoubleDouble left, double right) + public static DoubleDouble Subtract(DoubleDouble left, double right) { return PreciseMathHelper.AddScalar(left._high, left._low, -right); } /// Applies the expansion operation without discarding the low component. /// Preserves the component bits of subtraction with the equivalent zero-low scalar expansion. - 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. return PreciseMathHelper.AddScalar(-right._high, -right._low, left); } /// Applies the expansion operation without discarding the low component. - 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) { @@ -207,14 +207,14 @@ public readonly partial struct DoubleDouble : } /// Applies the expansion operation without discarding the low component. - public static DoubleDouble operator *(double left, DoubleDouble right) + public static DoubleDouble Multiply(double left, DoubleDouble right) { - return right * left; + return Multiply(right, left); } /// Divides by a scalar using a quotient estimate and two residual corrections. /// Preserves the component bits of division by the equivalent zero-low expansion. - 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) { @@ -246,7 +246,7 @@ public readonly partial struct DoubleDouble : /// Divides a scalar using a quotient estimate and two residual corrections. /// Preserves the component bits of division with the equivalent zero-low numerator. - 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) { diff --git a/0-source/Just.PreciseMath/DDMath.Log.cs b/0-source/Just.PreciseMath/DoubleDouble.Log.cs similarity index 88% rename from 0-source/Just.PreciseMath/DDMath.Log.cs rename to 0-source/Just.PreciseMath/DoubleDouble.Log.cs index 998d80c..9dc83f4 100644 --- a/0-source/Just.PreciseMath/DDMath.Log.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.Log.cs @@ -1,6 +1,6 @@ namespace Just.PreciseMath; -public static partial class DDMath +public readonly partial struct DoubleDouble { /// Returns the natural logarithm of a double-double value. /// @@ -17,7 +17,7 @@ public static partial class DDMath { if (double.IsNaN(value.High) || value.High < 0.0) { - return DoubleDouble.NaN; + return NaN; } if (value.High == 0.0) { @@ -27,9 +27,9 @@ public static partial class DDMath { return value; } - if (value == DoubleDouble.One) + if (value == One) { - return DoubleDouble.Zero; + return Zero; } return LogPositiveFinite(value); } @@ -49,9 +49,9 @@ public static partial class DDMath } 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)); - if (mantissa.High > DoubleDouble.Sqrt2.High) + if (mantissa.High > Sqrt2.High) { mantissa *= 0.5; ++exponent; @@ -72,6 +72,6 @@ public static partial class DDMath } // log(m) = 2*atanh(t). The omitted series tail is bounded by // 2*|t|^51/(51*(1-t^2)); rounding, not truncation, dominates. - return (sum * 2.0) + (DoubleDouble.Ln2 * exponent); + return (sum * 2.0) + (Ln2 * exponent); } } diff --git a/0-source/Just.PreciseMath/DoubleDouble.NumberBase.cs b/0-source/Just.PreciseMath/DoubleDouble.NumberBase.cs index 58feb7b..b0086d3 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.NumberBase.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.NumberBase.cs @@ -86,18 +86,6 @@ public readonly partial struct DoubleDouble return IsNaN(y) ? x : MinMagnitude(x, y); } - /// Adds one using double-double arithmetic, including its overflow and nonfinite behavior. - public static DoubleDouble operator ++(DoubleDouble value) - { - return value + 1.0; - } - - /// Subtracts one using double-double arithmetic, including its overflow and nonfinite behavior. - public static DoubleDouble operator --(DoubleDouble value) - { - return value - 1.0; - } - /// Tests normalization and canonical NaN and zero-residual bit patterns. public static bool IsCanonical(DoubleDouble value) { diff --git a/0-source/Just.PreciseMath/DoubleDouble.Operators.cs b/0-source/Just.PreciseMath/DoubleDouble.Operators.cs new file mode 100644 index 0000000..664c4f5 --- /dev/null +++ b/0-source/Just.PreciseMath/DoubleDouble.Operators.cs @@ -0,0 +1,100 @@ +namespace Just.PreciseMath; + +public readonly partial struct DoubleDouble +{ + /// + public static DoubleDouble operator +(DoubleDouble value) + { + return Plus(value); + } + + /// + public static DoubleDouble operator -(DoubleDouble value) + { + return Negate(value); + } + + /// + public static DoubleDouble operator +(DoubleDouble left, DoubleDouble right) + { + return Add(left, right); + } + + /// + public static DoubleDouble operator -(DoubleDouble left, DoubleDouble right) + { + return Subtract(left, right); + } + + /// + public static DoubleDouble operator *(DoubleDouble left, DoubleDouble right) + { + return Multiply(left, right); + } + + /// + public static DoubleDouble operator /(DoubleDouble left, DoubleDouble right) + { + return Divide(left, right); + } + + /// + public static DoubleDouble operator +(DoubleDouble left, double right) + { + return Add(left, right); + } + + /// + public static DoubleDouble operator +(double left, DoubleDouble right) + { + return Add(left, right); + } + + /// + public static DoubleDouble operator -(DoubleDouble left, double right) + { + return Subtract(left, right); + } + + /// + public static DoubleDouble operator -(double left, DoubleDouble right) + { + return Subtract(left, right); + } + + /// + public static DoubleDouble operator *(DoubleDouble left, double right) + { + return Multiply(left, right); + } + + /// + public static DoubleDouble operator *(double left, DoubleDouble right) + { + return Multiply(left, right); + } + + /// + public static DoubleDouble operator /(DoubleDouble left, double right) + { + return Divide(left, right); + } + + /// + public static DoubleDouble operator /(double left, DoubleDouble right) + { + return Divide(left, right); + } + + /// Adds one using double-double arithmetic, including its overflow and nonfinite behavior. + public static DoubleDouble operator ++(DoubleDouble value) + { + return Add(value, 1.0); + } + + /// Subtracts one using double-double arithmetic, including its overflow and nonfinite behavior. + public static DoubleDouble operator --(DoubleDouble value) + { + return Subtract(value, 1.0); + } +} diff --git a/0-source/Just.PreciseMath/DDMath.Pow.cs b/0-source/Just.PreciseMath/DoubleDouble.Pow.cs similarity index 83% rename from 0-source/Just.PreciseMath/DDMath.Pow.cs rename to 0-source/Just.PreciseMath/DoubleDouble.Pow.cs index 7ddfbdf..e888ac6 100644 --- a/0-source/Just.PreciseMath/DDMath.Pow.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.Pow.cs @@ -1,6 +1,6 @@ namespace Just.PreciseMath; -public static partial class DDMath +public readonly partial struct DoubleDouble { /// Raises a double-double value to an integer power. /// @@ -21,11 +21,11 @@ public static partial class DDMath { if (exponent == 0) { - return DoubleDouble.One; + return One; } if (double.IsNaN(value.High)) { - return DoubleDouble.NaN; + return NaN; } if (exponent == 1) { @@ -43,10 +43,10 @@ public static partial class DDMath // Widen BEFORE negation: abs(int.MinValue) is not an int. long remaining = exponent < 0 ? -(long)exponent : exponent; long factorExponent = Math.ILogB(value.High); - DoubleDouble factor = DoubleDouble.FromComponents( + DoubleDouble factor = FromComponents( Math.ScaleB(Math.Abs(value.High), -(int)factorExponent), Math.ScaleB(value.High < 0.0 ? -value.Low : value.Low, -(int)factorExponent)); - DoubleDouble result = DoubleDouble.One; + DoubleDouble result = One; long resultExponent = 0; while (remaining != 0) { @@ -105,35 +105,35 @@ public static partial class DDMath [Pure] 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); - if (DoubleDouble.IsInfinity(exponent)) + DoubleDouble magnitude = Abs(value); + if (IsInfinity(exponent)) { - if (magnitude == DoubleDouble.One) + if (magnitude == One) { - return DoubleDouble.One; + return One; } - bool grows = (magnitude > DoubleDouble.One) == (exponent.High > 0.0); - return grows ? new DoubleDouble(double.PositiveInfinity) : DoubleDouble.Zero; + bool grows = (magnitude > One) == (exponent.High > 0.0); + 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) { // A normalized integer in this range has a zero low component. // Reuse the integer algorithm, including exact identities and ties. return Pow(value, (int)exponent.High); } - bool negative = double.IsNegative(value.High) && DoubleDouble.IsOddInteger(exponent); - if (value.High == 0.0 || DoubleDouble.IsInfinity(value)) + bool negative = double.IsNegative(value.High) && IsOddInteger(exponent); + if (value.High == 0.0 || IsInfinity(value)) { bool infinite = (value.High == 0.0) == (exponent.High < 0.0); double result = infinite ? double.PositiveInfinity : 0.0; @@ -141,11 +141,11 @@ public static partial class DDMath } 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)); return negative ? -resultMagnitude : resultMagnitude; @@ -158,7 +158,7 @@ public static partial class DDMath // int.MinValue times any finite input's binary64 exponent. int adjustment = Math.ILogB(value.High); exponent += adjustment; - return DoubleDouble.FromComponents(Math.ScaleB(value.High, -adjustment), + return FromComponents(Math.ScaleB(value.High, -adjustment), Math.ScaleB(value.Low, -adjustment)); } } diff --git a/0-source/Just.PreciseMath/DoubleDouble.RootFunctions.cs b/0-source/Just.PreciseMath/DoubleDouble.RootFunctions.cs index 606594c..abf19ba 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.RootFunctions.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.RootFunctions.cs @@ -200,7 +200,7 @@ public readonly partial struct DoubleDouble : IRootFunctions case 1: return value; case -1: - return DDMath.Reciprocal(value); + return Reciprocal(value); case 2: return Sqrt(value); case -2: diff --git a/0-source/Just.PreciseMath/DoubleDoubleExtensions.cs b/0-source/Just.PreciseMath/DoubleDoubleExtensions.cs new file mode 100644 index 0000000..2a52c69 --- /dev/null +++ b/0-source/Just.PreciseMath/DoubleDoubleExtensions.cs @@ -0,0 +1,128 @@ +namespace Just.PreciseMath; + +/// Provides instance-style calls to the static DoubleDouble mathematical functions. +/// +/// Every extension returns a new result without modifying its receiver. Extensions +/// forward directly to DoubleDouble, retaining its numerical and special-value contracts. +/// +public static class DoubleDoubleExtensions +{ + /// + [Pure] + public static DoubleDouble Abs(this DoubleDouble value) + { + return DoubleDouble.Abs(value); + } + + /// + [Pure] + public static DoubleDouble Square(this DoubleDouble value) + { + return DoubleDouble.Square(value); + } + + /// + [Pure] + public static DoubleDouble Reciprocal(this DoubleDouble value) + { + return DoubleDouble.Reciprocal(value); + } + + /// + [Pure] + public static DoubleDouble Sqrt(this DoubleDouble value) + { + return DoubleDouble.Sqrt(value); + } + + /// + [Pure] + public static DoubleDouble InvSqrt(this DoubleDouble value) + { + return DoubleDouble.InvSqrt(value); + } + + /// + [Pure] + public static DoubleDouble Cbrt(this DoubleDouble value) + { + return DoubleDouble.Cbrt(value); + } + + /// + [Pure] + public static DoubleDouble Exp(this DoubleDouble value) + { + return DoubleDouble.Exp(value); + } + + /// + [Pure] + public static DoubleDouble Exp2(this DoubleDouble value) + { + return DoubleDouble.Exp2(value); + } + + /// + [Pure] + public static DoubleDouble Exp10(this DoubleDouble value) + { + return DoubleDouble.Exp10(value); + } + + /// + [Pure] + public static DoubleDouble ExpM1(this DoubleDouble value) + { + return DoubleDouble.ExpM1(value); + } + + /// + [Pure] + public static DoubleDouble Exp2M1(this DoubleDouble value) + { + return DoubleDouble.Exp2M1(value); + } + + /// + [Pure] + public static DoubleDouble Exp10M1(this DoubleDouble value) + { + return DoubleDouble.Exp10M1(value); + } + + /// + [Pure] + public static DoubleDouble Log(this DoubleDouble value) + { + return DoubleDouble.Log(value); + } + + /// + [Pure] + public static DoubleDouble RootN(this DoubleDouble value, int n) + { + return DoubleDouble.RootN(value, n); + } + + /// + [Pure] + public static DoubleDouble Pow(this DoubleDouble value, int exponent) + { + return DoubleDouble.Pow(value, exponent); + } + + /// + [Pure] + public static DoubleDouble Pow(this DoubleDouble value, double exponent) + { + return DoubleDouble.Pow(value, exponent); + } + + /// + [Pure] + public static DoubleDouble Pow(this DoubleDouble value, DoubleDouble exponent) + { + return DoubleDouble.Pow(value, exponent); + } +} diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs index 0a26ad7..2c4393c 100644 --- a/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleArithmeticTests.cs @@ -2,6 +2,55 @@ namespace Just.PreciseMath.Tests; 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] [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)] @@ -1050,7 +1099,7 @@ public class DoubleDoubleArithmeticTests 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(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 }) { CheckBits(DoubleDouble.Square(new DoubleDouble(value)), value * value); - CheckBits(DDMath.Square(new DoubleDouble(value)), value * value); } for (int exponent = -1074; exponent <= 1023; exponent++) { @@ -1139,7 +1187,6 @@ public class DoubleDoubleArithmeticTests DoubleDouble.IsZero(actual).ShouldBe((numerator << 1) <= denominator); } CheckBoundary(DoubleDouble.Square(-value), actual.High, actual.Low); - CheckBoundary(DDMath.Square(value), actual.High, actual.Low); } private static DoubleDouble[] MultiplicationCommutativityCases() diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleCbrtTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleCbrtTests.cs index f8d999f..cd10de7 100644 --- a/1-tests/Just.PreciseMath.Tests/DoubleDoubleCbrtTests.cs +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleCbrtTests.cs @@ -36,7 +36,6 @@ public class DoubleDoubleCbrtTests { AssertBound(input, DoubleDouble.Cbrt(input)); AssertBound(-input, DoubleDouble.Cbrt(-input)); - AssertBits(DoubleDouble.Cbrt(input), DDMath.Cbrt(input)); } } } diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleExponentialFunctionsTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleExponentialFunctionsTests.cs index 7e90f0e..e554aa8 100644 --- a/1-tests/Just.PreciseMath.Tests/DoubleDoubleExponentialFunctionsTests.cs +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleExponentialFunctionsTests.cs @@ -20,7 +20,6 @@ public class DoubleDoubleExponentialFunctionsTests DoubleDouble actual = Evaluate(operation, input); DoubleDouble.IsCanonical(actual).ShouldBeTrue(); AssertBits(actual, EvaluateGeneric(operation, input)); - AssertBits(actual, EvaluateFacade(operation, input)); if (overflow) { AssertBits(DoubleDouble.PositiveInfinity, actual); @@ -57,7 +56,7 @@ public class DoubleDoubleExponentialFunctionsTests } [Fact] - public void NaturalExponentialIsAvailableOnTheTypeAndMatchesFacadeBits() + public void NaturalExponentialMatchesGenericDispatchBits() { DoubleDouble[] values = [DoubleDouble.Zero, DoubleDouble.NegativeZero, 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))]; 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 expected = DoubleDouble.FromComponents(Math.ScaleB(1.0, exponent), sign * expectedLow); AssertBits(expected, DoubleDouble.Exp2(input)); - AssertBits(expected, DDMath.Exp2(input)); } } @@ -128,7 +126,7 @@ public class DoubleDoubleExponentialFunctionsTests [InlineData("ExpM1")] [InlineData("Exp2M1")] [InlineData("Exp10M1")] - public void CompleteContractHasGenericDispatchAndThinFacade(string operation) + public void CompleteContractHasGenericDispatch(string operation) { bool minusOne = operation.EndsWith("M1", StringComparison.Ordinal); DoubleDouble[] values = [DoubleDouble.Zero, DoubleDouble.NegativeZero, @@ -141,7 +139,6 @@ public class DoubleDoubleExponentialFunctionsTests DoubleDouble actual = Evaluate(operation, value); DoubleDouble.IsCanonical(actual).ShouldBeTrue(); AssertBits(actual, EvaluateGeneric(operation, value)); - AssertBits(actual, EvaluateFacade(operation, value)); if (DoubleDouble.IsZero(value)) { AssertBits(minusOne ? value : DoubleDouble.One, actual); @@ -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) { long bits = BitConverter.DoubleToInt64Bits(value); diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleExtensionsTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleExtensionsTests.cs new file mode 100644 index 0000000..c3fcb47 --- /dev/null +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleExtensionsTests.cs @@ -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 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)); + } +} diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleHypotTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleHypotTests.cs index d228c6b..e33efc1 100644 --- a/1-tests/Just.PreciseMath.Tests/DoubleDoubleHypotTests.cs +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleHypotTests.cs @@ -52,7 +52,6 @@ public class DoubleDoubleHypotTests DoubleDouble result = DoubleDouble.Hypot(x, y); AssertBound(x, y, result); 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); } - 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); AssertBits(result, DoubleDouble.Hypot(-z, x, -y)); 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(c, a, b)); 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); } AssertBits(result, DoubleDouble.Hypot(second, first)); - AssertBits(result, DDMath.Hypot(first, second)); } } } diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleRootFunctionsTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleRootFunctionsTests.cs index 54bbaee..e6d8e51 100644 --- a/1-tests/Just.PreciseMath.Tests/DoubleDoubleRootFunctionsTests.cs +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleRootFunctionsTests.cs @@ -3,7 +3,7 @@ namespace Just.PreciseMath.Tests; public class DoubleDoubleRootFunctionsTests { [Fact] - public void ExistingRootsAreAvailableOnTheTypeAndMatchFacadeBits() + public void RootsAreCanonicalAndSqrtMatchesGenericDispatchBits() { DoubleDouble[] values = [DoubleDouble.Zero, DoubleDouble.NegativeZero, DoubleDouble.NaN, DoubleDouble.PositiveInfinity, DoubleDouble.NegativeInfinity, @@ -11,11 +11,16 @@ public class DoubleDoubleRootFunctionsTests DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -1000))]; foreach (DoubleDouble value in values) { - AssertBits(DDMath.Sqrt(value), DoubleDouble.Sqrt(value)); - AssertBits(DDMath.InvSqrt(value), DoubleDouble.InvSqrt(value)); + AssertBits(DoubleDouble.Sqrt(value), GenericSqrt(value)); + DoubleDouble.IsCanonical(DoubleDouble.InvSqrt(value)).ShouldBeTrue(); } } + private static T GenericSqrt(T value) where T : IRootFunctions + { + return T.Sqrt(value); + } + private static void AssertBits(DoubleDouble expected, DoubleDouble actual) { BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(BitConverter.DoubleToInt64Bits(expected.High)); diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleRootNTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleRootNTests.cs index 1a4d8ed..0bf40f3 100644 --- a/1-tests/Just.PreciseMath.Tests/DoubleDoubleRootNTests.cs +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleRootNTests.cs @@ -68,7 +68,7 @@ public class DoubleDoubleRootNTests { 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.InvSqrt(value), DoubleDouble.RootN(value, -2)); AssertBits(DoubleDouble.Cbrt(value), DoubleDouble.RootN(value, 3)); @@ -153,7 +153,7 @@ public class DoubleDoubleRootNTests } [Fact] - public void GenericContractAndFacadeDispatchToTypeMembers() + public void GenericContractDispatchesToTypeMembers() { DoubleDouble[] values = [DoubleDouble.Zero, DoubleDouble.NegativeZero, DoubleDouble.NaN, DoubleDouble.PositiveInfinity, DoubleDouble.NegativeInfinity, new(-2.0), new(double.Epsilon), @@ -167,7 +167,6 @@ public class DoubleDoubleRootNTests { DoubleDouble result = DoubleDouble.RootN(value, degree); AssertBits(result, GenericRootN(value, degree)); - AssertBits(result, DDMath.RootN(value, degree)); } } } diff --git a/1-tests/Just.PreciseMath.Tests/PreciseMathExpTests.cs b/1-tests/Just.PreciseMath.Tests/PreciseMathExpTests.cs index 9b9a3dd..db73c66 100644 --- a/1-tests/Just.PreciseMath.Tests/PreciseMathExpTests.cs +++ b/1-tests/Just.PreciseMath.Tests/PreciseMathExpTests.cs @@ -16,7 +16,7 @@ public class PreciseMathExpTests // Inputs are exact stored binary64 sums, not rounded irrational constants. DoubleDouble input = DoubleDouble.FromComponents(high, 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(); if (overflow) { @@ -60,7 +60,7 @@ public class PreciseMathExpTests { // Independently round Decimal.exp(exact binary64 input)/2^-1074 to an // 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.Low).ShouldBe(0L); } @@ -70,7 +70,7 @@ public class PreciseMathExpTests [InlineData(709.781)] public void LegacyFiniteResultsDoNotOverflowPrematurely(double input) { - DoubleDouble actual = DDMath.Exp(new DoubleDouble(input)); + DoubleDouble actual = DoubleDouble.Exp(new DoubleDouble(input)); DoubleDouble.IsFinite(actual).ShouldBeTrue(); actual.Low.ShouldNotBe(0.0); // Component-aware accuracy for these exact inputs is pinned by ReferenceCases. @@ -83,7 +83,7 @@ public class PreciseMathExpTests double.MaxValue, double.MinValue]; foreach (double input in inputs) { - DoubleDouble actual = DDMath.Exp(new DoubleDouble(input)); + DoubleDouble actual = DoubleDouble.Exp(new DoubleDouble(input)); double expected = Math.Exp(input); BitConverter.DoubleToInt64Bits(actual.High).ShouldBe( BitConverter.DoubleToInt64Bits(double.IsNaN(expected) ? double.NaN : expected)); diff --git a/1-tests/Just.PreciseMath.Tests/PreciseMathInvSqrtTests.cs b/1-tests/Just.PreciseMath.Tests/PreciseMathInvSqrtTests.cs index 0f26469..021011b 100644 --- a/1-tests/Just.PreciseMath.Tests/PreciseMathInvSqrtTests.cs +++ b/1-tests/Just.PreciseMath.Tests/PreciseMathInvSqrtTests.cs @@ -11,7 +11,7 @@ public class PreciseMathInvSqrtTests public void IrrationalInverseRootsRetainMoreThanBinary64Precision(double value) { DoubleDouble input = new(value); - DoubleDouble actual = DDMath.InvSqrt(input); + DoubleDouble actual = DoubleDouble.InvSqrt(input); actual.Low.ShouldNotBe(0.0); AssertInvSqrtBound(input, actual); } @@ -23,7 +23,7 @@ public class PreciseMathInvSqrtTests double.NaN, -1.0, -double.Epsilon, double.MinValue]; 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. double expected = 1.0 / Math.Sqrt(value); BitConverter.DoubleToInt64Bits(actual.High).ShouldBe( @@ -33,7 +33,7 @@ public class PreciseMathInvSqrtTests } DoubleDouble negative = DoubleDouble.FromComponents(-1.0, Math.ScaleB(1.0, -54)); - DoubleDouble.IsNaN(DDMath.InvSqrt(negative)).ShouldBeTrue(); + DoubleDouble.IsNaN(DoubleDouble.InvSqrt(negative)).ShouldBeTrue(); } [Fact] @@ -42,7 +42,7 @@ public class PreciseMathInvSqrtTests // 1/sqrt(2^(2k)) = 2^-k exactly, including the minimum subnormal input. 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))); 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) }) { DoubleDouble input = DoubleDouble.FromComponents(1.0, low); - DoubleDouble actual = DDMath.InvSqrt(input); + DoubleDouble actual = DoubleDouble.InvSqrt(input); actual.High.ShouldBe(1.0); Math.Sign(actual.Low).ShouldBe(-Math.Sign(low)); AssertInvSqrtBound(input, actual); @@ -78,7 +78,7 @@ public class PreciseMathInvSqrtTests if (boundary > 0.0) { 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 }) @@ -86,7 +86,7 @@ public class PreciseMathInvSqrtTests DoubleDouble input = DoubleDouble.FromComponents(high, residual); 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)]; 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); 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 }) { 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.Low.ShouldBe(Math.ScaleB(-sign, exponent - 1)); DoubleDouble.IsCanonical(actual).ShouldBeTrue(); diff --git a/1-tests/Just.PreciseMath.Tests/PreciseMathLogTests.cs b/1-tests/Just.PreciseMath.Tests/PreciseMathLogTests.cs index 94c3308..e5b554c 100644 --- a/1-tests/Just.PreciseMath.Tests/PreciseMathLogTests.cs +++ b/1-tests/Just.PreciseMath.Tests/PreciseMathLogTests.cs @@ -17,10 +17,10 @@ public class PreciseMathLogTests foreach (double value in values) { 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)); - AssertBits(DDMath.Log(negative), double.NaN, 0.0); + AssertBits(DoubleDouble.Log(negative), double.NaN, 0.0); } [Theory] @@ -31,7 +31,7 @@ public class PreciseMathLogTests // 2200 digits and checked with Fraction. No binary64 log supplies the oracle. DoubleDouble input = DoubleDouble.FromComponents(high, 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})"); if (input == DoubleDouble.One) { @@ -51,7 +51,7 @@ public class PreciseMathLogTests for (int exponent = -1074; exponent <= 1023; ++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 // 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(); actual.Low.ShouldNotBe(0.0); DoubleDouble.IsCanonical(actual).ShouldBeTrue(); @@ -75,12 +75,12 @@ public class PreciseMathLogTests { 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); 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 // 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)); } } @@ -89,8 +89,8 @@ public class PreciseMathLogTests public void ArbitraryPairsAreNormalizedByThePublicFactoryBeforeLog() { DoubleDouble positive = DoubleDouble.FromComponents(0.0, 2.0); - AssertReferenceBound(DDMath.Log(positive), LogReferenceData.Ln2, 1, "zero-high public construction"); - AssertBits(DDMath.Log(DoubleDouble.FromComponents(0.0, -2.0)), double.NaN, 0.0); + AssertReferenceBound(DoubleDouble.Log(positive), LogReferenceData.Ln2, 1, "zero-high public construction"); + 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) diff --git a/1-tests/Just.PreciseMath.Tests/PreciseMathPowTests.cs b/1-tests/Just.PreciseMath.Tests/PreciseMathPowTests.cs index 9f76dc6..dd0f1f4 100644 --- a/1-tests/Just.PreciseMath.Tests/PreciseMathPowTests.cs +++ b/1-tests/Just.PreciseMath.Tests/PreciseMathPowTests.cs @@ -6,17 +6,17 @@ public class PreciseMathPowTests public void NegativePowerDoesNotOverflowBeforeReciprocation() { // 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); } [Fact] public void MinimumIntegerExponentDoesNotOverflowItsMagnitude() { - AssertBits(DDMath.Pow(DoubleDouble.One, int.MinValue), 1.0, 0.0); - AssertBits(DDMath.Pow(DoubleDouble.NegativeOne, int.MinValue), 1.0, 0.0); - AssertBits(DDMath.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(DoubleDouble.One, int.MinValue), 1.0, 0.0); + AssertBits(DoubleDouble.Pow(DoubleDouble.NegativeOne, int.MinValue), 1.0, 0.0); + AssertBits(DoubleDouble.Pow(new DoubleDouble(2.0), int.MinValue), 0.0, 0.0); + AssertBits(DoubleDouble.Pow(new DoubleDouble(0.5), int.MinValue), double.PositiveInfinity, 0.0); } [Fact] @@ -25,7 +25,7 @@ public class PreciseMathPowTests foreach (double sign in new[] { -1.0, 1.0 }) { 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 // special-value/sign oracle, not the finite precision oracle. 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); } } @@ -63,7 +63,7 @@ public class PreciseMathPowTests foreach (double sign in new[] { -1.0, 1.0 }) { 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); } } @@ -118,13 +118,13 @@ public class PreciseMathPowTests AssertExactPowerBound(input, -1); if (scale == -500 || scale == 500) { - DoubleDouble actual = DDMath.Pow(input, -2); + DoubleDouble actual = DoubleDouble.Pow(input, -2); actual.Low.ShouldNotBe(0.0); AssertExactPowerBound(input, -2); } 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 }) { 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)); - 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); - 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), 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); DoubleDouble overflowing = DoubleDouble.FromComponents(Math.ScaleB(1.0, scale), 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 }) { 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); int resultSign = sign < 0.0 && (exponent & 1) != 0 ? -1 : 1; // 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)); } - AssertRationalBound(DDMath.Pow(input, exponent), numerator, denominator, exponent, + AssertRationalBound(DoubleDouble.Pow(input, exponent), numerator, denominator, exponent, $"({input.High:R}, {input.Low:R})^{exponent}"); } diff --git a/1-tests/Just.PreciseMath.Tests/PreciseMathRealPowSpecialTests.cs b/1-tests/Just.PreciseMath.Tests/PreciseMathRealPowSpecialTests.cs index 08af8a8..ab99e45 100644 --- a/1-tests/Just.PreciseMath.Tests/PreciseMathRealPowSpecialTests.cs +++ b/1-tests/Just.PreciseMath.Tests/PreciseMathRealPowSpecialTests.cs @@ -10,10 +10,10 @@ public class PreciseMathRealPowSpecialTests DoubleDouble.FromComponents(double.MaxValue, double.Epsilon)]; foreach (DoubleDouble value in values) { - AssertBits(DDMath.Pow(value, 0.0), DoubleDouble.One); - AssertBits(DDMath.Pow(value, new DoubleDouble(-0.0)), DoubleDouble.One); - AssertBits(DDMath.Pow(value, 1.0), value); - AssertBits(DDMath.Pow(value, DoubleDouble.One), value); + AssertBits(DoubleDouble.Pow(value, 0.0), DoubleDouble.One); + AssertBits(DoubleDouble.Pow(value, new DoubleDouble(-0.0)), DoubleDouble.One); + AssertBits(DoubleDouble.Pow(value, 1.0), value); + AssertBits(DoubleDouble.Pow(value, DoubleDouble.One), value); } } @@ -37,8 +37,8 @@ public class PreciseMathRealPowSpecialTests continue; } DoubleDouble input = new(value); - AssertBits(DDMath.Pow(input, exponent), new DoubleDouble(expected)); - AssertBits(DDMath.Pow(input, new DoubleDouble(exponent)), new DoubleDouble(expected)); + AssertBits(DoubleDouble.Pow(input, 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 even = DoubleDouble.FromComponents(direction * high, 2.0); DoubleDouble fractional = DoubleDouble.FromComponents(direction * high, 0.5); - AssertBits(DDMath.Pow(DoubleDouble.NegativeOne, odd), DoubleDouble.NegativeOne); - AssertBits(DDMath.Pow(DoubleDouble.NegativeOne, even), DoubleDouble.One); - AssertBits(DDMath.Pow(DoubleDouble.NegativeOne, fractional), DoubleDouble.NaN); - AssertBits(DDMath.Pow(new DoubleDouble(-2.0), odd), + AssertBits(DoubleDouble.Pow(DoubleDouble.NegativeOne, odd), DoubleDouble.NegativeOne); + AssertBits(DoubleDouble.Pow(DoubleDouble.NegativeOne, even), DoubleDouble.One); + AssertBits(DoubleDouble.Pow(DoubleDouble.NegativeOne, fractional), DoubleDouble.NaN); + AssertBits(DoubleDouble.Pow(new DoubleDouble(-2.0), odd), 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)); - 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)); - 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)); - 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)); - 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)); - 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)); } } 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] @@ -84,9 +84,9 @@ public class PreciseMathRealPowSpecialTests foreach (double sign in new[] { -1.0, 1.0 }) { 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)); - 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)); } } @@ -102,11 +102,11 @@ public class PreciseMathRealPowSpecialTests { 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 // separately checks these algorithms against independent references. - AssertBits(DDMath.Pow(value, (double)exponent), expected); - AssertBits(DDMath.Pow(value, new DoubleDouble(exponent)), expected); + AssertBits(DoubleDouble.Pow(value, (double)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 // constant-reference generator. Compare components, not high+low in double. - foreach (DoubleDouble actual in new[] { DDMath.Pow(new DoubleDouble(2.0), 0.5), - DDMath.Pow(new DoubleDouble(2.0), new DoubleDouble(0.5)) }) + foreach (DoubleDouble actual in new[] { DoubleDouble.Pow(new DoubleDouble(2.0), 0.5), + DoubleDouble.Pow(new DoubleDouble(2.0), new DoubleDouble(0.5)) }) { actual.High.ShouldBe(1.4142135623730951); 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.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), - DDMath.Pow(new DoubleDouble(4.0), new DoubleDouble(exponent)) }) + foreach (DoubleDouble actual in new[] { DoubleDouble.Pow(new DoubleDouble(4.0), exponent), + DoubleDouble.Pow(new DoubleDouble(4.0), new DoubleDouble(exponent)) }) { BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(expectedBits); BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L); @@ -149,8 +149,8 @@ public class PreciseMathRealPowSpecialTests { bool grows = (value > 1.0) == (exponent > 0.0); DoubleDouble expected = new(grows ? double.PositiveInfinity : 0.0); - AssertBits(DDMath.Pow(new DoubleDouble(value), exponent), expected); - AssertBits(DDMath.Pow(new DoubleDouble(value), new DoubleDouble(exponent)), expected); + AssertBits(DoubleDouble.Pow(new DoubleDouble(value), exponent), expected); + AssertBits(DoubleDouble.Pow(new DoubleDouble(value), new DoubleDouble(exponent)), expected); } } } diff --git a/1-tests/Just.PreciseMath.Tests/PreciseMathRealPowTests.cs b/1-tests/Just.PreciseMath.Tests/PreciseMathRealPowTests.cs index 607c5a4..7224dcd 100644 --- a/1-tests/Just.PreciseMath.Tests/PreciseMathRealPowTests.cs +++ b/1-tests/Just.PreciseMath.Tests/PreciseMathRealPowTests.cs @@ -21,7 +21,7 @@ public class PreciseMathRealPowTests (Units(input.High) + Units(input.Low)).ShouldBe(Units(high) + Units(low)); BigInteger exponentUnits = Units(exponentHigh) + Units(scalar ? 0.0 : exponentLow); (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})"; AssertReferenceBound(actual, reference, exponentUnits, context); } @@ -36,9 +36,9 @@ public class PreciseMathRealPowTests // These references are exact rational identities. The general ln/exp // implementation is approximate, so do not require universal bit equality. 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"); - 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"); } @@ -55,17 +55,17 @@ public class PreciseMathRealPowTests { foreach (int exponent in new[] { -63, -7, -1, 0, 1, 7, 63 }) { - DoubleDouble expected = DDMath.Pow(input, exponent); - AssertBits(DDMath.Pow(input, (double)exponent), expected); - AssertBits(DDMath.Pow(input, new DoubleDouble(exponent)), expected); + DoubleDouble expected = DoubleDouble.Pow(input, exponent); + AssertBits(DoubleDouble.Pow(input, (double)exponent), expected); + AssertBits(DoubleDouble.Pow(input, new DoubleDouble(exponent)), expected); } } foreach (int exponent in new[] { int.MinValue, int.MaxValue }) { DoubleDouble input = DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -54)); - DoubleDouble expected = DDMath.Pow(input, exponent); - AssertBits(DDMath.Pow(input, (double)exponent), expected); - AssertBits(DDMath.Pow(input, new DoubleDouble(exponent)), expected); + DoubleDouble expected = DoubleDouble.Pow(input, exponent); + AssertBits(DoubleDouble.Pow(input, (double)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 }) { DoubleDouble input = DoubleDouble.FromComponents(high, low); - AssertBits(DDMath.Pow(input, 1.0), input); - AssertBits(DDMath.Pow(input, DoubleDouble.One), input); + AssertBits(DoubleDouble.Pow(input, 1.0), input); + AssertBits(DoubleDouble.Pow(input, DoubleDouble.One), input); } } } diff --git a/1-tests/Just.PreciseMath.Tests/PreciseMathReciprocalTests.cs b/1-tests/Just.PreciseMath.Tests/PreciseMathReciprocalTests.cs index ded089b..8c4c0cf 100644 --- a/1-tests/Just.PreciseMath.Tests/PreciseMathReciprocalTests.cs +++ b/1-tests/Just.PreciseMath.Tests/PreciseMathReciprocalTests.cs @@ -12,7 +12,7 @@ public class PreciseMathReciprocalTests foreach (double low in new[] { 0.0, Math.ScaleB(1.0, -54), -Math.ScaleB(1.0, -54) }) { DoubleDouble input = DoubleDouble.FromComponents(high, low); - DoubleDouble actual = DDMath.Reciprocal(input); + DoubleDouble actual = DoubleDouble.Reciprocal(input); actual.Low.ShouldNotBe(0.0); AssertReciprocalBound(input, actual); } @@ -25,7 +25,7 @@ public class PreciseMathReciprocalTests double.NaN, BitConverter.Int64BitsToDouble(0x7ff0000000000001L)]; 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. double expected = 1.0 / value; BitConverter.DoubleToInt64Bits(actual.High).ShouldBe( @@ -44,7 +44,7 @@ public class PreciseMathReciprocalTests { 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)); BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L); DoubleDouble.IsCanonical(actual).ShouldBeTrue(); @@ -104,7 +104,7 @@ public class PreciseMathReciprocalTests { double low = Math.ScaleB(lowSign, exponent); DoubleDouble input = DoubleDouble.FromComponents(sign, low); - DoubleDouble actual = DDMath.Reciprocal(input); + DoubleDouble actual = DoubleDouble.Reciprocal(input); actual.High.ShouldBe(sign); actual.Low.ShouldBe(-low); AssertReciprocalBound(input, actual); @@ -116,7 +116,7 @@ public class PreciseMathReciprocalTests private static void AssertMatchesDivision(DoubleDouble 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})"; BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(BitConverter.DoubleToInt64Bits(expected.High), context); BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(BitConverter.DoubleToInt64Bits(expected.Low), context); diff --git a/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs b/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs index 6282fa3..372ca74 100644 --- a/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs +++ b/1-tests/Just.PreciseMath.Tests/PreciseMathSqrtTests.cs @@ -10,7 +10,7 @@ public class PreciseMathSqrtTests [InlineData(2.25, 1.5)] 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); BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L); } @@ -21,7 +21,7 @@ public class PreciseMathSqrtTests // sqrt(2^(2k)) = 2^k exactly, including the minimum subnormal input. 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)); BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(0L); } @@ -36,7 +36,7 @@ public class PreciseMathSqrtTests public void IrrationalRootsRetainMoreThanBinary64Precision(double value) { DoubleDouble input = new(value); - DoubleDouble actual = DDMath.Sqrt(input); + DoubleDouble actual = DoubleDouble.Sqrt(input); actual.Low.ShouldNotBe(0.0); 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) }) { DoubleDouble input = DoubleDouble.FromComponents(1.0, low); - DoubleDouble actual = DDMath.Sqrt(input); + DoubleDouble actual = DoubleDouble.Sqrt(input); actual.High.ShouldBe(1.0); (Math.Sign(actual.Low) == Math.Sign(low)).ShouldBeTrue(); AssertSqrtBound(input, actual); @@ -61,7 +61,7 @@ public class PreciseMathSqrtTests double.NaN, -1.0, -double.Epsilon, double.MinValue]; foreach (double value in values) { - DoubleDouble actual = DDMath.Sqrt(new DoubleDouble(value)); + DoubleDouble actual = DoubleDouble.Sqrt(new DoubleDouble(value)); double expected = Math.Sqrt(value); // Construction canonicalizes all NaNs to double.NaN, including // 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.IsNaN(DDMath.Sqrt(negative)).ShouldBeTrue(); + DoubleDouble.IsNaN(DoubleDouble.Sqrt(negative)).ShouldBeTrue(); } [Fact] @@ -92,7 +92,7 @@ public class PreciseMathSqrtTests if (boundary > 0.0) { 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 }) @@ -100,7 +100,7 @@ public class PreciseMathSqrtTests DoubleDouble input = DoubleDouble.FromComponents(high, residual); 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))]; 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)); 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 }) { 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] @@ -177,7 +177,7 @@ public class PreciseMathSqrtTests DoubleDouble input = DoubleDouble.FromComponents(high, low); 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 }) { 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.Low.ShouldBe(Math.ScaleB(sign, exponent - 1)); DoubleDouble.IsCanonical(actual).ShouldBeTrue(); @@ -358,9 +358,6 @@ public class PreciseMathSqrtTests DoubleDouble actual = DoubleDouble.Sqrt(input); actual.High.ShouldBe(high); 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); } diff --git a/1-tests/Just.PreciseMath.Tests/PreciseMathTests.cs b/1-tests/Just.PreciseMath.Tests/PreciseMathTests.cs index ba030f8..94f5792 100644 --- a/1-tests/Just.PreciseMath.Tests/PreciseMathTests.cs +++ b/1-tests/Just.PreciseMath.Tests/PreciseMathTests.cs @@ -10,7 +10,7 @@ public class PreciseMathTests foreach (double scalar in scalars) { // 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); DoubleDouble.IsCanonical(actual).ShouldBeTrue(); } @@ -19,8 +19,8 @@ public class PreciseMathTests { DoubleDouble positive = DoubleDouble.FromComponents(1.0, low); DoubleDouble negative = DoubleDouble.FromComponents(-1.0, -low); - AssertBits(DDMath.Abs(positive), 1.0, low); - AssertBits(DDMath.Abs(negative), 1.0, low); + AssertBits(DoubleDouble.Abs(positive), 1.0, low); + AssertBits(DoubleDouble.Abs(negative), 1.0, low); } } diff --git a/README.md b/README.md index 356eed5..eab9880 100644 --- a/README.md +++ b/README.md @@ -22,8 +22,11 @@ already-normalized results. Mathematical constants use precomputed high/low pair checked against independently computed high-precision values; accessing them does not perform double-double arithmetic or allocate on the heap. -- Arithmetic: unary `+`/`-`, binary `+`, `-`, `*`, `/`, and both operand orders with - a `double`. Addition retains residuals under cancellation; multiplication uses +- Arithmetic: static `Add`, `Subtract`, `Multiply`, and `Divide` methods accept + `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 use specialized scalar paths rather than promoting the scalar to `DoubleDouble`. 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 approximate DD arithmetic, so conversions are not guaranteed exact round trips. 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 @@ -102,16 +105,40 @@ The additional `InvSqrt` helper and three-argument `Hypot` overload live alongsi them. Exponential functions live in `DoubleDouble.ExponentialFunctions.cs`: `Exp`, `Exp2`, `Exp10`, `ExpM1`, `Exp2M1`, and `Exp10M1` complete `IExponentialFunctions`, including cancellation-safe overrides of -the interface's default minus-one methods. `DDMath` exposes thin forwarding -wrappers for roots and exponentials; logarithm, power, and reciprocal kernels -remain in `DDMath`. +the interface's default minus-one methods. Natural logarithm and power kernels +live in `DoubleDouble.Log.cs` and `DoubleDouble.Pow.cs`; reciprocal lives in +`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 positive zero and either infinity to positive infinity, and returns canonical NaN. - It shares the existing `DoubleDouble.Abs` implementation. -- `Square(DoubleDouble)`: forwards to the specialized `DoubleDouble.Square` kernel. +- `Square(DoubleDouble)`: uses a specialized double-double squaring kernel. It combines the equal cross terms with FMA, normalizes, then incorporates the low-component square without losing selected exactly representable tails. 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 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 reciprocal = DDMath.Reciprocal(new DoubleDouble(3.0)); -DoubleDouble magnitude = DDMath.Abs(-root); -DoubleDouble smallPower = DDMath.Pow(new DoubleDouble(2.0), -1024); +DoubleDouble reciprocal = DoubleDouble.Reciprocal(new DoubleDouble(3.0)); +DoubleDouble magnitude = DoubleDouble.Abs(-root); +DoubleDouble smallPower = DoubleDouble.Pow(new DoubleDouble(2.0), -1024); DoubleDouble exponential = DoubleDouble.Exp(new DoubleDouble(1.0)); DoubleDouble binaryPower = DoubleDouble.Exp2(new DoubleDouble(-1000.0)); DoubleDouble decimalPower = DoubleDouble.Exp10(new DoubleDouble(0.5)); DoubleDouble tinyChange = DoubleDouble.ExpM1(new DoubleDouble(1e-30)); -DoubleDouble logarithm = DDMath.Log(new DoubleDouble(10.0)); -DoubleDouble fractionalPower = DDMath.Pow(new DoubleDouble(2.0), 0.5); +DoubleDouble logarithm = DoubleDouble.Log(new DoubleDouble(10.0)); +DoubleDouble fractionalPower = DoubleDouble.Pow(new DoubleDouble(2.0), 0.5); 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 @@ -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. Hypotenuse tests use exact squared inequalities with `2^-100` relative error plus 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 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 @@ -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 are checked at 450/650 digits, including adjacent lows at true range boundaries, series transitions, and tiny inputs. Additional component checks pin sparse -base-two corrections; constrained generic calls and facade calls match direct -calls bit-for-bit. Integer-power tests use exact rational -references and high-precision fixtures for large exponents; their tested absolute +base-two corrections; constrained generic calls match direct calls bit-for-bit. +Integer-power tests use exact rational references and high-precision fixtures for +large exponents; their tested absolute error bound is `|exact result| * (|exponent| + 1) * 2^-100 + 2^-1074`, not uniform relative accuracy independent of the exponent. These functions are approximate; 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 -Natural `DDMath.Log`, the complete exponential family on `DoubleDouble`, and all -three `DDMath.Pow` overloads are implemented. +Natural `DoubleDouble.Log`, the complete exponential family on `DoubleDouble`, and all +three `DoubleDouble.Pow` overloads are implemented. Logarithms in other bases, generic-math interfaces beyond `ISignedNumber`, `IFloatingPointConstants`, `IRootFunctions`, and `IExponentialFunctions`, additional text formats/shortest- round-trip formatting, and non-arithmetic performance benchmarks remain deferred.