From 4a7ce1b8ed00f62e1616a802a95c16538eb09778 Mon Sep 17 00:00:00 2001 From: just Date: Thu, 24 Sep 2026 02:25:35 +0400 Subject: [PATCH] modulus implementation --- .../DoubleDouble.Arithmetic.cs | 163 ++++++++- .../DoubleDouble.NumberBase.cs | 6 + .../DoubleDouble.Operators.cs | 9 + 0-source/Just.PreciseMath/DoubleDouble.cs | 2 +- .../DoubleDoubleModulusTests.cs | 337 ++++++++++++++++++ .../DoubleDoubleNumberTests.cs | 99 +++++ README.md | 38 +- 7 files changed, 646 insertions(+), 8 deletions(-) create mode 100644 1-tests/Just.PreciseMath.Tests/DoubleDoubleModulusTests.cs create mode 100644 1-tests/Just.PreciseMath.Tests/DoubleDoubleNumberTests.cs diff --git a/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs b/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs index 56f2da3..e6b5213 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs @@ -4,7 +4,8 @@ public readonly partial struct DoubleDouble : IAdditionOperators, ISubtractionOperators, IMultiplyOperators, - IDivisionOperators + IDivisionOperators, + IModulusOperators { /// Returns the square using a specialized double-double product. /// @@ -152,6 +153,29 @@ public readonly partial struct DoubleDouble : return PreciseMathHelper.AddScalarFinite(partial._high, partial._low, finalCorrection); } + /// Returns the remainder after division with the quotient truncated toward zero. + /// + /// Uses both components to select the exact integer quotient, even when that + /// quotient cannot be represented by a floating-point value. The remainder is + /// rounded to a normalized pair, not computed from an approximate division. + /// Its sign, including an exact zero, follows the dividend; this is C# remainder, + /// not Euclidean modulo or IEEE nearest-integer remainder. A zero divisor, + /// infinite dividend, or NaN operand returns canonical NaN. A finite dividend + /// modulo infinity is unchanged. Ordinary finite cases use error-free floating + /// expansions on the stack. Cases outside that kernel's safe quotient/range + /// bounds fall back to bounded BigInteger intermediates. + /// + [Pure] + public static DoubleDouble Modulus(DoubleDouble left, DoubleDouble right) => ModulusCore(left._high, left._low, right._high, right._low); + + /// + [Pure] + public static DoubleDouble Modulus(DoubleDouble left, double right) => ModulusCore(left._high, left._low, right, 0.0); + + /// + [Pure] + public static DoubleDouble Modulus(double left, DoubleDouble right) => ModulusCore(left, 0.0, right._high, right._low); + /// 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 Add(DoubleDouble left, double right) => PreciseMathHelper.AddScalar(left._high, left._low, right); @@ -240,4 +264,141 @@ public readonly partial struct DoubleDouble : return PreciseMathHelper.DivideScalarFinite(left, right); } + + // Each pair is normalized/canonical, except that a scalar high may be any NaN. + // Sharing components avoids promoting mixed operands to the DD/DD overload. + private static DoubleDouble ModulusCore(double leftHigh, double leftLow, double rightHigh, double rightLow) + { + if (!double.IsFinite(leftHigh) || double.IsNaN(rightHigh) || rightHigh == 0.0) + { + return NaN; + } + if (leftHigh == 0.0 || double.IsInfinity(rightHigh)) + { + return new DoubleDouble(leftHigh, leftLow); + } + if (leftLow == 0.0 && rightLow == 0.0) + { + return new DoubleDouble(leftHigh % rightHigh); + } + + bool negative = double.IsNegative(leftHigh); + double aHigh = Math.Abs(leftHigh); + double aLow = negative ? -leftLow : leftLow; + double bHigh = Math.Abs(rightHigh); + double bLow = double.IsNegative(rightHigh) ? -rightLow : rightLow; + if (aHigh < bHigh || (aHigh == bHigh && aLow < bLow)) + { + return new DoubleDouble(leftHigh, leftLow); + } + // Equal highs imply a quotient of one (including equal operands). + // Subtract just the lows, avoiding large-high transforms at the range limit. + if (aHigh == bHigh) + { + DoubleDouble difference = FromComponents(aLow, -bLow); + return negative ? -difference : difference; + } + if (TryModulusFinite(aHigh, aLow, bHigh, bLow, out DoubleDouble remainder)) + { + return negative ? -remainder : remainder; + } + + return ModulusBoundary(leftHigh, leftLow, rightHigh, rightLow); + } + + // Positive normalized operands, a > b and distinct highs. Use the quotient + // only as an integer estimate: exact expansions certify 0 <= a - q*b < b + // BEFORE rounding the result. Never infer the quotient from a rounded DD remainder. + private static bool TryModulusFinite(double aHigh, double aLow, double bHigh, double bLow, out DoubleDouble result) + { + result = default; + double quotient = Math.Truncate(aHigh / bHigh); + // Integer products have no bits below epsilon, so their FMA residuals + // remain exact even when subnormal. Limit q so +/-1 corrections are exact. + // The high guard leaves headroom for every TwoSum intermediate, including + // comparison with b; the sum of input magnitudes stays below 2^1023. + if (quotient > 4503599627370496.0 || !PreciseMathHelper.IsAdditionWithinFastRange(aHigh)) + { + return false; + } + + Span expansion = stackalloc double[8]; + Span comparison = stackalloc double[8]; + for (int attempt = 0; attempt < 3; ++attempt) + { + (double productHigh, double errorHigh) = PreciseMathHelper.TwoMultiply(bHigh, quotient); + (double productLow, double errorLow) = PreciseMathHelper.TwoMultiply(bLow, quotient); + int length = 0; + length = GrowModulusExpansion(expansion, length, aHigh); + length = GrowModulusExpansion(expansion, length, -productHigh); + length = GrowModulusExpansion(expansion, length, aLow); + length = GrowModulusExpansion(expansion, length, -errorHigh); + length = GrowModulusExpansion(expansion, length, -productLow); + length = GrowModulusExpansion(expansion, length, -errorLow); + if (expansion[length - 1] < 0.0) + { + quotient -= 1.0; + continue; + } + + expansion[..length].CopyTo(comparison); + int comparisonLength = GrowModulusExpansion(comparison, length, -bHigh); + comparisonLength = GrowModulusExpansion(comparison, comparisonLength, -bLow); + if (comparison[comparisonLength - 1] >= 0.0) + { + quotient += 1.0; + continue; + } + + // Cancellation has already been resolved exactly. The largest term + // dominates this nonoverlapping expansion; sum the residual terms + // before normalizing to the public approximate two-component result. + double low = 0.0; + for (int i = 0; i < length - 1; ++i) + { + low += expansion[i]; + } + result = FromComponents(expansion[length - 1], low); + return true; + } + return false; + } + + // Zero-eliminating grow-expansion (TwoSum, smallest component first). + // Preserve the exact sum and increasing nonoverlap: the last term determines + // its sign. The caller reserves one extra slot for each added scalar. Writing + // in place is safe because the output cursor never overtakes the input cursor. + private static int GrowModulusExpansion(Span expansion, int length, double value) + { + int count = 0; + for (int i = 0; i < length; ++i) + { + (double sum, double error) = PreciseMathHelper.TwoAdd(value, expansion[i]); + if (error != 0.0) + { + expansion[count++] = error; + } + value = sum; + } + if (value != 0.0 || count == 0) + { + expansion[count++] = value == 0.0 ? 0.0 : value; + } + return count; + } + + // Finite normalized component pairs, with nonzero highs. All components are + // integer multiples of 2^-1074. Integer remainder therefore gives the exact + // truncating remainder without rounding/overflowing a floating-point quotient. + // Keep conversion and allocation outside the ordinary dispatch body. + [MethodImpl(MethodImplOptions.NoInlining)] + private static DoubleDouble ModulusBoundary(double leftHigh, double leftLow, double rightHigh, double rightLow) + { + BigInteger dividend = PreciseMathHelper.ArithmeticUnits(leftHigh) + PreciseMathHelper.ArithmeticUnits(leftLow); + BigInteger divisor = PreciseMathHelper.ArithmeticUnits(rightHigh) + PreciseMathHelper.ArithmeticUnits(rightLow); + BigInteger remainder = dividend % divisor; + return remainder.IsZero + ? new DoubleDouble(Math.CopySign(0.0, leftHigh)) + : PreciseMathHelper.ArithmeticFromRatio(remainder, BigInteger.One << 1074); + } } diff --git a/0-source/Just.PreciseMath/DoubleDouble.NumberBase.cs b/0-source/Just.PreciseMath/DoubleDouble.NumberBase.cs index eb3b63a..56f030e 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.NumberBase.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.NumberBase.cs @@ -42,6 +42,12 @@ public readonly partial struct DoubleDouble return IsNegative(value) ? -value : value; } + /// Returns -1 for negative values, 0 for either zero, and 1 for positive values. + /// The value is NaN. + /// Normalization makes the high component's sign authoritative. Overrides runtime-dependent INumber defaults. + [Pure] + public static int Sign(DoubleDouble value) => Math.Sign(value._high); + /// Returns the operand with greater magnitude; ties prefer positive values and NaN propagates. public static DoubleDouble MaxMagnitude(DoubleDouble x, DoubleDouble y) { diff --git a/0-source/Just.PreciseMath/DoubleDouble.Operators.cs b/0-source/Just.PreciseMath/DoubleDouble.Operators.cs index bc52eca..3b87045 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.Operators.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.Operators.cs @@ -20,6 +20,9 @@ public readonly partial struct DoubleDouble /// public static DoubleDouble operator /(DoubleDouble left, DoubleDouble right) => Divide(left, right); + /// + public static DoubleDouble operator %(DoubleDouble left, DoubleDouble right) => Modulus(left, right); + /// public static DoubleDouble operator +(DoubleDouble left, double right) => Add(left, right); @@ -44,6 +47,12 @@ public readonly partial struct DoubleDouble /// public static DoubleDouble operator /(double left, DoubleDouble right) => Divide(left, right); + /// + public static DoubleDouble operator %(DoubleDouble left, double right) => Modulus(left, right); + + /// + public static DoubleDouble operator %(double left, DoubleDouble right) => Modulus(left, right); + /// Adds one using double-double arithmetic, including its overflow and nonfinite behavior. public static DoubleDouble operator ++(DoubleDouble value) => Add(value, 1.0); diff --git a/0-source/Just.PreciseMath/DoubleDouble.cs b/0-source/Just.PreciseMath/DoubleDouble.cs index 868fbf1..b805987 100644 --- a/0-source/Just.PreciseMath/DoubleDouble.cs +++ b/0-source/Just.PreciseMath/DoubleDouble.cs @@ -12,7 +12,7 @@ namespace Just.PreciseMath; /// rounded. Equals treats NaNs as equal for collections, while operators do not. /// public readonly partial struct DoubleDouble : - ISignedNumber + ISignedNumber, INumber { internal readonly double _high; internal readonly double _low; diff --git a/1-tests/Just.PreciseMath.Tests/DoubleDoubleModulusTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleModulusTests.cs new file mode 100644 index 0000000..004e661 --- /dev/null +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleModulusTests.cs @@ -0,0 +1,337 @@ +namespace Just.PreciseMath.Tests; + +public class DoubleDoubleModulusTests +{ + [Fact] + public void ModulusRetainsTheDividendResidualThroughNamedOperatorAndGenericCalls() + { + // (2^53 + 1) % 2 = 1 exactly; using only High would return zero. + DoubleDouble left = DoubleDouble.FromComponents(9007199254740992.0, 1.0); + DoubleDouble right = new(2.0); + AssertBits(DoubleDouble.Modulus(left, right), DoubleDouble.One); + AssertBits(left % right, DoubleDouble.One); + AssertBits(GenericModulus(left, right), DoubleDouble.One); + AssertBits(DoubleDouble.Modulus(left, 2.0), DoubleDouble.One); + AssertBits(left % 2.0, DoubleDouble.One); + AssertBits(GenericScalarModulus(left, 2.0), DoubleDouble.One); + } + + [Fact] + public void ModulusRetainsTheDivisorResidualWithAScalarDividend() + { + // 2 % (1 + 2^-54) = 1 - 2^-54 exactly, with truncated quotient 1. + DoubleDouble right = DoubleDouble.FromComponents(1.0, Math.ScaleB(1.0, -54)); + DoubleDouble expected = DoubleDouble.FromComponents(1.0, -Math.ScaleB(1.0, -54)); + AssertBits(DoubleDouble.Modulus(2.0, right), expected); + AssertBits(2.0 % right, expected); + AssertBits(new DoubleDouble(2.0) % right, expected); + } + + [Fact] + public void ScalarSpecialValuesAndRemainderSignsMatchBinary64() + { + double[] values = [0.0, -0.0, 1.0, -1.0, 2.0, -2.0, 5.5, -5.5, + double.Epsilon, -double.Epsilon, double.MaxValue, -double.MaxValue, + double.PositiveInfinity, double.NegativeInfinity, double.NaN]; + foreach (double left in values) + { + foreach (double right in values) + { + // The BCL is an independent oracle for zero-low inputs. NaNs + // are canonicalized by this library, regardless of the CPU payload. + DoubleDouble expected = new(left % right); + AssertBits(DoubleDouble.Modulus(new DoubleDouble(left), new DoubleDouble(right)), expected); + AssertBits(new DoubleDouble(left) % new DoubleDouble(right), expected); + AssertBits(new DoubleDouble(left) % right, expected); + AssertBits(left % new DoubleDouble(right), expected); + } + } + } + + [Fact] + public void ExactCasesIncludeSparseLowsCancellationAndUnrepresentableQuotients() + { + // Exact dyadic derivations; no expected result is calculated with DD arithmetic. + (DoubleDouble Left, DoubleDouble Right, DoubleDouble Expected)[] cases = + [ + // Tiny perturbations must decide the integer quotient before rounding. + (new(2.0, -double.Epsilon), new(1.0), new(1.0, -double.Epsilon)), + (new(2.0, double.Epsilon), new(1.0), new(double.Epsilon)), + (new(2.0), new(1.0, double.Epsilon), new(1.0, -double.Epsilon)), + (new(2.0), new(1.0, -double.Epsilon), new(2.0 * double.Epsilon)), + // Huge quotient exceeds both binary64 range and DD integer precision. + // MaxValue/epsilon = (2^53 - 1)*2^2045, congruent to 2 modulo 3. + (new(double.MaxValue, double.Epsilon), new(3.0 * double.Epsilon), new(0.0)), + (new(Math.ScaleB(1.0, 1000), 1.0), new(3.0), new(2.0)), + // Dividend smaller than the divisor, and an exact multiple with a low. + (new(1.0, -double.Epsilon), new(1.0), new(1.0, -double.Epsilon)), + (new(2.0, Math.ScaleB(1.0, -53)), new(1.0, Math.ScaleB(1.0, -54)), new(0.0)), + // Exact component sums satisfy left = 5*right, but High/High + // rounds below 5. The estimate must be corrected upward, not accepted. + (new(7.350727155000126, 2.4501029207264783e-16), new(1.4701454310000253, -8.422470454048922e-17), new(0.0)), + // Subnormal integer units: 5 % 3 = 2. + (new(5.0 * double.Epsilon), new(3.0 * double.Epsilon), new(2.0 * double.Epsilon)) + ]; + foreach ((DoubleDouble left, DoubleDouble right, DoubleDouble expected) in cases) + { + foreach (int leftSign in new[] { -1, 1 }) + { + foreach (int rightSign in new[] { -1, 1 }) + { + DoubleDouble a = leftSign < 0 ? -left : left; + DoubleDouble b = rightSign < 0 ? -right : right; + DoubleDouble result = leftSign < 0 ? -expected : expected; + AssertBits(a % b, result); + AssertBits(DoubleDouble.Modulus(a, b), result); + AssertBits(GenericModulus(a, b), result); + if (b.Low == 0.0) + { + AssertBits(a % b.High, result); + } + if (a.Low == 0.0) + { + AssertBits(a.High % b, result); + } + } + } + } + } + + [Fact] + public void NonfiniteDispatchDoesNotDiscardAFiniteDividendLow() + { + DoubleDouble value = new(1.0, double.Epsilon); + foreach (DoubleDouble left in new[] { value, -value }) + { + AssertBits(left % DoubleDouble.PositiveInfinity, left); + AssertBits(left % DoubleDouble.NegativeInfinity, left); + AssertBits(left % double.PositiveInfinity, left); + AssertBits(left % DoubleDouble.Zero, DoubleDouble.NaN); + AssertBits(left % DoubleDouble.NegativeZero, DoubleDouble.NaN); + AssertBits(left % DoubleDouble.NaN, DoubleDouble.NaN); + AssertBits(DoubleDouble.PositiveInfinity % left, DoubleDouble.NaN); + AssertBits(DoubleDouble.NaN % left, DoubleDouble.NaN); + } + } + + [Fact] + public void ModulusMeetsExactRationalBoundAcrossExponentsAndOverloads() + { + DoubleDouble[] values = FiniteCases(); + foreach (DoubleDouble left in values) + { + foreach (DoubleDouble right in values) + { + AssertRemainder(left % right, left, right); + AssertRemainder(left % right.High, left, new DoubleDouble(right.High)); + AssertRemainder(left.High % right, new DoubleDouble(left.High), right); + AssertBits(DoubleDouble.Modulus(left, right), left % right); + AssertBits(DoubleDouble.Modulus(left, right.High), left % right.High); + AssertBits(DoubleDouble.Modulus(left.High, right), left.High % right); + AssertBits(left % right.High, left % new DoubleDouble(right.High)); + AssertBits(left.High % right, new DoubleDouble(left.High) % right); + } + } + } + + [Fact] + public void ModulusSamplesEveryFiniteBinaryExponent() + { + for (int exponent = -1074; exponent <= 1023; ++exponent) + { + double high = Math.ScaleB(1.5, exponent); + DoubleDouble left = DoubleDouble.FromComponents(high, Math.ScaleB(high, -54)); + DoubleDouble right = DoubleDouble.FromComponents(3.0, -double.Epsilon); + AssertRemainder(left % right, left, right); + AssertRemainder(right % left, right, left); + AssertRemainder(-left % right, -left, right); + } + } + + [Theory] + [InlineData(0)] + [InlineData(1)] + [InlineData(2)] + public void OrdinaryTwoComponentModulusDoesNotAllocate(int overload) + { + // Include ordinary/dense operands, sparse quotient-boundary corrections, + // and common power-of-two rescalings; assertions/setup are not measured. + (DoubleDouble Left, DoubleDouble Right)[] cases = + [ + (new(5.5, Math.ScaleB(1.0, -54)), new(2.0, Math.ScaleB(1.0, -53))), + (new(2.0, -double.Epsilon), new(1.0, double.Epsilon)), + (new(2.0, double.Epsilon), new(1.0, -double.Epsilon)), + (new(2.0, Math.ScaleB(1.0, -53)), new(1.0, Math.ScaleB(1.0, -54))), + (new(7.350727155000126, 2.4501029207264783e-16), new(1.4701454310000253, -8.422470454048922e-17)), + (new(Math.ScaleB(5.5, -1000), Math.ScaleB(1.0, -1054)), new(Math.ScaleB(2.0, -1000), Math.ScaleB(1.0, -1053))), + (new(Math.ScaleB(5.5, 1000), Math.ScaleB(1.0, 946)), new(Math.ScaleB(2.0, 1000), Math.ScaleB(1.0, 947))), + // Identity/equality need no arithmetic even outside the kernel range. + (new(1.0, -double.Epsilon), new(double.MaxValue, double.Epsilon)), + (new(double.MaxValue, double.Epsilon), new(double.MaxValue, double.Epsilon)) + ]; + foreach ((DoubleDouble left, DoubleDouble right) in cases) + { + foreach (int sign in new[] { -1, 1 }) + { + DoubleDouble dividend = sign < 0 ? -left : left; + _ = MeasureModulusAllocations(overload, dividend, right, out _); + long allocated = MeasureModulusAllocations(overload, dividend, right, out DoubleDouble result); + AssertRemainder(result, overload == 2 ? new DoubleDouble(dividend.High) : dividend, + overload == 1 ? new DoubleDouble(right.High) : right); + allocated.ShouldBe(0L, $"overload {overload}: ({dividend.High:R}, {dividend.Low:R}) % ({right.High:R}, {right.Low:R})"); + } + } + } + + [Fact] + public void ModulusHandlesKernelRangeAndQuotientGuardNeighbors() + { + foreach (int exponent in new[] { -1000, -969, -54, 0, 450, 968, 1018, 1019 }) + { + DoubleDouble left = DoubleDouble.FromComponents(Math.ScaleB(5.5, exponent), Math.ScaleB(1.0, exponent - 54)); + DoubleDouble right = DoubleDouble.FromComponents(Math.ScaleB(2.0, exponent), -Math.ScaleB(1.0, exponent - 53)); + AssertRemainder(left % right, left, right); + AssertRemainder(left % right.High, left, new DoubleDouble(right.High)); + AssertRemainder(left.High % right, new DoubleDouble(left.High), right); + } + double limit = Math.ScaleB(1.0, 52); + foreach (double high in new[] { Math.BitDecrement(limit), limit, Math.BitIncrement(limit), limit + 1.0 }) + { + foreach (double low in new[] { -double.Epsilon, 0.0, double.Epsilon, -0.25, 0.25 }) + { + DoubleDouble left = DoubleDouble.FromComponents(high, low); + foreach (DoubleDouble right in new[] { new DoubleDouble(1.0, -double.Epsilon), + DoubleDouble.One, new DoubleDouble(1.0, double.Epsilon) }) + { + AssertRemainder(left % right, left, right); + AssertRemainder(-left % right, -left, right); + } + } + } + } + + [Fact] + public void RandomSignificandsMeetTheIndependentRemainderBound() + { + Random random = new(0x51a7); + int[] gaps = [0, 1, 2, 10, 51, 52, 53, 54, 500]; + for (int i = 0; i < 5000; ++i) + { + int exponent = random.Next(-1000, 501); + int gap = gaps[i % gaps.Length]; + DoubleDouble left = DoubleDouble.FromComponents(Math.ScaleB(1.0 + random.NextDouble(), exponent + gap), + Math.ScaleB(random.NextDouble() - 0.5, exponent + gap - 53)); + DoubleDouble right = DoubleDouble.FromComponents(Math.ScaleB(1.0 + random.NextDouble(), exponent), + Math.ScaleB(random.NextDouble() - 0.5, exponent - 53)); + AssertRemainder(left % right, left, right); + AssertRemainder(-left % -right, -left, -right); + AssertRemainder(left % right.High, left, new DoubleDouble(right.High)); + AssertRemainder(left.High % right, new DoubleDouble(left.High), right); + AssertBits(left % right.High, left % new DoubleDouble(right.High)); + AssertBits(left.High % right, new DoubleDouble(left.High) % right); + } + } + + [Fact] + public void UnsafeQuotientAndUpperRangeRetainTheExactAllocatingFallback() + { + (DoubleDouble Left, DoubleDouble Right)[] cases = + [ + (new(Math.ScaleB(1.0, 1000), 1.0), new(3.0)), + (new(Math.ScaleB(1.0, 1022), Math.ScaleB(1.0, 968)), new(Math.ScaleB(3.0, 1020), Math.ScaleB(1.0, 966))) + ]; + foreach ((DoubleDouble left, DoubleDouble right) in cases) + { + _ = MeasureModulusAllocations(0, left, right, out _); + long allocated = MeasureModulusAllocations(0, left, right, out DoubleDouble result); + allocated.ShouldBeGreaterThan(0L); + AssertRemainder(result, left, right); + } + } + + private static long MeasureModulusAllocations(int overload, DoubleDouble left, DoubleDouble right, out DoubleDouble result) + { + result = default; + long before = GC.GetAllocatedBytesForCurrentThread(); + for (int i = 0; i < 128; ++i) + { + result = overload switch + { + 0 => left % right, + 1 => left % right.High, + 2 => left.High % right, + _ => throw new ArgumentOutOfRangeException(nameof(overload)) + }; + } + return GC.GetAllocatedBytesForCurrentThread() - before; + } + + private static DoubleDouble[] FiniteCases() + { + List values = []; + foreach (int exponent in new[] { -1074, -1022, -969, -450, -54, 0, 54, 450, 969, 1023 }) + { + double high = Math.ScaleB(1.5, exponent); + foreach (double low in new[] { 0.0, Math.ScaleB(high, -54), -Math.ScaleB(high, -54), double.Epsilon, -double.Epsilon }) + { + DoubleDouble value = DoubleDouble.FromComponents(high, low); + if (!DoubleDouble.IsZero(value)) + { + values.Add(value); + values.Add(-value); + } + } + } + values.Add(DoubleDouble.FromComponents(double.MaxValue, Math.BitDecrement(Math.ScaleB(1.0, 970)))); + values.Add(new DoubleDouble(1.0, -Math.ScaleB(1.0, -54))); + values.Add(new DoubleDouble(2.0, -double.Epsilon)); + return values.ToArray(); + } + + private static void AssertRemainder(DoubleDouble actual, DoubleDouble left, DoubleDouble right) + { + // Independent exact dyadic oracle: extract significands with BCL scaling, + // then truncate the exact rational quotient as an integer. Never divide, + // multiply or subtract using the implementation under test. + BigInteger dividend = Units(left.High) + Units(left.Low); + BigInteger divisor = Units(right.High) + Units(right.Low); + BigInteger expected = dividend - ((dividend / divisor) * divisor); + DoubleDouble.IsFinite(actual).ShouldBeTrue(); + DoubleDouble.IsCanonical(actual).ShouldBeTrue(); + DoubleDouble.IsNegative(actual).ShouldBe(DoubleDouble.IsNegative(left)); + if (expected.IsZero) + { + AssertBits(actual, new DoubleDouble(Math.CopySign(0.0, left.High))); + return; + } + + actual.High.ShouldNotBe(0.0); + BigInteger observed = Units(actual.High) + Units(actual.Low); + BigInteger error = BigInteger.Abs(observed - expected); + // The usual arithmetic bound, in units of epsilon: relative 2^-100 + // plus one minimum subnormal. Keep both components in the comparison. + (error << 100).ShouldBeLessThanOrEqualTo(BigInteger.Abs(expected) + (BigInteger.One << 100)); + } + + private static BigInteger Units(double value) + { + if (value == 0.0) + { + return BigInteger.Zero; + } + int exponent = Math.ILogB(value); + BigInteger significand = new(Math.ScaleB(Math.Abs(value), -exponent) * 4503599627370496.0); + BigInteger units = significand << (exponent + 1022); + return double.IsNegative(value) ? -units : units; + } + + private static T GenericModulus(T left, T right) where T : INumber => left % right; + + private static T GenericScalarModulus(T left, double right) where T : IModulusOperators => left % right; + + private static void AssertBits(DoubleDouble actual, DoubleDouble expected) + { + 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/DoubleDoubleNumberTests.cs b/1-tests/Just.PreciseMath.Tests/DoubleDoubleNumberTests.cs new file mode 100644 index 0000000..f934977 --- /dev/null +++ b/1-tests/Just.PreciseMath.Tests/DoubleDoubleNumberTests.cs @@ -0,0 +1,99 @@ +namespace Just.PreciseMath.Tests; + +public class DoubleDoubleNumberTests +{ + [Fact] + public void NumberDefaultsUseBothComponents() + { + DoubleDouble lower = new(1.0, -double.Epsilon); + DoubleDouble upper = new(1.0, double.Epsilon); + AssertBits(Max(lower, upper), upper); + AssertBits(Max(upper, lower), upper); + AssertBits(MaxNumber(lower, upper), upper); + AssertBits(Min(lower, upper), lower); + AssertBits(Min(upper, lower), lower); + AssertBits(MinNumber(lower, upper), lower); + AssertBits(Clamp(DoubleDouble.Zero, lower, upper), lower); + AssertBits(Clamp(DoubleDouble.One, lower, upper), DoubleDouble.One); + AssertBits(Clamp(new DoubleDouble(2.0), lower, upper), upper); + Should.Throw(() => Clamp(DoubleDouble.One, upper, lower)); + AssertBits(CopySign(upper, DoubleDouble.NegativeZero), new DoubleDouble(-1.0, -double.Epsilon)); + AssertBits(CopySign(-upper, DoubleDouble.Zero), upper); + Sign(lower).ShouldBe(1); + Sign(-upper).ShouldBe(-1); + } + + [Fact] + public void NumberDefaultsHandleSpecialValuesAndSignedZeros() + { + double[] values = [0.0, -0.0, 1.0, -1.0, double.Epsilon, -double.Epsilon, + double.MaxValue, double.PositiveInfinity, double.NegativeInfinity, double.NaN]; + foreach (double left in values) + { + DoubleDouble x = new(left); + foreach (double right in values) + { + DoubleDouble y = new(right); + AssertBits(Max(x, y), new DoubleDouble(double.Max(left, right))); + AssertBits(Min(x, y), new DoubleDouble(double.Min(left, right))); + AssertBits(MaxNumber(x, y), new DoubleDouble(double.MaxNumber(left, right))); + AssertBits(MinNumber(x, y), new DoubleDouble(double.MinNumber(left, right))); + AssertBits(CopySign(x, y), new DoubleDouble(double.CopySign(left, right))); + } + AssertBits(Clamp(x, DoubleDouble.NegativeOne, DoubleDouble.One), new DoubleDouble(double.Clamp(left, -1.0, 1.0))); + if (double.IsNaN(left)) + { + Should.Throw(() => Sign(x)); + Should.Throw(() => DoubleDouble.Sign(x)); + } + else + { + Sign(x).ShouldBe(double.Sign(left)); + DoubleDouble.Sign(x).ShouldBe(double.Sign(left)); + } + } + } + +#if NET10_0_OR_GREATER + [Fact] + public void NativeDefaultsRetainFullComponentOrdering() + { + DoubleDouble lower = new(1.0, -double.Epsilon); + DoubleDouble upper = new(1.0, double.Epsilon); + AssertBits(MaxNative(lower, upper), upper); + AssertBits(MinNative(lower, upper), lower); + AssertBits(ClampNative(DoubleDouble.Zero, lower, upper), lower); + AssertBits(ClampNative(DoubleDouble.One, lower, upper), DoubleDouble.One); + AssertBits(ClampNative(new DoubleDouble(2.0), lower, upper), upper); + Should.Throw(() => ClampNative(DoubleDouble.One, upper, lower)); + } +#endif + + private static T Clamp(T value, T min, T max) where T : INumber => T.Clamp(value, min, max); + + private static T CopySign(T value, T sign) where T : INumber => T.CopySign(value, sign); + + private static T Max(T x, T y) where T : INumber => T.Max(x, y); + + private static T MaxNumber(T x, T y) where T : INumber => T.MaxNumber(x, y); + + private static T Min(T x, T y) where T : INumber => T.Min(x, y); + + private static T MinNumber(T x, T y) where T : INumber => T.MinNumber(x, y); + + private static int Sign(T value) where T : INumber => T.Sign(value); + +#if NET10_0_OR_GREATER + private static T ClampNative(T value, T min, T max) where T : INumber => T.ClampNative(value, min, max); + + private static T MaxNative(T x, T y) where T : INumber => T.MaxNative(x, y); + + private static T MinNative(T x, T y) where T : INumber => T.MinNative(x, y); +#endif + + private static void AssertBits(DoubleDouble actual, DoubleDouble expected) + { + BitConverter.DoubleToInt64Bits(actual.High).ShouldBe(BitConverter.DoubleToInt64Bits(expected.High)); + BitConverter.DoubleToInt64Bits(actual.Low).ShouldBe(BitConverter.DoubleToInt64Bits(expected.Low)); + } +} diff --git a/README.md b/README.md index 29b1888..8dcacbd 100644 --- a/README.md +++ b/README.md @@ -22,20 +22,37 @@ 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: static `Add`, `Subtract`, `Multiply`, and `Divide` methods accept +- Arithmetic: static `Add`, `Subtract`, `Multiply`, `Divide`, and `Modulus` methods accept `DoubleDouble` pairs and both operand orders with a `double`. `Plus` and `Negate` - provide unary operations. The arithmetic operators (`+`, `-`, `*`, `/`, unary + provide unary operations. The arithmetic operators (`+`, `-`, `*`, `/`, `%`, unary `+`/`-`, and `++`/`--`) are thin wrappers around these methods. Addition retains residuals under cancellation; multiplication uses - fused multiply-add; division uses residual corrections. Mixed `double` operators + fused multiply-add; division uses residual corrections. Mixed `double` overloads use specialized scalar paths rather than promoting the scalar to `DoubleDouble`. Division retains the complete first remainder and applies two quotient corrections, normalizing before the final correction. Mixed addition and subtraction also retain intermediate sum residuals through final normalization. Mixed addition, subtraction, and division produce the same component bits as their `DoubleDouble` operations with the scalar represented as a zero-low pair. +- Remainder: `Modulus` and `%` use C# truncation-toward-zero semantics, not + Euclidean modulo or IEEE nearest-integer remainder. Both components determine + the exact integer quotient; the remainder is then rounded to a normalized pair. + Its sign, including zero, follows the dividend. A zero divisor, infinite dividend, + or NaN operand produces canonical NaN; a finite dividend modulo infinity is + unchanged. Zero-low operands use binary64 `%`; smaller dividends and equal-high + operands also have non-allocating shortcuts. The general finite path uses FMA + residuals and stack-based error-free sums to verify the exact integer quotient + before rounding the remainder. Its initial quotient estimate is limited to + `2^52` and the dividend's high exponent to 1020, leaving safe transform headroom. + Outside those bounds, or if quotient correction cannot be certified, bounded + `BigInteger` arithmetic is the fallback. Warmed allocation tests pin zero managed + allocations for representative two-component cases in all three overloads; + no timing improvement is claimed. Mixed overloads retain the component bits of the + equivalent DD/DD operation. Exact-dyadic tests and independent rational checks + cover sparse lows, quotient boundaries, subnormals and extreme exponents; this + does not promise universally correctly rounded component pairs. - Exponent boundaries: bounded `BigInteger` calculations avoid intermediate - overflow and underflow on the exceptional finite path. Ordinary arithmetic uses + overflow and underflow on the exceptional finite path. Ordinary `+`, `-`, `*`, `/` use floating-point transforms without allocations, although sparse division correction products can also reach a boundary path. The stored value remains two doubles; this is not an arbitrary-precision API. @@ -389,13 +406,22 @@ sampled error bounds does not guarantee correct range decisions for every input. to recover the components. This is not a shortest-round-trip format and does not change the standard formats' 999-digit precision limit or the parser's length limit. -`DoubleDouble` implements `ISignedNumber`, including the inherited +`DoubleDouble` implements `INumber` and `ISignedNumber`, including the inherited `INumberBase` contracts: binary radix, classification, absolute value, magnitude selection, increment/decrement, and generic numeric conversions. Integer/parity tests and magnitude comparisons retain both components. Magnitude ties prefer positive values for maximum and negative values for minimum, including signed zero; the `Number` variants prefer a number over NaN. +`INumber` provides comparison and remainder dispatch as well as constrained generic +access to `Clamp`, `CopySign`, `Max`, `Min`, `MaxNumber`, and `MinNumber` through +the interface defaults. The .NET 10 native variants also use their defaults. +These defaults do not create directly callable static methods on `DoubleDouble`. +`DoubleDouble.Sign` is explicitly implemented: it returns -1, 0, or 1 and throws +`ArithmeticException` for NaN consistently on both target frameworks (unlike the +.NET 8 interface default). Comparisons and sign copying preserve low components +and signed-zero behavior; NaN results retain the library's canonical representation. + `CreateChecked`, `CreateSaturating`, and `CreateTruncating` support built-in numeric types and `BigInteger`. Floating overflow produces signed infinity in all modes. Finite integer output truncates the exact sum, then throws on overflow, clamps, @@ -459,7 +485,7 @@ bool success = DoubleDouble.TryParse("1.25e-2".AsSpan(), CultureInfo.InvariantCu The complete logarithmic and exponential families on `DoubleDouble`, and all three `DoubleDouble.Pow` overloads are implemented. -Generic-math interfaces beyond `ISignedNumber`, `IFloatingPointConstants`, +Generic-math interfaces beyond `INumber`, `ISignedNumber`, `IFloatingPointConstants`, `IRootFunctions`, `IExponentialFunctions`, `ILogarithmicFunctions`, and `IPowerFunctions`, additional text formats/shortest- round-trip formatting, and non-arithmetic performance benchmarks remain deferred. Replacing allocating arithmetic boundary fallbacks is also deferred; the current