diff --git a/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs b/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs index 90d091e..ca86a5f 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs @@ -1,17 +1,9 @@ namespace Just.PreciseMath; public readonly partial struct DoubleDouble : - IAdditiveIdentity, - IMultiplicativeIdentity, - IUnaryPlusOperators, - IUnaryNegationOperators, - IAdditionOperators, IAdditionOperators, - ISubtractionOperators, ISubtractionOperators, - IMultiplyOperators, IMultiplyOperators, - IDivisionOperators, IDivisionOperators { /// Returns the operand unchanged. diff --git a/0-source/Just.PreciseMath/DoubleDouble.Comparison.cs b/0-source/Just.PreciseMath/DoubleDouble.Comparison.cs index c4f439c..2d0a8d9 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.Comparison.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.Comparison.cs @@ -1,13 +1,25 @@ namespace Just.PreciseMath; public readonly partial struct DoubleDouble : IComparable, IComparable, + IEquatable, + IEqualityOperators, IComparisonOperators { - /// Returns the normalized components. - public void Decompose(out double high, out double low) + /// + /// Tests numerical equality; NaN operands are never equal. + /// + [Pure, MethodImpl(MethodImplOptions.AggressiveInlining)] + public static bool operator ==(DoubleDouble left, DoubleDouble right) { - high = _high; - low = _low; + return left._high == right._high && left._low == right._low; + } + /// + /// Tests numerical inequality; NaN operands are always unequal. + /// + [Pure, MethodImpl(MethodImplOptions.AggressiveInlining)] + public static bool operator !=(DoubleDouble left, DoubleDouble right) + { + return !(left == right); } /// Orders NaN before other values, and compares finite values using both components. diff --git a/0-source/Just.PreciseMath/DoubleDouble.Constants.cs b/0-source/Just.PreciseMath/DoubleDouble.Constants.cs new file mode 100644 index 0000000..a39a63b --- /dev/null +++ b/0-source/Just.PreciseMath/DoubleDouble.Constants.cs @@ -0,0 +1,144 @@ +namespace Just.PreciseMath; + +public readonly partial struct DoubleDouble : + IFloatingPointConstants +{ + /// + /// Represents an additive identity value. + /// + public static DoubleDouble AdditiveIdentity => Zero; + /// + /// Represents a multiplicative identity value. + /// + public static DoubleDouble MultiplicativeIdentity => One; + + /// + /// Represents a value that is not a number (NaN). + /// + public static DoubleDouble NaN => new(double.NaN); + + /// Gets positive infinity with a positive-zero low component. + public static DoubleDouble PositiveInfinity => new(double.PositiveInfinity); + + /// Gets negative infinity with a positive-zero low component. + public static DoubleDouble NegativeInfinity => new(double.NegativeInfinity); + + /// + /// Represents a unit value. + /// + public static DoubleDouble One => new(1.0, 0); + /// + /// Represents a negative unit value. + /// + public static DoubleDouble NegativeOne => new(-1.0, 0); + /// + /// Represents a zero value. + /// + public static DoubleDouble Zero => new(); + + /// Gets negative zero with a positive-zero low component. + /// Compares equal to while preserving the high component's sign bit. + public static DoubleDouble NegativeZero => new(-0.0); + + /// Gets the smallest positive representable value, 2^-1074. + /// + /// The high component is the minimum positive binary64 subnormal and the low + /// component is positive zero. This is not machine epsilon or a relative-error tolerance. + /// + public static DoubleDouble Epsilon => new(double.Epsilon); + + // Mathematical constants store normalized, precomputed binary64 pairs: the + // nearest high and then the nearest residual of the high-precision value. + // No DD arithmetic, parsing, normalization or heap allocation occurs on access. + // Reproduce with the tests' ReferenceData/generate_constants.py. + /// + /// Represents the ratio of the circumference of a circle to its diameter, specified by the constant, π. + /// + public static DoubleDouble Pi => new(3.141592653589793, 1.2246467991473532e-16); + /// + /// Represents the natural logarithmic base, specified by the constant, e. + /// + public static DoubleDouble E => new(2.718281828459045, 1.4456468917292502e-16); + /// + /// Represents the natural logarithm of value 2. + /// + public static DoubleDouble Ln2 => new(0.6931471805599453, 2.3190468138462996e-17); + + /// Gets τ = 2π, the angle of one full turn in radians. + public static DoubleDouble Tau => new(6.283185307179586, 2.4492935982947064e-16); + + /// Gets π/2, the angle of 90 degrees in radians. + public static DoubleDouble PiOver2 => new(1.5707963267948966, 6.123233995736766e-17); + + /// Gets π/3, the angle of 60 degrees in radians. + public static DoubleDouble PiOver3 => new(1.0471975511965979, -1.072081766451091e-16); + + /// Gets π/4, the angle of 45 degrees in radians. + public static DoubleDouble PiOver4 => new(0.7853981633974483, 3.061616997868383e-17); + + /// Gets π/6, the angle of 30 degrees in radians. + public static DoubleDouble PiOver6 => new(0.5235987755982989, -5.360408832255455e-17); + + /// Gets 1/π. + public static DoubleDouble InvPi => new(0.3183098861837907, -1.9678676675182486e-17); + + /// Gets 1/(2π), the factor for converting radians to turns. + public static DoubleDouble InvTau => new(0.15915494309189535, -9.839338337591243e-18); + + /// Gets π/180. Multiply an angle in degrees by this value to obtain radians. + public static DoubleDouble DegToRad => new(0.017453292519943295, 2.9486522708701687e-19); + + /// Gets 180/π. Multiply an angle in radians by this value to obtain degrees. + public static DoubleDouble RadToDeg => new(57.29577951308232, -1.9878495670576283e-15); + + /// Gets 1/e = exp(-1). + public static DoubleDouble InvE => new(0.36787944117144233, -1.2428753672788363e-17); + + /// Gets ln(10), the factor for converting base-10 logarithms to natural logarithms. + public static DoubleDouble Ln10 => new(2.302585092994046, -2.1707562233822494e-16); + + /// Gets log₂(e) = 1/ln(2), the factor for converting natural logarithms to base 2. + public static DoubleDouble Log2E => new(1.4426950408889634, 2.0355273740931033e-17); + + /// Gets log₁₀(e) = 1/ln(10), the factor for converting natural logarithms to base 10. + public static DoubleDouble Log10E => new(0.4342944819032518, 1.098319650216765e-17); + + /// Gets log₂(10), the factor for converting base-10 logarithms to base 2. + public static DoubleDouble Log2Of10 => new(3.321928094887362, 1.661617516973592e-16); + + /// Gets log₁₀(2), the factor for converting base-2 logarithms to base 10. + public static DoubleDouble Log10Of2 => new(0.3010299956639812, -2.8037281277851704e-18); + + /// Gets √2. + public static DoubleDouble Sqrt2 => new(1.4142135623730951, -9.667293313452913e-17); + + /// Gets √3. + public static DoubleDouble Sqrt3 => new(1.7320508075688772, 1.0035084221806903e-16); + + /// Gets √5. + public static DoubleDouble Sqrt5 => new(2.23606797749979, -1.0864230407365012e-16); + + /// Gets 1/√2, also the sine and cosine of π/4. + public static DoubleDouble InvSqrt2 => new(0.7071067811865476, -4.833646656726457e-17); + + /// Gets 1/√3, also the tangent of π/6. + public static DoubleDouble InvSqrt3 => new(0.5773502691896257, 3.3450280739356345e-17); + + /// Gets √π, the Gaussian integral over the real line for exp(-x²). + public static DoubleDouble SqrtPi => new(1.772453850905516, -7.666586499825799e-17); + + /// Gets 1/√π, a Gaussian normalization factor. + public static DoubleDouble InvSqrtPi => new(0.5641895835477563, 7.66772980658294e-18); + + /// Gets 2/√π, the normalization factor in the error-function integral. + public static DoubleDouble TwoInvSqrtPi => new(1.1283791670955126, 1.533545961316588e-17); + + /// Gets √(2π), used in Gaussian integrals and Stirling's approximation. + public static DoubleDouble SqrtTau => new(2.5066282746310007, -1.8328579980459167e-16); + + /// Gets 1/√(2π), the standard normal probability density's normalization factor. + public static DoubleDouble InvSqrtTau => new(0.3989422804014327, -2.49232720227773e-17); + + /// Gets the golden ratio φ = (1 + √5)/2. + public static DoubleDouble GoldenRatio => new(1.618033988749895, -5.432115203682506e-17); +} diff --git a/0-source/Just.PreciseMath/DoubleDouble.cs b/0-source/Just.PreciseMath/DoubleDouble.cs index 4a955a2..39ac13e 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.cs @@ -12,13 +12,12 @@ namespace Just.PreciseMath; /// rounded. Equals treats NaNs as equal for collections, while operators do not. /// public readonly partial struct DoubleDouble : - IEquatable, - IEqualityOperators, ISignedNumber { internal readonly double _high; internal readonly double _low; + /// /// Stores trusted components without normalization or validation. /// @@ -72,129 +71,6 @@ public readonly partial struct DoubleDouble : { } - #region Static constants - /// - /// Represents an additive identity value. - /// - public static DoubleDouble AdditiveIdentity => Zero; - /// - /// Represents a multiplicative identity value. - /// - public static DoubleDouble MultiplicativeIdentity => One; - - /// - /// Represents a value that is not a number (NaN). - /// - public static DoubleDouble NaN => new(double.NaN); - /// - /// Represents a unit value. - /// - public static DoubleDouble One => new(1.0, 0); - /// - /// Represents a negative unit value. - /// - public static DoubleDouble NegativeOne => new(-1.0, 0); - /// - /// Represents a zero value. - /// - public static DoubleDouble Zero => new(); - - // Mathematical constants store normalized, precomputed binary64 pairs: the - // nearest high and then the nearest residual of the high-precision value. - // No DD arithmetic, parsing, normalization or heap allocation occurs on access. - // Reproduce with the tests' ReferenceData/generate_constants.py. - /// - /// Represents the ratio of the circumference of a circle to its diameter, specified by the constant, π. - /// - public static DoubleDouble Pi => new(3.141592653589793, 1.2246467991473532e-16); - /// - /// Represents the natural logarithmic base, specified by the constant, e. - /// - public static DoubleDouble E => new(2.718281828459045, 1.4456468917292502e-16); - /// - /// Represents the natural logarithm of value 2. - /// - public static DoubleDouble Ln2 => new(0.6931471805599453, 2.3190468138462996e-17); - - /// Gets τ = 2π, the angle of one full turn in radians. - public static DoubleDouble Tau => new(6.283185307179586, 2.4492935982947064e-16); - - /// Gets π/2, the angle of 90 degrees in radians. - public static DoubleDouble PiOver2 => new(1.5707963267948966, 6.123233995736766e-17); - - /// Gets π/3, the angle of 60 degrees in radians. - public static DoubleDouble PiOver3 => new(1.0471975511965979, -1.072081766451091e-16); - - /// Gets π/4, the angle of 45 degrees in radians. - public static DoubleDouble PiOver4 => new(0.7853981633974483, 3.061616997868383e-17); - - /// Gets π/6, the angle of 30 degrees in radians. - public static DoubleDouble PiOver6 => new(0.5235987755982989, -5.360408832255455e-17); - - /// Gets 1/π. - public static DoubleDouble InvPi => new(0.3183098861837907, -1.9678676675182486e-17); - - /// Gets 1/(2π), the factor for converting radians to turns. - public static DoubleDouble InvTau => new(0.15915494309189535, -9.839338337591243e-18); - - /// Gets π/180. Multiply an angle in degrees by this value to obtain radians. - public static DoubleDouble DegToRad => new(0.017453292519943295, 2.9486522708701687e-19); - - /// Gets 180/π. Multiply an angle in radians by this value to obtain degrees. - public static DoubleDouble RadToDeg => new(57.29577951308232, -1.9878495670576283e-15); - - /// Gets 1/e = exp(-1). - public static DoubleDouble InvE => new(0.36787944117144233, -1.2428753672788363e-17); - - /// Gets ln(10), the factor for converting base-10 logarithms to natural logarithms. - public static DoubleDouble Ln10 => new(2.302585092994046, -2.1707562233822494e-16); - - /// Gets log₂(e) = 1/ln(2), the factor for converting natural logarithms to base 2. - public static DoubleDouble Log2E => new(1.4426950408889634, 2.0355273740931033e-17); - - /// Gets log₁₀(e) = 1/ln(10), the factor for converting natural logarithms to base 10. - public static DoubleDouble Log10E => new(0.4342944819032518, 1.098319650216765e-17); - - /// Gets log₂(10), the factor for converting base-10 logarithms to base 2. - public static DoubleDouble Log2Of10 => new(3.321928094887362, 1.661617516973592e-16); - - /// Gets log₁₀(2), the factor for converting base-2 logarithms to base 10. - public static DoubleDouble Log10Of2 => new(0.3010299956639812, -2.8037281277851704e-18); - - /// Gets √2. - public static DoubleDouble Sqrt2 => new(1.4142135623730951, -9.667293313452913e-17); - - /// Gets √3. - public static DoubleDouble Sqrt3 => new(1.7320508075688772, 1.0035084221806903e-16); - - /// Gets √5. - public static DoubleDouble Sqrt5 => new(2.23606797749979, -1.0864230407365012e-16); - - /// Gets 1/√2, also the sine and cosine of π/4. - public static DoubleDouble InvSqrt2 => new(0.7071067811865476, -4.833646656726457e-17); - - /// Gets 1/√3, also the tangent of π/6. - public static DoubleDouble InvSqrt3 => new(0.5773502691896257, 3.3450280739356345e-17); - - /// Gets √π, the Gaussian integral over the real line for exp(-x²). - public static DoubleDouble SqrtPi => new(1.772453850905516, -7.666586499825799e-17); - - /// Gets 1/√π, a Gaussian normalization factor. - public static DoubleDouble InvSqrtPi => new(0.5641895835477563, 7.66772980658294e-18); - - /// Gets 2/√π, the normalization factor in the error-function integral. - public static DoubleDouble TwoInvSqrtPi => new(1.1283791670955126, 1.533545961316588e-17); - - /// Gets √(2π), used in Gaussian integrals and Stirling's approximation. - public static DoubleDouble SqrtTau => new(2.5066282746310007, -1.8328579980459167e-16); - - /// Gets 1/√(2π), the standard normal probability density's normalization factor. - public static DoubleDouble InvSqrtTau => new(0.3989422804014327, -2.49232720227773e-17); - - /// Gets the golden ratio φ = (1 + √5)/2. - public static DoubleDouble GoldenRatio => new(1.618033988749895, -5.432115203682506e-17); - #endregion - /// /// High part of DoubleDouble /// @@ -224,20 +100,10 @@ public readonly partial struct DoubleDouble : return HashCode.Combine(_high, _low); } - /// - /// Tests numerical equality; NaN operands are never equal. - /// - [Pure, MethodImpl(MethodImplOptions.AggressiveInlining)] - public static bool operator ==(DoubleDouble left, DoubleDouble right) + /// Returns the normalized components. + public void Decompose(out double high, out double low) { - return left._high == right._high && left._low == right._low; - } - /// - /// Tests numerical inequality; NaN operands are always unequal. - /// - [Pure, MethodImpl(MethodImplOptions.AggressiveInlining)] - public static bool operator !=(DoubleDouble left, DoubleDouble right) - { - return !(left == right); + high = _high; + low = _low; } } diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleTests.cs index 9691a30..790f0e7 100644 --- a/1-tests/Just.PreciseMath.Tests/DoubleDoubleTests.cs +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleTests.cs @@ -1,3 +1,4 @@ +using System.Numerics; using Shouldly; using Xunit; @@ -5,6 +6,71 @@ namespace Just.PreciseMath.Tests; public class DoubleDoubleTests { + [Fact] + public void EpsilonIsTheMinimumPositiveBinary64Subnormal() + { + // Each finite component is an integer multiple of 2^-1074, so their + // exact sum cannot have a smaller positive quantum. This is not a + // relative-error tolerance or a fixed significand spacing near one. + DoubleDouble value = DoubleDouble.Epsilon; + BitConverter.DoubleToInt64Bits(value.High).ShouldBe(1L); + BitConverter.DoubleToInt64Bits(value.Low).ShouldBe(0L); + DoubleDouble.IsCanonical(value).ShouldBeTrue(); + DoubleDouble.IsFinite(value).ShouldBeTrue(); + DoubleDouble.IsSubnormal(value).ShouldBeTrue(); + DoubleDouble.IsNormal(value).ShouldBeFalse(); + DoubleDouble.IsPositive(value).ShouldBeTrue(); + (value > DoubleDouble.Zero).ShouldBeTrue(); + } + + [Theory] + [InlineData("NegativeZero", -0.0)] + [InlineData("PositiveInfinity", double.PositiveInfinity)] + [InlineData("NegativeInfinity", double.NegativeInfinity)] + [InlineData("NaN", double.NaN)] + public void SpecialConstantsHaveCanonicalComponentBits(string name, double expectedHigh) + { + DoubleDouble value = name switch + { + "NegativeZero" => DoubleDouble.NegativeZero, + "PositiveInfinity" => DoubleDouble.PositiveInfinity, + "NegativeInfinity" => DoubleDouble.NegativeInfinity, + "NaN" => DoubleDouble.NaN, + _ => throw new ArgumentOutOfRangeException(nameof(name)), + }; + BitConverter.DoubleToInt64Bits(value.High).ShouldBe(BitConverter.DoubleToInt64Bits(expectedHigh)); + BitConverter.DoubleToInt64Bits(value.Low).ShouldBe(0L); + DoubleDouble.IsCanonical(value).ShouldBeTrue(); + DoubleDouble.IsNaN(value).ShouldBe(double.IsNaN(expectedHigh)); + DoubleDouble.IsFinite(value).ShouldBe(double.IsFinite(expectedHigh)); + DoubleDouble.IsPositiveInfinity(value).ShouldBe(double.IsPositiveInfinity(expectedHigh)); + DoubleDouble.IsNegativeInfinity(value).ShouldBe(double.IsNegativeInfinity(expectedHigh)); + DoubleDouble.IsNegative(value).ShouldBe(double.IsNegative(expectedHigh)); + DoubleDouble.IsZero(value).ShouldBe(expectedHigh == 0.0); + } + + [Fact] + public void NegativeZeroConstantPreservesZeroEqualityAndHashing() + { + DoubleDouble value = DoubleDouble.NegativeZero; + (value == DoubleDouble.Zero).ShouldBeTrue(); + value.Equals(DoubleDouble.Zero).ShouldBeTrue(); + value.CompareTo(DoubleDouble.Zero).ShouldBe(0); + value.GetHashCode().ShouldBe(DoubleDouble.Zero.GetHashCode()); + BitConverter.DoubleToInt64Bits((-value).High).ShouldBe(0L); + BitConverter.DoubleToInt64Bits((-value).Low).ShouldBe(0L); + BitConverter.DoubleToInt64Bits(DoubleDouble.Zero.High).ShouldBe(0L); + } + + [Fact] + public void FloatingPointConstantsSupportGenericDispatch() + { + (DoubleDouble e, DoubleDouble pi, DoubleDouble tau) = GetFloatingPointConstants(); + e.ShouldBe(DoubleDouble.E); + pi.ShouldBe(DoubleDouble.Pi); + tau.ShouldBe(DoubleDouble.Tau); + } + [Theory] [InlineData("Pi", 3.141592653589793, 1.2246467991473532e-16)] [InlineData("E", 2.718281828459045, 1.4456468917292502e-16)] @@ -184,4 +250,9 @@ public class DoubleDoubleTests result.High.ShouldBe(high); result.Low.ShouldBe(low); } + + private static (T E, T Pi, T Tau) GetFloatingPointConstants() where T : IFloatingPointConstants + { + return (T.E, T.Pi, T.Tau); + } } diff --git a/README.md b/README.md index 4d8ba4b..25fa923 100644 --- a/README.md +++ b/README.md @@ -45,6 +45,13 @@ with exact component checks for selected representable cases. Near underflow, extended precision necessarily decreases; overflow produces infinity. Performance of the allocating exponent-boundary path is not covered by the basic benchmarks. +Named value constants include `Zero`, `NegativeZero`, `One`, `NegativeOne`, `NaN`, +`PositiveInfinity`, `NegativeInfinity`, and `Epsilon`. `Epsilon` is the smallest +positive representable value, `2^-1074` (the same as `double.Epsilon`), **not** a +relative-error tolerance or machine epsilon. These constants have a positive-zero +low component; `NegativeZero` preserves the high sign bit while comparing and +hashing equal to `Zero`. + ## Predefined mathematical constants All constants below are static `DoubleDouble` properties. Each stores the nearest @@ -52,6 +59,10 @@ binary64 high component followed by the nearest binary64 residual, rather than calculating a ratio, root, or logarithm on access. Names use PascalCase, including `Pi`, `E`, and `Ln2`. +`DoubleDouble` implements `IFloatingPointConstants` for generic +access to `E`, `Pi`, and `Tau`; this does not imply support for the full +`IFloatingPointIeee754` interface. + | Group | Properties and values | |---|---| | Circle and common angles | `Pi` (π), `Tau` (2π), `PiOver2`, `PiOver3`, `PiOver4`, `PiOver6` | @@ -280,9 +291,9 @@ bool success = DoubleDouble.TryParse("1.25e-2".AsSpan(), CultureInfo.InvariantCu ## Deferred scope Natural `DDMath.Log`, `Exp`, and all three `Pow` overloads are implemented. -Logarithms in other bases, generic-math interfaces beyond `ISignedNumber`, -additional text formats/general round-trip formatting, and non-arithmetic performance -benchmarks remain deferred. +Logarithms in other bases, generic-math interfaces beyond `ISignedNumber` and +`IFloatingPointConstants`, additional text formats/general round-trip formatting, +and non-arithmetic performance benchmarks remain deferred. Replacing allocating arithmetic boundary fallbacks is also deferred; the current `BigInteger` paths remain in place. That optimization does not require removing `BigInteger` from conversions, parsing, formatting, or independent test oracles.