modulus implementation

0-source/Just.PreciseMath/DoubleDouble.Arithmetic.cs

@@ -4,7 +4,8 @@ public readonly partial struct DoubleDouble :
     IAdditionOperators<DoubleDouble, double, DoubleDouble>,
     ISubtractionOperators<DoubleDouble, double, DoubleDouble>,
     IMultiplyOperators<DoubleDouble, double, DoubleDouble>,
-    IDivisionOperators<DoubleDouble, double, DoubleDouble>
+    IDivisionOperators<DoubleDouble, double, DoubleDouble>,
+    IModulusOperators<DoubleDouble, double, DoubleDouble>
 {
     /// <summary>Returns the square using a specialized double-double product.</summary>
     /// <remarks>
@@ -152,6 +153,29 @@ public readonly partial struct DoubleDouble :
         return PreciseMathHelper.AddScalarFinite(partial._high, partial._low, finalCorrection);
     }
 
+    /// <summary>Returns the remainder after division with the quotient truncated toward zero.</summary>
+    /// <remarks>
+    /// 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.
+    /// </remarks>
+    [Pure]
+    public static DoubleDouble Modulus(DoubleDouble left, DoubleDouble right) => ModulusCore(left._high, left._low, right._high, right._low);
+
+    /// <inheritdoc cref="Modulus(DoubleDouble, DoubleDouble)"/>
+    [Pure]
+    public static DoubleDouble Modulus(DoubleDouble left, double right) => ModulusCore(left._high, left._low, right, 0.0);
+
+    /// <inheritdoc cref="Modulus(DoubleDouble, DoubleDouble)"/>
+    [Pure]
+    public static DoubleDouble Modulus(double left, DoubleDouble right) => ModulusCore(left, 0.0, right._high, right._low);
+
     /// <summary>Applies the expansion operation without discarding the low component.</summary>
     /// <remarks>Preserves the component bits of addition with the equivalent zero-low scalar expansion.</remarks>
     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<double> expansion = stackalloc double[8];
+        Span<double> 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<double> 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);
+    }
 }

0-source/Just.PreciseMath/DoubleDouble.NumberBase.cs

@@ -42,6 +42,12 @@ public readonly partial struct DoubleDouble
         return IsNegative(value) ? -value : value;
     }
 
+    /// <summary>Returns -1 for negative values, 0 for either zero, and 1 for positive values.</summary>
+    /// <exception cref="ArithmeticException">The value is NaN.</exception>
+    /// <remarks>Normalization makes the high component's sign authoritative. Overrides runtime-dependent INumber defaults.</remarks>
+    [Pure]
+    public static int Sign(DoubleDouble value) => Math.Sign(value._high);
+
     /// <summary>Returns the operand with greater magnitude; ties prefer positive values and NaN propagates.</summary>
     public static DoubleDouble MaxMagnitude(DoubleDouble x, DoubleDouble y)
     {

0-source/Just.PreciseMath/DoubleDouble.Operators.cs

@@ -20,6 +20,9 @@ public readonly partial struct DoubleDouble
     /// <inheritdoc cref="Divide(DoubleDouble, DoubleDouble)"/>
     public static DoubleDouble operator /(DoubleDouble left, DoubleDouble right) => Divide(left, right);
 
+    /// <inheritdoc cref="Modulus(DoubleDouble, DoubleDouble)"/>
+    public static DoubleDouble operator %(DoubleDouble left, DoubleDouble right) => Modulus(left, right);
+
     /// <inheritdoc cref="Add(DoubleDouble, double)"/>
     public static DoubleDouble operator +(DoubleDouble left, double right) => Add(left, right);
 
@@ -44,6 +47,12 @@ public readonly partial struct DoubleDouble
     /// <inheritdoc cref="Divide(double, DoubleDouble)"/>
     public static DoubleDouble operator /(double left, DoubleDouble right) => Divide(left, right);
 
+    /// <inheritdoc cref="Modulus(DoubleDouble, double)"/>
+    public static DoubleDouble operator %(DoubleDouble left, double right) => Modulus(left, right);
+
+    /// <inheritdoc cref="Modulus(double, DoubleDouble)"/>
+    public static DoubleDouble operator %(double left, DoubleDouble right) => Modulus(left, right);
+
     /// <summary>Adds one using double-double arithmetic, including its overflow and nonfinite behavior.</summary>
     public static DoubleDouble operator ++(DoubleDouble value) => Add(value, 1.0);
 

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.
 /// </remarks>
 public readonly partial struct DoubleDouble :
-    ISignedNumber<DoubleDouble>
+    ISignedNumber<DoubleDouble>, INumber<DoubleDouble>
 {
     internal readonly double _high;
     internal readonly double _low;