diff --git a/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.AddSub.cs b/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.AddSub.cs index a151c983bd5d03..c088f325321982 100644 --- a/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.AddSub.cs +++ b/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.AddSub.cs @@ -71,13 +71,13 @@ public static void AddSelf(Span left, ReadOnlySpan right) // Same as above, but we're writing the result directly to a and // stop execution, if we're out of b and c is already 0. - for ( ; i < right.Length; i++) + for (; i < right.Length; i++) { long digit = (Unsafe.Add(ref leftPtr, i) + carry) + right[i]; Unsafe.Add(ref leftPtr, i) = unchecked((uint)digit); carry = digit >> 32; } - for ( ; carry != 0 && i < left.Length; i++) + for (; carry != 0 && i < left.Length; i++) { long digit = left[i] + carry; left[i] = (uint)digit; @@ -100,7 +100,7 @@ public static void Subtract(ReadOnlySpan left, ReadOnlySpan right, S { Debug.Assert(right.Length >= 1); Debug.Assert(left.Length >= right.Length); - Debug.Assert(Compare(left, right) >= 0); + Debug.Assert(CompareActual(left, right) >= 0); Debug.Assert(bits.Length == left.Length); // Switching to managed references helps eliminating @@ -132,8 +132,9 @@ public static void Subtract(ReadOnlySpan left, ReadOnlySpan right, S public static void SubtractSelf(Span left, ReadOnlySpan right) { Debug.Assert(left.Length >= right.Length); + // Assertion failing per https://github.com/dotnet/runtime/issues/97780 - // Debug.Assert(Compare(left, right) >= 0); + // Debug.Assert(CompareActual(left, right) >= 0); int i = 0; long carry = 0L; diff --git a/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.DivRem.cs b/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.DivRem.cs index bddeb75298f481..899691bf0ace91 100644 --- a/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.DivRem.cs +++ b/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.DivRem.cs @@ -8,35 +8,39 @@ namespace System.Numerics { internal static partial class BigIntegerCalculator { +#if DEBUG + // Mutable for unit testing... + internal static +#else + internal const +#endif + int DivideBurnikelZieglerThreshold = 32; + public static void Divide(ReadOnlySpan left, uint right, Span quotient, out uint remainder) { - Debug.Assert(left.Length >= 1); - Debug.Assert(quotient.Length == left.Length); - - // Executes the division for one big and one 32-bit integer. - // Thus, we've similar code than below, but there is no loop for - // processing the 32-bit integer, since it's a single element. - + InitializeForDebug(quotient); ulong carry = 0UL; - - for (int i = left.Length - 1; i >= 0; i--) - { - ulong value = (carry << 32) | left[i]; - ulong digit = value / right; - quotient[i] = (uint)digit; - carry = value - digit * right; - } + Divide(left, right, quotient, ref carry); remainder = (uint)carry; } public static void Divide(ReadOnlySpan left, uint right, Span quotient) + { + InitializeForDebug(quotient); + ulong carry = 0UL; + Divide(left, right, quotient, ref carry); + } + + private static void Divide(ReadOnlySpan left, uint right, Span quotient, ref ulong carry) { Debug.Assert(left.Length >= 1); Debug.Assert(quotient.Length == left.Length); + InitializeForDebug(quotient); - // Same as above, but only computing the quotient. + // Executes the division for one big and one 32-bit integer. + // Thus, we've similar code than below, but there is no loop for + // processing the 32-bit integer, since it's a single element. - ulong carry = 0UL; for (int i = left.Length - 1; i >= 0; i--) { ulong value = (carry << 32) | left[i]; @@ -68,9 +72,18 @@ public static void Divide(ReadOnlySpan left, ReadOnlySpan right, Spa Debug.Assert(left.Length >= right.Length); Debug.Assert(quotient.Length == left.Length - right.Length + 1); Debug.Assert(remainder.Length == left.Length); + InitializeForDebug(quotient); + InitializeForDebug(remainder); - left.CopyTo(remainder); - Divide(remainder, right, quotient); + if (right.Length < DivideBurnikelZieglerThreshold || left.Length - right.Length < DivideBurnikelZieglerThreshold) + { + left.CopyTo(remainder); + DivideGrammarSchool(remainder, right, quotient); + } + else + { + DivideBurnikelZiegler(left, right, quotient, remainder); + } } public static void Divide(ReadOnlySpan left, ReadOnlySpan right, Span quotient) @@ -79,22 +92,30 @@ public static void Divide(ReadOnlySpan left, ReadOnlySpan right, Spa Debug.Assert(right.Length >= 1); Debug.Assert(left.Length >= right.Length); Debug.Assert(quotient.Length == left.Length - right.Length + 1); + InitializeForDebug(quotient); - // Same as above, but only returning the quotient. + if (right.Length < DivideBurnikelZieglerThreshold || left.Length - right.Length < DivideBurnikelZieglerThreshold) + { + // Same as above, but only returning the quotient. - uint[]? leftCopyFromPool = null; + uint[]? leftCopyFromPool = null; - // NOTE: left will get overwritten, we need a local copy - // However, mutated left is not used afterwards, so use array pooling or stack alloc - Span leftCopy = (left.Length <= StackAllocThreshold ? - stackalloc uint[StackAllocThreshold] - : leftCopyFromPool = ArrayPool.Shared.Rent(left.Length)).Slice(0, left.Length); - left.CopyTo(leftCopy); + // NOTE: left will get overwritten, we need a local copy + // However, mutated left is not used afterwards, so use array pooling or stack alloc + Span leftCopy = (left.Length <= StackAllocThreshold ? + stackalloc uint[StackAllocThreshold] + : leftCopyFromPool = ArrayPool.Shared.Rent(left.Length)).Slice(0, left.Length); + left.CopyTo(leftCopy); - Divide(leftCopy, right, quotient); + DivideGrammarSchool(leftCopy, right, quotient); - if (leftCopyFromPool != null) - ArrayPool.Shared.Return(leftCopyFromPool); + if (leftCopyFromPool != null) + ArrayPool.Shared.Return(leftCopyFromPool); + } + else + { + DivideBurnikelZiegler(left, right, quotient, default); + } } public static void Remainder(ReadOnlySpan left, ReadOnlySpan right, Span remainder) @@ -102,21 +123,97 @@ public static void Remainder(ReadOnlySpan left, ReadOnlySpan right, Debug.Assert(left.Length >= 1); Debug.Assert(right.Length >= 1); Debug.Assert(left.Length >= right.Length); - Debug.Assert(remainder.Length >= left.Length); + Debug.Assert(remainder.Length == left.Length); + InitializeForDebug(remainder); - // Same as above, but only returning the remainder. + if (right.Length < DivideBurnikelZieglerThreshold || left.Length - right.Length < DivideBurnikelZieglerThreshold) + { + // Same as above, but only returning the remainder. - left.CopyTo(remainder); - Divide(remainder, right, default); + left.CopyTo(remainder); + DivideGrammarSchool(remainder, right, default); + } + else + { + int quotientLength = left.Length - right.Length + 1; + uint[]? quotientFromPool = null; + + Span quotient = (quotientLength <= StackAllocThreshold ? + stackalloc uint[StackAllocThreshold] + : quotientFromPool = ArrayPool.Shared.Rent(quotientLength)).Slice(0, quotientLength); + + DivideBurnikelZiegler(left, right, quotient, remainder); + + if (quotientFromPool != null) + ArrayPool.Shared.Return(quotientFromPool); + } } - private static void Divide(Span left, ReadOnlySpan right, Span bits) + /// + /// Logically equivalent to the following code. + /// + /// quotient = left / right; + /// left %= right; + /// + /// + private static void DivRem(Span left, ReadOnlySpan right, Span quotient) { Debug.Assert(left.Length >= 1); Debug.Assert(right.Length >= 1); Debug.Assert(left.Length >= right.Length); - Debug.Assert(bits.Length == left.Length - right.Length + 1 - || bits.Length == 0); + Debug.Assert(quotient.Length == left.Length - right.Length + 1 + || quotient.Length == 0); + InitializeForDebug(quotient); + + if (right.Length < DivideBurnikelZieglerThreshold || left.Length - right.Length < DivideBurnikelZieglerThreshold) + { + DivideGrammarSchool(left, right, quotient); + } + else + { + uint[]? leftCopyFromPool = null; + + // NOTE: left will get overwritten, we need a local copy + // However, mutated left is not used afterwards, so use array pooling or stack alloc + Span leftCopy = (left.Length <= StackAllocThreshold ? + stackalloc uint[StackAllocThreshold] + : leftCopyFromPool = ArrayPool.Shared.Rent(left.Length)).Slice(0, left.Length); + left.CopyTo(leftCopy); + + uint[]? quotientActualFromPool = null; + scoped Span quotientActual; + + if (quotient.Length == 0) + { + int quotientLength = left.Length - right.Length + 1; + + quotientActual = (quotientLength <= StackAllocThreshold ? + stackalloc uint[StackAllocThreshold] + : quotientActualFromPool = ArrayPool.Shared.Rent(quotientLength)).Slice(0, quotientLength); + } + else + { + quotientActual = quotient; + } + + DivideBurnikelZiegler(leftCopy, right, quotientActual, left); + + if (quotientActualFromPool != null) + ArrayPool.Shared.Return(quotientActualFromPool); + if (leftCopyFromPool != null) + ArrayPool.Shared.Return(leftCopyFromPool); + } + } + + private static void DivideGrammarSchool(Span left, ReadOnlySpan right, Span quotient) + { + Debug.Assert(left.Length >= 1); + Debug.Assert(right.Length >= 1); + Debug.Assert(left.Length >= right.Length); + Debug.Assert( + quotient.Length == 0 + || quotient.Length == left.Length - right.Length + 1 + || (CompareActual(left.Slice(left.Length - right.Length), right) < 0 && quotient.Length == left.Length - right.Length)); // Executes the "grammar-school" algorithm for computing q = a / b. // Before calculating q_i, we get more bits into the highest bit @@ -185,8 +282,8 @@ private static void Divide(Span left, ReadOnlySpan right, Span } // We have the digit! - if ((uint)n < (uint)bits.Length) - bits[n] = (uint)digit; + if ((uint)n < (uint)quotient.Length) + quotient[n] = (uint)digit; if ((uint)i < (uint)left.Length) left[i] = 0; @@ -230,7 +327,7 @@ private static uint SubtractDivisor(Span left, ReadOnlySpan right, u ref uint leftElement = ref left[i]; if (leftElement < digit) ++carry; - leftElement = unchecked(leftElement - digit); + leftElement -= digit; } return (uint)carry; @@ -250,19 +347,352 @@ private static bool DivideGuessTooBig(ulong q, ulong valHi, uint valLo, ulong chkLo = divLo * q; chkHi += (chkLo >> 32); - chkLo &= 0xFFFFFFFF; + uint chkLoUInt32 = (uint)(chkLo); + + return (chkHi > valHi) || ((chkHi == valHi) && (chkLoUInt32 > valLo)); + } + + private static void DivideBurnikelZiegler(ReadOnlySpan left, ReadOnlySpan right, Span quotient, Span remainder) + { + Debug.Assert(left.Length >= 1); + Debug.Assert(right.Length >= 1); + Debug.Assert(left.Length >= right.Length); + Debug.Assert(quotient.Length == left.Length - right.Length + 1); + Debug.Assert(remainder.Length == left.Length + || remainder.Length == 0); + + // Executes the Burnikel-Ziegler algorithm for computing q = a / b. + // + // Burnikel, C., Ziegler, J.: Fast recursive division. Research Report MPI-I-98-1-022, MPI Saarbrücken, 1998 + + + // Fast recursive division: Algorithm 3 + int n; + { + // m = min{1< right.Length} + int m = (int)BitOperations.RoundUpToPowerOf2((uint)right.Length / (uint)DivideBurnikelZieglerThreshold + 1); + + int j = (right.Length + m - 1) / m; // Ceil(right.Length/m) + n = j * m; + } + + int sigmaDigit = n - right.Length; + int sigmaSmall = BitOperations.LeadingZeroCount(right[^1]); + + uint[]? bFromPool = null; + + Span b = (n <= StackAllocThreshold ? + stackalloc uint[StackAllocThreshold] + : bFromPool = ArrayPool.Shared.Rent(n)).Slice(0, n); + + int aLength = left.Length + sigmaDigit; + + // if: BitOperations.LeadingZeroCount(left[^1]) < sigmaSmall, requires one more digit obviously. + // if: BitOperations.LeadingZeroCount(left[^1]) == sigmaSmall, requires one more digit, because the leftmost bit of a must be 0. + + if (BitOperations.LeadingZeroCount(left[^1]) <= sigmaSmall) + ++aLength; + + uint[]? aFromPool = null; + + Span a = (aLength <= StackAllocThreshold ? + stackalloc uint[StackAllocThreshold] + : aFromPool = ArrayPool.Shared.Rent(aLength)).Slice(0, aLength); + + // 4. normalize + static void Normalize(ReadOnlySpan src, int sigmaDigit, int sigmaSmall, Span bits) + { + Debug.Assert((uint)sigmaSmall <= 32); + Debug.Assert(src.Length + sigmaDigit <= bits.Length); + + bits.Slice(0, sigmaDigit).Clear(); + Span dst = bits.Slice(sigmaDigit); + src.CopyTo(dst); + dst.Slice(src.Length).Clear(); + + if (sigmaSmall != 0) + { + // Left shift + int carryShift = 32 - sigmaSmall; + uint carry = 0; + + for (int i = 0; i < bits.Length; i++) + { + uint carryTmp = bits[i] >> carryShift; + bits[i] = bits[i] << sigmaSmall | carry; + carry = carryTmp; + } - if (chkHi < valHi) - return false; - if (chkHi > valHi) - return true; + Debug.Assert(carry == 0); + } + } + + Normalize(left, sigmaDigit, sigmaSmall, a); + Normalize(right, sigmaDigit, sigmaSmall, b); + + + int t = Math.Max(2, (a.Length + n - 1) / n); // Max(2, Ceil(a.Length/n)) + Debug.Assert(t < a.Length || (t == a.Length && (int)a[^1] >= 0)); + + uint[]? rFromPool = null; + Span r = ((n + 1) <= StackAllocThreshold ? + stackalloc uint[StackAllocThreshold] + : rFromPool = ArrayPool.Shared.Rent(n + 1)).Slice(0, n + 1); + + uint[]? zFromPool = null; + Span z = (2 * n <= StackAllocThreshold ? + stackalloc uint[StackAllocThreshold] + : zFromPool = ArrayPool.Shared.Rent(2 * n)).Slice(0, 2 * n); + a.Slice((t - 2) * n).CopyTo(z); + z.Slice(a.Length - (t - 2) * n).Clear(); + + Span quotientUpper = quotient.Slice((t - 2) * n); + if (quotientUpper.Length < n) + { + uint[]? qFromPool = null; + Span q = (n <= StackAllocThreshold ? + stackalloc uint[StackAllocThreshold] + : qFromPool = ArrayPool.Shared.Rent(n)).Slice(0, n); + + BurnikelZieglerD2n1n(z, b, q, r); + + Debug.Assert(!q.Slice(quotientUpper.Length).ContainsAnyExcept(0u)); + q.Slice(0, quotientUpper.Length).CopyTo(quotientUpper); + + if (qFromPool != null) + ArrayPool.Shared.Return(qFromPool); + } + else + { + BurnikelZieglerD2n1n(z, b, quotientUpper.Slice(0, n), r); + quotientUpper.Slice(n).Clear(); + } + + for (int i = t - 3; i >= 0; i--) + { + a.Slice(i * n, n).CopyTo(z); + r.Slice(0, n).CopyTo(z.Slice(n)); + BurnikelZieglerD2n1n(z, b, quotient.Slice(i * n, n), r); + } + + if (zFromPool != null) + ArrayPool.Shared.Return(zFromPool); + if (bFromPool != null) + ArrayPool.Shared.Return(bFromPool); + if (aFromPool != null) + ArrayPool.Shared.Return(aFromPool); + + Debug.Assert(r[^1] == 0); + Debug.Assert(!r.Slice(0, sigmaDigit).ContainsAnyExcept(0u)); + if (remainder.Length != 0) + { + Span rt = r.Slice(sigmaDigit); + remainder.Slice(rt.Length).Clear(); + + if (sigmaSmall != 0) + { + // Right shift + Debug.Assert((uint)sigmaSmall <= 32); - if (chkLo < valLo) - return false; - if (chkLo > valLo) - return true; + int carryShift = 32 - sigmaSmall; + uint carry = 0; - return false; + for (int i = rt.Length - 1; i >= 0; i--) + { + remainder[i] = rt[i] >> sigmaSmall | carry; + carry = rt[i] << carryShift; + } + + Debug.Assert(carry == 0); + } + else + { + rt.CopyTo(remainder); + } + } + + if (rFromPool != null) + ArrayPool.Shared.Return(rFromPool); + } + + private static void BurnikelZieglerFallback(ReadOnlySpan left, ReadOnlySpan right, Span quotient, Span remainder) + { + // Fast recursive division: Algorithm 1 + // 1. If n is odd or smaller than some convenient constant + + Debug.Assert(left.Length == 2 * right.Length); + Debug.Assert(CompareActual(left.Slice(right.Length), right) < 0); + Debug.Assert(quotient.Length == right.Length); + Debug.Assert(remainder.Length >= right.Length + 1); + Debug.Assert(right[^1] > 0); + Debug.Assert(right.Length < DivideBurnikelZieglerThreshold); + + left = left.Slice(0, ActualLength(left)); + + if (left.Length < right.Length) + { + quotient.Clear(); + left.CopyTo(remainder); + remainder.Slice(left.Length).Clear(); + } + else if (right.Length == 1) + { + ulong carry; + + if (quotient.Length < left.Length) + { + Debug.Assert(quotient.Length + 1 == left.Length); + Debug.Assert(left[^1] < right[0]); + + carry = left[^1]; + Divide(left.Slice(0, quotient.Length), right[0], quotient, ref carry); + } + else + { + carry = 0; + quotient.Slice(left.Length).Clear(); + Divide(left, right[0], quotient, ref carry); + } + + if (remainder.Length != 0) + { + remainder.Slice(1).Clear(); + remainder[0] = (uint)carry; + } + } + else + { + uint[]? r1FromPool = null; + Span r1 = (left.Length <= StackAllocThreshold ? + stackalloc uint[StackAllocThreshold] + : r1FromPool = ArrayPool.Shared.Rent(left.Length)).Slice(0, left.Length); + + left.CopyTo(r1); + int quotientLength = Math.Min(left.Length - right.Length + 1, quotient.Length); + + quotient.Slice(quotientLength).Clear(); + DivideGrammarSchool(r1, right, quotient.Slice(0, quotientLength)); + + if (r1.Length < remainder.Length) + { + remainder.Slice(r1.Length).Clear(); + r1.CopyTo(remainder); + } + else + { + Debug.Assert(!r1.Slice(remainder.Length).ContainsAnyExcept(0u)); + r1.Slice(0, remainder.Length).CopyTo(remainder); + } + + if (r1FromPool != null) + ArrayPool.Shared.Return(r1FromPool); + } + } + + private static void BurnikelZieglerD2n1n(ReadOnlySpan left, ReadOnlySpan right, Span quotient, Span remainder) + { + // Fast recursive division: Algorithm 1 + Debug.Assert(left.Length == 2 * right.Length); + Debug.Assert(CompareActual(left.Slice(right.Length), right) < 0); + Debug.Assert(quotient.Length == right.Length); + Debug.Assert(remainder.Length >= right.Length + 1); + Debug.Assert(right[^1] > 0); + + if ((right.Length & 1) != 0 || right.Length < DivideBurnikelZieglerThreshold) + { + BurnikelZieglerFallback(left, right, quotient, remainder); + return; + } + + int halfN = right.Length >> 1; + + uint[]? r1FromPool = null; + Span r1 = ((right.Length + 1) <= StackAllocThreshold ? + stackalloc uint[StackAllocThreshold] + : r1FromPool = ArrayPool.Shared.Rent(right.Length + 1)).Slice(0, right.Length + 1); + + BurnikelZieglerD3n2n(left.Slice(right.Length), left.Slice(halfN, halfN), right, quotient.Slice(halfN), r1); + BurnikelZieglerD3n2n(r1.Slice(0, right.Length), left.Slice(0, halfN), right, quotient.Slice(0, halfN), remainder); + + if (r1FromPool != null) + ArrayPool.Shared.Return(r1FromPool); + } + + private static void BurnikelZieglerD3n2n(ReadOnlySpan left12, ReadOnlySpan left3, ReadOnlySpan right, Span quotient, Span remainder) + { + // Fast recursive division: Algorithm 2 + Debug.Assert(right.Length % 2 == 0); + Debug.Assert(left12.Length == right.Length); + Debug.Assert(2 * left3.Length == right.Length); + Debug.Assert(2 * quotient.Length == right.Length); + Debug.Assert(remainder.Length >= right.Length + 1); + Debug.Assert(right[^1] > 0); + Debug.Assert(CompareActual(left12, right) < 0); + + int n = right.Length >> 1; + + ReadOnlySpan a1 = left12.Slice(n); + ReadOnlySpan b1 = right.Slice(n); + ReadOnlySpan b2 = right.Slice(0, n); + Span r1 = remainder.Slice(n); + uint[]? dFromPool = null; + Span d = (right.Length <= StackAllocThreshold ? + stackalloc uint[StackAllocThreshold] + : dFromPool = ArrayPool.Shared.Rent(right.Length)).Slice(0, right.Length); + + if (CompareActual(a1, b1) < 0) + { + BurnikelZieglerD2n1n(left12, b1, quotient, r1); + + d.Clear(); + MultiplyActual(quotient, b2, d); + } + else + { + Debug.Assert(CompareActual(a1, b1) == 0); + quotient.Fill(uint.MaxValue); + + ReadOnlySpan a2 = left12.Slice(0, n); + Add(a2, b1, r1); + + d.Slice(0, n).Clear(); + b2.CopyTo(d.Slice(n)); + SubtractSelf(d, b2); + } + + // R = [R1, A3] + left3.CopyTo(remainder.Slice(0, n)); + + Span rr = remainder.Slice(0, d.Length + 1); + + while (CompareActual(rr, d) < 0) + { + AddSelf(rr, right); + int qi = -1; + while (quotient[++qi] == 0) ; + Debug.Assert((uint)qi < (uint)quotient.Length); + --quotient[qi]; + quotient.Slice(0, qi).Fill(uint.MaxValue); + } + + SubtractSelf(rr, d); + + if (dFromPool != null) + ArrayPool.Shared.Return(dFromPool); + + static void MultiplyActual(ReadOnlySpan left, ReadOnlySpan right, Span bits) + { + Debug.Assert(bits.Length == left.Length + right.Length); + + left = left.Slice(0, ActualLength(left)); + right = right.Slice(0, ActualLength(right)); + bits = bits.Slice(0, left.Length + right.Length); + + if (left.Length < right.Length) + Multiply(right, left, bits); + else + Multiply(left, right, bits); + } } } } diff --git a/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.FastReducer.cs b/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.FastReducer.cs index 0b56d2b52e2458..01acaf675b144b 100644 --- a/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.FastReducer.cs +++ b/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.FastReducer.cs @@ -32,7 +32,7 @@ public FastReducer(ReadOnlySpan modulus, Span r, Span mu, Span r[r.Length - 1] = 1; // Let mu = 4^k / m - Divide(r, modulus, mu); + DivRem(r, modulus, mu); _modulus = modulus; _q1 = q1; diff --git a/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.GcdInv.cs b/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.GcdInv.cs index 43995f88965a82..e17ff2e363e995 100644 --- a/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.GcdInv.cs +++ b/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.GcdInv.cs @@ -58,7 +58,7 @@ public static void Gcd(ReadOnlySpan left, ReadOnlySpan right, Span= 2); Debug.Assert(right.Length >= 2); - Debug.Assert(Compare(left, right) >= 0); + Debug.Assert(CompareActual(left, right) >= 0); Debug.Assert(result.Length == left.Length); left.CopyTo(result); diff --git a/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.PowMod.cs b/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.PowMod.cs index a0fc78df61e24a..3836fa08293005 100644 --- a/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.PowMod.cs +++ b/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.PowMod.cs @@ -227,7 +227,7 @@ stackalloc uint[StackAllocThreshold] if (value.Length > modulus.Length) { - Remainder(value, modulus, valueCopy); + Remainder(value, modulus, valueCopy.Slice(0, value.Length)); } else { @@ -276,7 +276,7 @@ stackalloc uint[StackAllocThreshold] if (value.Length > modulus.Length) { - Remainder(value, modulus, valueCopy); + Remainder(value, modulus, valueCopy.Slice(0, value.Length)); } else { diff --git a/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.Utils.cs b/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.Utils.cs index 2af543d113f26e..fc1b6b1011a81a 100644 --- a/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.Utils.cs +++ b/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.Utils.cs @@ -31,6 +31,26 @@ public static int Compare(ReadOnlySpan left, ReadOnlySpan right) return left[iv] < right[iv] ? -1 : 1; } + private static int CompareActual(ReadOnlySpan left, ReadOnlySpan right) + { + if (left.Length != right.Length) + { + if (left.Length < right.Length) + { + if (ActualLength(right.Slice(left.Length)) > 0) + return -1; + right = right.Slice(0, left.Length); + } + else + { + if (ActualLength(left.Slice(right.Length)) > 0) + return +1; + left = left.Slice(0, right.Length); + } + } + return Compare(left, right); + } + public static int ActualLength(ReadOnlySpan value) { // Since we're reusing memory here, the actual length @@ -49,11 +69,18 @@ private static int Reduce(Span bits, ReadOnlySpan modulus) if (bits.Length >= modulus.Length) { - Divide(bits, modulus, default); + DivRem(bits, modulus, default); return ActualLength(bits.Slice(0, modulus.Length)); } return bits.Length; } + + [Conditional("DEBUG")] + public static void InitializeForDebug(Span bits) + { + // Reproduce the case where the return value of `stackalloc uint` is not initialized to zero. + bits.Fill(0xCD); + } } } diff --git a/src/libraries/System.Runtime.Numerics/tests/BigInteger/divide.cs b/src/libraries/System.Runtime.Numerics/tests/BigInteger/divide.cs index 92cd28c39f7ed9..8b2bdce1348d1e 100644 --- a/src/libraries/System.Runtime.Numerics/tests/BigInteger/divide.cs +++ b/src/libraries/System.Runtime.Numerics/tests/BigInteger/divide.cs @@ -25,6 +25,22 @@ public static void RunDivideTwoLargeBI() } } + [Fact] + public static void RunDivideOneLargeOneHalfBI() + { + byte[] tempByteArray1 = new byte[0]; + byte[] tempByteArray2 = new byte[0]; + + // Divide Method - Two Large BigIntegers + for (int i = -1; i <= 1; i++) + for (int j = -1; j <= 1; j++) + { + tempByteArray1 = GetRandomByteArray(s_random, 512 + i); + tempByteArray2 = GetRandomByteArray(s_random, 256 + j); + VerifyDivideString(Print(tempByteArray1) + Print(tempByteArray2) + "bDivide"); + } + } + [Fact] public static void RunDivideTwoSmallBI() { @@ -149,6 +165,31 @@ public static void RunDivideBoundary() VerifyDivideString(Math.Pow(2, 33) + " 2 bDivide"); } + [Fact] + public static void RunDivideLarge() + { + byte[] tempByteArray1 = new byte[0]; + byte[] tempByteArray2 = new byte[0]; + + for (int i = 1; i < 15; i += 3) + { + var num = 1 << i; + tempByteArray1 = GetRandomByteArray(s_random, num); + tempByteArray2 = GetRandomByteArray(s_random, num / 2); + VerifyDivideString(Print(tempByteArray1) + Print(tempByteArray2) + "bDivide"); + } + + // Divide Method - Two Large BigIntegers + for (int i = -1; i <= 1; i++) + for (int j = -1; j <= 1; j++) + { + int num = (1 << 13) + i; + tempByteArray1 = GetRandomByteArray(s_random, num); + tempByteArray2 = GetRandomByteArray(s_random, num / 3 + j); + VerifyDivideString(Print(tempByteArray1) + Print(tempByteArray2) + "bDivide"); + } + } + [Fact] public static void RunOverflow() { @@ -161,6 +202,21 @@ public static void RunOverflow() Assert.Equal(z, BigInteger.Divide(x, y)); } + [Fact] + public void D3n2nBound() + { + var right = (BigInteger.One << (BigIntegerCalculator.DivideBurnikelZieglerThreshold * 4 * 32 - 1)) + + (BigInteger.One << (BigIntegerCalculator.DivideBurnikelZieglerThreshold * 2 * 32)) - 1; + var rem = right - 1; + + var qi = BigIntegerCalculator.DivideBurnikelZieglerThreshold * 8 * 32 * 4 - 1; + var q = (BigInteger.One << qi) - 1; + var left = q * right + rem; + + var q2 = BigInteger.Divide(left, right); + Assert.Equal(q, q2); + } + private static void VerifyDivideString(string opstring) { StackCalc sc = new StackCalc(opstring); diff --git a/src/libraries/System.Runtime.Numerics/tests/BigInteger/divrem.cs b/src/libraries/System.Runtime.Numerics/tests/BigInteger/divrem.cs index 2b740e1efabeaa..14c08f0aa9aea9 100644 --- a/src/libraries/System.Runtime.Numerics/tests/BigInteger/divrem.cs +++ b/src/libraries/System.Runtime.Numerics/tests/BigInteger/divrem.cs @@ -59,6 +59,22 @@ public static void RunDivRem_OneSmallOneLargeBI() } } + [Fact] + public static void RunDivRem_OneLargeOneHalfBI() + { + byte[] tempByteArray1 = new byte[0]; + byte[] tempByteArray2 = new byte[0]; + + // Divide Method - Two Large BigIntegers + for (int i = -1; i <= 1; i++) + for (int j = -1; j <= 1; j++) + { + tempByteArray1 = GetRandomByteArray(s_random, 512 + i); + tempByteArray2 = GetRandomByteArray(s_random, 256 + j); + VerifyDivRemString(Print(tempByteArray1) + Print(tempByteArray2) + "bDivRem"); + } + } + [Fact] public static void RunDivRem_OneLargeOne0BI() { @@ -114,6 +130,73 @@ public static void Boundary() VerifyDivRemString(Math.Pow(2, 33) + " 2 bDivRem"); } + [Fact] + public static void RunDivRemMedium() + { + byte[] tempByteArray1 = new byte[0]; + byte[] tempByteArray2 = new byte[0]; + + for (int i = 1; i < 8; i++) + { + var num = 1 << i; + tempByteArray1 = GetRandomByteArray(s_random, num); + tempByteArray2 = GetRandomByteArray(s_random, num / 2); + VerifyDivRemString(Print(tempByteArray2) + Print(tempByteArray1) + "bDivRem"); + } + + // Divide Method - Two Large BigIntegers + for (int i = -1; i <= 1; i++) + for (int j = -1; j <= 1; j++) + { + int num = (1 << 7) + i; + tempByteArray1 = GetRandomByteArray(s_random, num); + tempByteArray2 = GetRandomByteArray(s_random, num / 3 + j); + VerifyDivRemString(Print(tempByteArray2) + Print(tempByteArray1) + "bDivRem"); + } + } + + [Fact] + [OuterLoop] + public static void RunDivRemLarge() + { + byte[] tempByteArray1 = new byte[0]; + byte[] tempByteArray2 = new byte[0]; + + for (int i = 8; i < 10; i++) + { + var num = 1 << i; + tempByteArray1 = GetRandomByteArray(s_random, num); + tempByteArray2 = GetRandomByteArray(s_random, num / 2); + VerifyDivRemString(Print(tempByteArray2) + Print(tempByteArray1) + "bDivRem"); + } + + // Divide Method - Two Large BigIntegers + for (int i = -1; i <= 1; i++) + for (int j = -1; j <= 1; j++) + { + int num = (1 << 10) + i; + tempByteArray1 = GetRandomByteArray(s_random, num); + tempByteArray2 = GetRandomByteArray(s_random, num / 3 + j); + VerifyDivRemString(Print(tempByteArray2) + Print(tempByteArray1) + "bDivRem"); + } + } + + [Fact] + public void D3n2nBound() + { + var right = (BigInteger.One << (BigIntegerCalculator.DivideBurnikelZieglerThreshold * 4 * 32 - 1)) + + (BigInteger.One << (BigIntegerCalculator.DivideBurnikelZieglerThreshold * 2 * 32)) - 1; + var rem = right - 1; + + var qi = BigIntegerCalculator.DivideBurnikelZieglerThreshold * 8 * 32 * 4 - 1; + var q = (BigInteger.One << qi) - 1; + var left = q * right + rem; + + var (q2, r2) = BigInteger.DivRem(left, right); + Assert.Equal(q, q2); + Assert.Equal(rem, r2); + } + [Fact] public static void RunDivRemTests() { @@ -211,4 +294,16 @@ private static string Print(byte[] bytes) return MyBigIntImp.Print(bytes); } } + + + [Collection(nameof(DisableParallelization))] + public class divremTestThreshold + { + [Fact] + public void RunDivRemMedium() + { + BigIntTools.Utils.RunWithFakeThreshold(BigIntegerCalculator.DivideBurnikelZieglerThreshold, 8, divremTest.RunDivRemMedium); + } + + } } diff --git a/src/libraries/System.Runtime.Numerics/tests/BigInteger/modpow.cs b/src/libraries/System.Runtime.Numerics/tests/BigInteger/modpow.cs index f2296da7583424..5f840a9a520990 100644 --- a/src/libraries/System.Runtime.Numerics/tests/BigInteger/modpow.cs +++ b/src/libraries/System.Runtime.Numerics/tests/BigInteger/modpow.cs @@ -156,6 +156,31 @@ public static void ModPow2Large1SmallInt() } } + // InlineData randomly generated using a new Random(0) and the same logic as is used in MyBigIntImp + // When using the VerifyModPowString approach, these tests were taking over 100s to execute. + [Theory] + [OuterLoop] + [InlineData("16152447169612934532253858872860101823992426489763666485380167221632010974596350526903901586786157942589749704137347759619657187784260082814484793173551647870373927196771852", "81750440863317178198977948957223170296434385667401518194303765341584160931262254891905182291359611229458113456891057547210603179481274334549407149567413965138338569615186383", "-7196805437625531929619339290119232497987723600413951949478998858079935202962418871178150151495255147489121599998392939293913217145771279946557062800866475934948778", "5616380469802854111883204613624062553120851412450654399190717359559254034938493684003207209263083893212739965150170342564431867389596229480294680122528026782956958")] + [InlineData("-19804882257271871615998867413310154569092691610863183900352723872936597192840601413521317416195627481021453523098672122338092864102790996172", "5005126326020660905045117276428380928616741884628331700680598638081087432348216495484429178211470872172882902036752474804904132643", "-103139435102033366056363338063498713727200190198435", "-4559679593639479333979462122256288215243720737403")] + [InlineData("-2263742881720266366742488789295767051326353176981812961528848101545308514727342989711034207353102864275894987436688194776201313772226059035143797121004935923223042331190890429211181925543749113890129660515170733848725", "313166794469483345944045915670847620773229708424183728974404367769353433443313247319899209581239311758905801464781585268041623664425", "3984523819293796257818294503433333109083365267654089093971156331037874112339941681299823291492643701164442964256327008451060278818307268485187524722068240", "-1502346344475556570236985614111713440763176490652894928852811056060906839905964408583918853958610145894621840382970373570196361549098246030413222124498085")] + [InlineData("-342701069914551895507647608205424602441707687704978754182486401529587201154652163656221404036731562708712821963465111719289193200432064875769386559645346920", "247989781302056934762335026151076445832166884867735502165354252207668083157132317377069604102049233014316799294014890817943261246892765157594412634897440785353202366563028", "121555428622583377664148832113619145387775383377032217138159544127299380518157949963314283123957062240152711187509503414343", "87578369862034238407481381238856926729112161247036763388287150463197193460326629595765471822752579542310337770783772458710")] + [InlineData("-282593950368131775131360433563237877977879464854725217398276355042086422366377452969365517205334520940629323311057746859", "5959258935361466196732788139899933762310441595693546573755448590100882267274831199165902682626769861372370905838130200967035", "6598019436100687108279703832125132357070343951841815860721553173215685978621505459125000339496396952653080051757197617932586296524960251609958919233462530", "-4035534917213670830138661334307123840766321272945538719146336835321463334263841831803093764143165039453996360590495570418141622764990299884920213157241339")] + [InlineData("-283588760164723199492097446398514238808996680579814508529658395835708678571214111525627362048679162021949222615325057958783388520044013697149834530411072380435126870273057157745943859", "1042611904427950637176782337251399378305726352782300537151930702574076654415735617544217054055762002017753284951033975382044655538090809873113604", "11173562248848713955826639654969725554069867375462328112124015145073186375237811117727665778232780449525476829715663254692620996726130822561707626585790443774041565237684936409844925424596571418502946", "6662129352542591544500713232459850446949913817909326081510506555775206102887692404112984526760120407457772674917650956873499346517965094865621779695963015030158124625116211104048940313058280680420919")] + [InlineData("683399436848090272429540515929404372035986606104192913128573936597145312908480564700882440819526500604918037963508399983297602351856208190565527", "223751878996658298590720798129833935741775535718632942242965085592936259576946732440406302671204348239437556817289012497483482656", "1204888420642360606571457515385663952017382888522547766961071698778167608427220546474854934298311882921224875807375321847274969073309433075566441363244101914489524538123410250010519308563731930923389473186", "1136631484875680074951300738594593722168933395227758228596156355418704717505715681950525129323209331607163560404958604424924870345828742295978850144844693079191828839673460389985036424301333691983679427765")] + [InlineData("736513799451968530811005754031332418210960966881742655756522735504778110620671049112529346250333710388060811959329786494662578020803", "2461175085563866950903873687720858523536520498137697316698237108626602445202960480677695918813575265778826908481129155012799", "-4722693720735888562993277045098354134891725536023070176847814685098361292027040929352405620815883795027263132404351040", "4351573186631261607388198896754285562669240685903971199359912143458682154189588696264319780329366022294935204028039787")] + [InlineData("1596188639947986471148999794547338", "685242191738212089917782567856594513073397739443", "41848166029740752457613562518205823134173790454761036532025758411484449588176128053901271638836032557551179866133091058357374964041641117585422447497779410336188602585660372002644517668041207657383104649333118253", "39246949850380693159338034407642149926180988060650630387722725303281343126585456713282439764667310808891687831648451269002447916277601468727040185218264602698691553232132525542650964722093335105211816394635493987")] + [InlineData("-1506852741293695463963822334869441845197951776565891060639754936248401744065969556756496718308248025911048010080232290368562210204958094544173209793990218122", "64905085725614938357105826012272472070910693443851911667721848542473785070747281964799379996923261457185847459711", "2740467233603031668807697475486217767705051", "-1905434239471820365929630558127219204166613")] + public static void ModPow3LargeInt(string value, string exponent, string modulus, string expected) + { + BigInteger valueInt = BigInteger.Parse(value); + BigInteger exponentInt = BigInteger.Parse(exponent); + BigInteger modulusInt = BigInteger.Parse(modulus); + BigInteger resultInt = BigInteger.Parse(expected); + + // Once with default threshold + Assert.Equal(resultInt, BigInteger.ModPow(valueInt, exponentInt, modulusInt)); + } + [Fact] public static void ModPow0Power() { diff --git a/src/libraries/System.Runtime.Numerics/tests/BigInteger/op_divide.cs b/src/libraries/System.Runtime.Numerics/tests/BigInteger/op_divide.cs index 5d5714c306d7e0..89ba3ce10b6be7 100644 --- a/src/libraries/System.Runtime.Numerics/tests/BigInteger/op_divide.cs +++ b/src/libraries/System.Runtime.Numerics/tests/BigInteger/op_divide.cs @@ -32,6 +32,15 @@ public static void RunDividePositivenonZero() VerifyDivideString(Print(tempByteArray1) + Print(tempByteArray2) + "b/"); } + // Divide Method - One large and one half BigIntegers + for (int i = -1; i <= 1; i++) + for (int j = -1; j <= 1; j++) + { + tempByteArray1 = GetRandomByteArray(s_random, 512 + i); + tempByteArray2 = GetRandomByteArray(s_random, 256 + j); + VerifyDivideString(Print(tempByteArray1) + Print(tempByteArray2) + "b/"); + } + // Divide Method - One large and one small BigIntegers for (int i = 0; i < s_samples; i++) { @@ -140,6 +149,21 @@ public static void RunOverflow() Assert.Equal(z, x / y); } + [Fact] + public void D3n2nBound() + { + var right = (BigInteger.One << (BigIntegerCalculator.DivideBurnikelZieglerThreshold * 4 * 32 - 1)) + + (BigInteger.One << (BigIntegerCalculator.DivideBurnikelZieglerThreshold * 2 * 32)) - 1; + var rem = right - 1; + + var qi = BigIntegerCalculator.DivideBurnikelZieglerThreshold * 8 * 32 * 4 - 1; + var q = (BigInteger.One << qi) - 1; + var left = q * right + rem; + + var q2 = left / right; + Assert.Equal(q, q2); + } + private static void VerifyDivideString(string opstring) { StackCalc sc = new StackCalc(opstring); diff --git a/src/libraries/System.Runtime.Numerics/tests/BigInteger/op_modulus.cs b/src/libraries/System.Runtime.Numerics/tests/BigInteger/op_modulus.cs index 97f561dce0577f..3481e2414c976e 100644 --- a/src/libraries/System.Runtime.Numerics/tests/BigInteger/op_modulus.cs +++ b/src/libraries/System.Runtime.Numerics/tests/BigInteger/op_modulus.cs @@ -33,6 +33,15 @@ public static void RunRemainderTestsPositive() VerifyRemainderString(Print(tempByteArray1) + Print(tempByteArray2) + "b%"); } + // Divide Method - One large and one half BigIntegers + for (int i = -1; i <= 1; i++) + for (int j = -1; j <= 1; j++) + { + tempByteArray1 = GetRandomByteArray(s_random, 512 + i); + tempByteArray2 = GetRandomByteArray(s_random, 256 + j); + VerifyRemainderString(Print(tempByteArray1) + Print(tempByteArray2) + "b%"); + } + // Remainder Method - One large and one small BigIntegers for (int i = 0; i < s_samples; i++) { @@ -168,6 +177,21 @@ public static void RunRemainderAxiom0ModX() } } + [Fact] + public void D3n2nBound() + { + var right = (BigInteger.One << (BigIntegerCalculator.DivideBurnikelZieglerThreshold * 4 * 32 - 1)) + + (BigInteger.One << (BigIntegerCalculator.DivideBurnikelZieglerThreshold * 2 * 32)) - 1; + var rem = right - 1; + + var qi = BigIntegerCalculator.DivideBurnikelZieglerThreshold * 8 * 32 * 4 - 1; + var q = (BigInteger.One << qi) - 1; + var left = q * right + rem; + + var r2 = left % right; + Assert.Equal(rem, r2); + } + private static void VerifyRemainderString(string opstring) { StackCalc sc = new StackCalc(opstring); diff --git a/src/libraries/System.Runtime.Numerics/tests/BigInteger/remainder.cs b/src/libraries/System.Runtime.Numerics/tests/BigInteger/remainder.cs index 9068e32dcfe0eb..902a54f90ef709 100644 --- a/src/libraries/System.Runtime.Numerics/tests/BigInteger/remainder.cs +++ b/src/libraries/System.Runtime.Numerics/tests/BigInteger/remainder.cs @@ -33,6 +33,15 @@ public static void RunRemainderPositive() VerifyRemainderString(Print(tempByteArray1) + Print(tempByteArray2) + "bRemainder"); } + // Divide Method - One large and one half BigIntegers + for (int i = -1; i <= 1; i++) + for (int j = -1; j <= 1; j++) + { + tempByteArray1 = GetRandomByteArray(s_random, 512 + i); + tempByteArray2 = GetRandomByteArray(s_random, 256 + j); + VerifyRemainderString(Print(tempByteArray1) + Print(tempByteArray2) + "bRemainder"); + } + // Remainder Method - One large and one small BigIntegers for (int i = 0; i < s_samples; i++) {