diff --git a/src/libraries/System.Private.CoreLib/src/System.Private.CoreLib.Shared.projitems b/src/libraries/System.Private.CoreLib/src/System.Private.CoreLib.Shared.projitems
index d297affcabc559..bbac4e6d58572d 100644
--- a/src/libraries/System.Private.CoreLib/src/System.Private.CoreLib.Shared.projitems
+++ b/src/libraries/System.Private.CoreLib/src/System.Private.CoreLib.Shared.projitems
@@ -633,6 +633,7 @@
+
diff --git a/src/libraries/System.Private.CoreLib/src/System/Number.BigInteger.cs b/src/libraries/System.Private.CoreLib/src/System/Number.BigInteger.cs
index 9b043a01c2842f..e117f1a0f99720 100644
--- a/src/libraries/System.Private.CoreLib/src/System/Number.BigInteger.cs
+++ b/src/libraries/System.Private.CoreLib/src/System/Number.BigInteger.cs
@@ -10,7 +10,7 @@ namespace System
{
internal static partial class Number
{
- [StructLayout(LayoutKind.Sequential, Pack = 1)]
+ [StructLayout(LayoutKind.Sequential)]
internal ref struct BigInteger
{
// The longest binary mantissa requires: explicit mantissa bits + abs(min exponent)
@@ -30,26 +30,192 @@ internal ref struct BigInteger
// We require BitsPerBlock additional bits for shift space used during the pre-division preparation
private const int MaxBits = BitsForLongestBinaryMantissa + BitsForLongestDigitSequence + BitsPerBlock;
- private const int BitsPerBlock = sizeof(int) * 8;
+ // Blocks are the native word width so the algorithms scale with the machine.
+#if TARGET_64BIT
+ internal const int BitsPerBlock = 64;
+ private const int BlockShift = 6;
+#else
+ internal const int BitsPerBlock = 32;
+ private const int BlockShift = 5;
+#endif
// We need one extra block to make our shift left algorithm significantly simpler
private const int MaxBlockCount = ((MaxBits + (BitsPerBlock - 1)) / BitsPerBlock) + 1;
- private static ReadOnlySpan Pow10UInt32Table =>
+ // Single-block powers of 10, shared with the primitive types. On 64-bit a block holds
+ // 10^0..10^19; on 32-bit, 10^0..10^9. MultiplyPow10 folds any exponent that fits a
+ // single block through here, deferring only larger exponents to Pow10BigNumTable.
+#if TARGET_64BIT
+ private static ReadOnlySpan Pow10Table => ulong.PowersOf10;
+#else
+ private static ReadOnlySpan Pow10Table => uint.PowersOf10;
+#endif
+
+ // Pow10 splits the exponent into a low part applied via a single Pow10Table lookup+multiply
+ // and a high part accumulated from Pow10BigNumTable. The low part covers 10^0..10^(2^N - 1),
+ // the largest power-of-two range whose values each fit in a single block: N == BlockShift - 2
+ // (3 on 32-bit -> Pow10BigNumTable starts at 10^8; 4 on 64-bit -> it starts at 10^16).
+ private const int Pow10SmallExpBits = BlockShift - 2;
+ private const uint Pow10SmallExpMask = (1u << Pow10SmallExpBits) - 1;
+
+ // Pow10BigNumTable stores each power as a valid BigInteger memory image: a `_length` (in
+ // blocks) immediately followed by that many packed blocks. Each bitness uses a blob whose
+ // element type matches the native word -- ulong on 64-bit, uint on 32-bit -- so both `_length`
+ // (a nint) and each block occupy exactly one element. That keeps the blocks naturally aligned
+ // and, since CreateSpan byte-swaps per element, endian-correct when reinterpreted as a
+ // BigInteger. The 64-bit variant is generated + verified by src\libraries tooling in history.
+#if TARGET_64BIT
+ private static ReadOnlySpan Pow10BigNumTableIndices =>
[
- 1, // 10^0
- 10, // 10^1
- 100, // 10^2
- 1000, // 10^3
- 10000, // 10^4
- 100000, // 10^5
- 1000000, // 10^6
- 10000000, // 10^7
- // These last two are accessed only by MultiplyPow10.
- 100000000, // 10^8
- 1000000000 // 10^9
+ 0, // 10^16
+ 2, // 10^32
+ 5, // 10^64
+ 10, // 10^128
+ 18, // 10^256
+ 33, // 10^512
+ 61, // 10^1024
];
+ private static ReadOnlySpan Pow10BigNumTable =>
+ [
+ // 10^16
+ 1, // _length
+ 0x002386F26FC10000, // _blocks
+
+ // 10^32
+ 2, // _length
+ 0x85ACEF8100000000, // _blocks
+ 0x000004EE2D6D415B,
+
+ // 10^64
+ 4, // _length
+ 0x0000000000000000, // _blocks
+ 0x6E38ED64BF6A1F01,
+ 0xE93FF9F4DAA797ED,
+ 0x0000000000184F03,
+
+ // 10^128
+ 7, // _length
+ 0x0000000000000000, // _blocks
+ 0x0000000000000000,
+ 0x03DF99092E953E01,
+ 0x2374E42F0F1538FD,
+ 0xC404DC08D3CFF5EC,
+ 0xA6337F19BCCDB0DA,
+ 0x0000024EE91F2603,
+
+ // 10^256
+ 14, // _length
+ 0x0000000000000000, // _blocks
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0xBED3875B982E7C01,
+ 0x12152F87D8D99F72,
+ 0xCF4A6E706BDE50C6,
+ 0x26B2716ED595D80F,
+ 0x1D153624ADC666B0,
+ 0x63FF540E3C42D35A,
+ 0x65F9EF17CC5573C0,
+ 0x80DCC7F755BC28F2,
+ 0x5FDCEFCEF46EEDDC,
+ 0x00000000000553F7,
+
+ // 10^512
+ 27, // _length
+ 0x0000000000000000, // _blocks
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x77F27267FC6CF801,
+ 0x5D96976F8F9546DC,
+ 0xC31E1AD9B83A8A97,
+ 0x94E6574746C40513,
+ 0x4475B579C88976C1,
+ 0xAA1DA1BF28F8733B,
+ 0x1E25CFEA703ED321,
+ 0xBC51FB2EB21A2F22,
+ 0xBFA3EDAC96E14F5D,
+ 0xE7FC7153329C57AE,
+ 0x85A91924C3FC0695,
+ 0xB2908EE0F95F635E,
+ 0x1366732A93ABADE4,
+ 0x69BE5B0E9449775C,
+ 0xB099BC817343AFAC,
+ 0xA269974845A71D46,
+ 0x8A0B1F138CB07303,
+ 0xC1D238D98CAB8A97,
+ 0x0000001C633415D4,
+
+ // 10^1024
+ 54, // _length
+ 0x0000000000000000, // _blocks
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0xF55B2B722919F001,
+ 0x1EC29F866E7C215B,
+ 0x15C51A88991C4E87,
+ 0x4C7D1E1A140AC535,
+ 0x0ED1440ECC2CD819,
+ 0x7DE16CFB896634EE,
+ 0x9FCE837D1E43F61F,
+ 0x233E55C7231D2B9C,
+ 0xF451218B65DC60D7,
+ 0xC96359861C5CD134,
+ 0xA7E89431922BBB9F,
+ 0x62BE695A9F9F2A07,
+ 0x045B7A748E1042C4,
+ 0x8AD822A51ABE1DE3,
+ 0xD814B505BA34C411,
+ 0x8FC51A16BF3FDEB3,
+ 0xF56DEEECB1B896BC,
+ 0xB6F4654B31FB6BFD,
+ 0x6B7595FB101A3616,
+ 0x80D98089DC1A47FE,
+ 0x9A20288280BDA5A5,
+ 0xFC8F1F9031EB0F66,
+ 0xE26A7B7E976A3310,
+ 0x3CE3A0B8DF68368A,
+ 0x75A351A28E4262CE,
+ 0x445975836CB0B6C9,
+ 0xC356E38A31B5653F,
+ 0x0190FBA035FAABA6,
+ 0x88BC491B9FC4ED52,
+ 0x005B80411640114A,
+ 0x1E8D4649F4F3235E,
+ 0x73C5534936A8DE06,
+ 0xC1A6970CA7E6BD2A,
+ 0xD2DB49EF47187094,
+ 0xAE6209D4926C3F5B,
+ 0x34F4A3C62D433949,
+ 0xD9D61A05D4305D94,
+ 0x0000000000000325,
+
+ // 5 trailing zero blocks so the last entry can be reinterpreted as a full BigInteger
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ 0x0000000000000000,
+ ];
+#else
private static ReadOnlySpan Pow10BigNumTableIndices =>
[
0, // 10^8
@@ -316,114 +482,53 @@ internal ref struct BigInteger
0x00000000,
0x00000000,
];
+#endif
- private int _length;
+ private nint _length;
private BlocksBuffer _blocks;
- public static void Add(scoped ref BigInteger lhs, scoped ref BigInteger rhs, out BigInteger result)
+ public static void Add(scoped ref readonly BigInteger lhs, scoped ref readonly BigInteger rhs, out BigInteger result)
{
Unsafe.SkipInit(out result);
- // determine which operand has the smaller length
- ref BigInteger large = ref (lhs._length < rhs._length) ? ref rhs : ref lhs;
- ref BigInteger small = ref (lhs._length < rhs._length) ? ref lhs : ref rhs;
-
- int largeLength = large._length;
- int smallLength = small._length;
-
- // The output will be at least as long as the largest input
- result._length = largeLength;
+ // Order operands so the longer one drives the loop.
+ ref readonly BigInteger large = ref (lhs._length < rhs._length ? ref rhs : ref lhs);
+ ref readonly BigInteger small = ref (lhs._length < rhs._length ? ref lhs : ref rhs);
- // Add each block and add carry the overflow to the next block
- ulong carry = 0;
+ int largeLength = (int)large._length;
+ int smallLength = (int)small._length;
- int largeIndex = 0;
- int smallIndex = 0;
- int resultIndex = 0;
-
- while (smallIndex < smallLength)
- {
- ulong sum = carry + large._blocks[largeIndex] + small._blocks[smallIndex];
- carry = sum >> 32;
- result._blocks[resultIndex] = (uint)sum;
-
- largeIndex++;
- smallIndex++;
- resultIndex++;
- }
-
- // Add the carry to any blocks that only exist in the large operand
- while (largeIndex < largeLength)
+ if (smallLength == 0)
{
- ulong sum = carry + large._blocks[largeIndex];
- carry = sum >> 32;
- result._blocks[resultIndex] = (uint)sum;
-
- largeIndex++;
- resultIndex++;
+ SetValue(out result, in large);
+ return;
}
- int resultLength = largeLength;
+ Debug.Assert(unchecked((uint)largeLength) < MaxBlockCount);
- // If there's still a carry, append a new block
- if (carry != 0)
+ if (unchecked((uint)largeLength) >= MaxBlockCount)
{
- Debug.Assert(carry == 1);
- Debug.Assert(resultIndex == resultLength);
- Debug.Assert(unchecked((uint)resultLength) < MaxBlockCount);
+ // We shouldn't reach here, and the above assert will help flag this
+ // during testing, but we'll ensure that we return a safe value of
+ // zero in the case we end up overflowing in any way.
- if (unchecked((uint)resultLength) >= MaxBlockCount)
- {
- // We shouldn't reach here, and the above assert will help flag this
- // during testing, but we'll ensure that we return a safe value of
- // zero in the case we end up overflowing in any way.
+ SetZero(out result);
+ return;
+ }
- SetZero(out result);
- return;
- }
+ // The shared kernel writes largeLength + 1 blocks; the extra holds the carry-out.
+ Span resultBlocks = result._blocks;
+ BigIntegerCalculator.Add(large._blocks[..largeLength], small._blocks[..smallLength], resultBlocks[..(largeLength + 1)]);
- result._blocks[resultIndex] = 1;
- result._length++;
- }
+ result._length = (resultBlocks[largeLength] != 0) ? (largeLength + 1) : largeLength;
}
- public static int Compare(scoped ref BigInteger lhs, scoped ref BigInteger rhs)
+ public static int Compare(scoped ref readonly BigInteger lhs, scoped ref readonly BigInteger rhs)
{
Debug.Assert(unchecked((uint)lhs._length) <= MaxBlockCount);
Debug.Assert(unchecked((uint)rhs._length) <= MaxBlockCount);
- int lhsLength = lhs._length;
- int rhsLength = rhs._length;
-
- int lengthDelta = lhsLength - rhsLength;
-
- if (lengthDelta != 0)
- {
- return lengthDelta;
- }
-
- if (lhsLength == 0)
- {
- Debug.Assert(rhsLength == 0);
- return 0;
- }
-
- for (int index = lhsLength - 1; index >= 0; index--)
- {
- long delta = (long)lhs._blocks[index] - rhs._blocks[index];
-
- if (delta != 0)
- {
- return delta > 0 ? 1 : -1;
- }
- }
-
- return 0;
- }
-
- public static int CountSignificantBits(uint value)
- {
- return 32 - BitOperations.LeadingZeroCount(value);
+ return BigIntegerCalculator.Compare(lhs._blocks[..(int)lhs._length], rhs._blocks[..(int)rhs._length]);
}
public static int CountSignificantBits(ulong value)
@@ -431,7 +536,7 @@ public static int CountSignificantBits(ulong value)
return 64 - BitOperations.LeadingZeroCount(value);
}
- public static int CountSignificantBits(ref BigInteger value)
+ public static int CountSignificantBits(scoped ref readonly BigInteger value)
{
if (value.IsZero())
{
@@ -441,11 +546,11 @@ public static int CountSignificantBits(ref BigInteger value)
// We don't track any unused blocks, so we only need to do a BSR on the
// last index and add that to the number of bits we skipped.
- int lastIndex = value._length - 1;
- return (lastIndex * BitsPerBlock) + CountSignificantBits(value._blocks[lastIndex]);
+ int lastIndex = (int)value._length - 1;
+ return (lastIndex * BitsPerBlock) + (BitsPerBlock - BitOperations.LeadingZeroCount(value._blocks[lastIndex]));
}
- public static void DivRem(scoped ref BigInteger lhs, scoped ref BigInteger rhs, out BigInteger quo, out BigInteger rem)
+ public static void DivRem(scoped ref readonly BigInteger lhs, scoped ref readonly BigInteger rhs, out BigInteger quo, out BigInteger rem)
{
Unsafe.SkipInit(out quo);
@@ -461,14 +566,14 @@ public static void DivRem(scoped ref BigInteger lhs, scoped ref BigInteger rhs,
return;
}
- int lhsLength = lhs._length;
- int rhsLength = rhs._length;
+ int lhsLength = (int)lhs._length;
+ int rhsLength = (int)rhs._length;
if ((lhsLength == 1) && (rhsLength == 1))
{
- (uint quotient, uint remainder) = Math.DivRem(lhs._blocks[0], rhs._blocks[0]);
- SetUInt32(out quo, quotient);
- SetUInt32(out rem, remainder);
+ (nuint quotient, nuint remainder) = Math.DivRem(lhs._blocks[0], rhs._blocks[0]);
+ SetBlock(out quo, quotient);
+ SetBlock(out rem, remainder);
return;
}
@@ -478,14 +583,13 @@ public static void DivRem(scoped ref BigInteger lhs, scoped ref BigInteger rhs,
int quoLength = lhsLength;
- ulong rhsValue = rhs._blocks[0];
- ulong carry = 0;
+ nuint rhsValue = rhs._blocks[0];
+ nuint carry = 0;
for (int i = quoLength - 1; i >= 0; i--)
{
- ulong value = (carry << 32) | lhs._blocks[i];
- ulong digit;
- (digit, carry) = Math.DivRem(value, rhsValue);
+ // carry is the running remainder, always < rhsValue, so the quotient fits a block.
+ nuint digit = BigIntegerCalculator.DivRem(carry, lhs._blocks[i], rhsValue, out carry);
if ((digit == 0) && (i == (quoLength - 1)))
{
@@ -493,12 +597,12 @@ public static void DivRem(scoped ref BigInteger lhs, scoped ref BigInteger rhs,
}
else
{
- quo._blocks[i] = (uint)digit;
+ quo._blocks[i] = digit;
}
}
quo._length = quoLength;
- SetUInt32(out rem, (uint)carry);
+ SetBlock(out rem, carry);
return;
}
@@ -506,142 +610,27 @@ public static void DivRem(scoped ref BigInteger lhs, scoped ref BigInteger rhs,
{
// Handle the case where we have no quotient
SetZero(out quo);
- SetValue(out rem, ref lhs);
+ SetValue(out rem, in lhs);
return;
}
else
{
int quoLength = lhsLength - rhsLength + 1;
- SetValue(out rem, ref lhs);
- int remLength = lhsLength;
+ SetValue(out rem, in lhs);
- // Executes the "grammar-school" algorithm for computing q = a / b.
- // Before calculating q_i, we get more bits into the highest bit
- // block of the divisor. Thus, guessing digits of the quotient
- // will be more precise. Additionally we'll get r = a % b.
-
- uint divHi = rhs._blocks[rhsLength - 1];
- uint divLo = rhs._blocks[rhsLength - 2];
-
- // We measure the leading zeros of the divisor
- int shiftLeft = BitOperations.LeadingZeroCount(divHi);
- int shiftRight = 32 - shiftLeft;
-
- // And, we make sure the most significant bit is set
- if (shiftLeft > 0)
- {
- divHi = (divHi << shiftLeft) | (divLo >> shiftRight);
- divLo <<= shiftLeft;
-
- if (rhsLength > 2)
- {
- divLo |= rhs._blocks[rhsLength - 3] >> shiftRight;
- }
- }
+ // The shared kernel divides in place: rem holds the dividend on entry and the
+ // remainder (in its low rhsLength blocks) on exit, with quo receiving the quotient.
+ Span quoBlocks = quo._blocks;
+ BigIntegerCalculator.DivideGrammarSchool(rem._blocks[..lhsLength], rhs._blocks[..rhsLength], quoBlocks[..quoLength]);
- // Then, we divide all of the bits as we would do it using
- // pen and paper: guessing the next digit, subtracting, ...
- for (int i = lhsLength; i >= rhsLength; i--)
- {
- int n = i - rhsLength;
- uint t = i < lhsLength ? rem._blocks[i] : 0;
-
- ulong valHi = ((ulong)t << 32) | rem._blocks[i - 1];
- uint valLo = i > 1 ? rem._blocks[i - 2] : 0;
-
- // We shifted the divisor, we shift the dividend too
- if (shiftLeft > 0)
- {
- valHi = (valHi << shiftLeft) | (valLo >> shiftRight);
- valLo <<= shiftLeft;
-
- if (i > 2)
- {
- valLo |= rem._blocks[i - 3] >> shiftRight;
- }
- }
-
- // First guess for the current digit of the quotient,
- // which naturally must have only 32 bits...
- ulong digit = valHi / divHi;
-
- if (digit > uint.MaxValue)
- {
- digit = uint.MaxValue;
- }
-
- // Our first guess may be a little bit to big
- while (DivideGuessTooBig(digit, valHi, valLo, divHi, divLo))
- {
- digit--;
- }
-
- if (digit > 0)
- {
- // rem and rhs have different lifetimes here and compiler is warning
- // about potential for one to copy into the other. This is a place
- // ref scoped parameters would alleviate.
- // https://github.com/dotnet/roslyn/issues/64393
-#pragma warning disable CS9080
- // Now it's time to subtract our current quotient
- uint carry = SubtractDivisor(ref rem, n, ref rhs, digit);
-
- if (carry != t)
- {
- Debug.Assert(carry == t + 1);
-
- // Our guess was still exactly one too high
- carry = AddDivisor(ref rem, n, ref rhs);
- digit--;
-
- Debug.Assert(carry == 1);
- }
-#pragma warning restore CS9080
- }
-
- // We have the digit!
- if (quoLength != 0)
- {
- if ((digit == 0) && (n == (quoLength - 1)))
- {
- quoLength--;
- }
- else
- {
- quo._blocks[n] = (uint)digit;
- }
- }
-
- if (i < remLength)
- {
- remLength--;
- }
- }
-
- quo._length = quoLength;
-
- // We need to check for the case where remainder is zero
-
- for (int i = remLength - 1; i >= 0; i--)
- {
- if (rem._blocks[i] == 0)
- {
- remLength--;
- }
- else
- {
- // As soon as we find a non-zero block, the rest of remainder is significant
- break;
- }
- }
-
- rem._length = remLength;
+ quo._length = BigIntegerCalculator.ActualLength(quoBlocks[..quoLength]);
+ rem._length = BigIntegerCalculator.ActualLength(rem._blocks[..rhsLength]);
}
}
- public static uint HeuristicDivide(ref BigInteger dividend, ref BigInteger divisor)
+ public static uint HeuristicDivide(scoped ref BigInteger dividend, scoped ref readonly BigInteger divisor)
{
- int divisorLength = divisor._length;
+ int divisorLength = (int)divisor._length;
if (dividend._length < divisorLength)
{
@@ -652,83 +641,74 @@ public static uint HeuristicDivide(ref BigInteger dividend, ref BigInteger divis
// Reference inequality:
// a/b - floor(floor(a)/(floor(b) + 1)) < 2
int lastIndex = divisorLength - 1;
- uint quotient = dividend._blocks[lastIndex] / (divisor._blocks[lastIndex] + 1);
+ nuint quotient = dividend._blocks[lastIndex] / (divisor._blocks[lastIndex] + 1);
if (quotient != 0)
{
- // Now we use our estimated quotient to update each block of dividend.
- // dividend = dividend - divisor * quotient
- int index = 0;
-
- ulong borrow = 0;
- ulong carry = 0;
-
- do
+ // dividend -= divisor * quotient. This runs once per output digit in Dragon4's hot
+ // loop, so the widening multiply-subtract is kept fully inline: routing through the
+ // shared span kernels (SubMul1) regressed ToString("G50") ~40%+ (out-of-line call +
+ // Span setup spilling the block pointers), and even an inlinable scalar helper left
+ // ~8% on the table because the UInt128 body did not inline. The `nint.Size` branch
+ // constant-folds so only the native-width loop is emitted. The estimate never
+ // overshoots enough to underflow, so the trailing borrow is discarded.
+ nuint borrow = 0;
+ if (nint.Size == 8)
{
- ulong product = ((ulong)divisor._blocks[index] * quotient) + carry;
- carry = product >> 32;
-
- ulong difference = (ulong)dividend._blocks[index] - (uint)product - borrow;
- borrow = (difference >> 32) & 1;
-
- dividend._blocks[index] = (uint)difference;
-
- index++;
+ for (int i = 0; i < divisorLength; i++)
+ {
+ UInt128 product = (UInt128)(ulong)divisor._blocks[i] * quotient + borrow;
+ nuint low = (nuint)product.Lower;
+ nuint orig = dividend._blocks[i];
+ dividend._blocks[i] = orig - low;
+ borrow = (nuint)product.Upper + ((orig < low) ? (nuint)1 : 0);
+ }
}
- while (index < divisorLength);
-
- // Remove all leading zero blocks from dividend
- while ((divisorLength > 0) && (dividend._blocks[divisorLength - 1] == 0))
+ else
{
- divisorLength--;
+ for (int i = 0; i < divisorLength; i++)
+ {
+ ulong product = (ulong)(uint)divisor._blocks[i] * quotient + borrow;
+ nuint low = (nuint)(uint)product;
+ nuint orig = dividend._blocks[i];
+ dividend._blocks[i] = orig - low;
+ borrow = (nuint)(uint)(product >> 32) + ((orig < low) ? (nuint)1 : 0);
+ }
}
+ // Remove all leading zero blocks from dividend
+ divisorLength = BigIntegerCalculator.ActualLength(dividend._blocks[..divisorLength]);
dividend._length = divisorLength;
}
// If the dividend is still larger than the divisor, we overshot our estimate quotient. To correct,
// we increment the quotient and subtract one more divisor from the dividend (Because we guaranteed the error range).
- if (Compare(ref dividend, ref divisor) >= 0)
+ if (Compare(ref dividend, in divisor) >= 0)
{
quotient++;
- // dividend = dividend - divisor
- int index = 0;
- ulong borrow = 0;
-
- do
+ // dividend -= divisor. This is the cold correction path (only on overshoot), so it
+ // reuses the shared SubWithBorrow leaf, which inlines into this local loop.
+ nuint borrow = 0;
+ for (int i = 0; i < divisorLength; i++)
{
- ulong difference = (ulong)dividend._blocks[index] - divisor._blocks[index] - borrow;
- borrow = (difference >> 32) & 1;
-
- dividend._blocks[index] = (uint)difference;
-
- index++;
+ dividend._blocks[i] = BigIntegerCalculator.SubWithBorrow(dividend._blocks[i], divisor._blocks[i], borrow, out nuint nextBorrow);
+ borrow = nextBorrow;
}
- while (index < divisorLength);
// Remove all leading zero blocks from dividend
- while ((divisorLength > 0) && (dividend._blocks[divisorLength - 1] == 0))
- {
- divisorLength--;
- }
-
+ divisorLength = BigIntegerCalculator.ActualLength(dividend._blocks[..divisorLength]);
dividend._length = divisorLength;
}
- return quotient;
+ Debug.Assert(quotient < 10);
+ return (uint)quotient;
}
- public static void Multiply(scoped ref BigInteger lhs, uint value, out BigInteger result)
+ public static void Multiply(scoped ref readonly BigInteger lhs, nuint value, out BigInteger result)
{
Unsafe.SkipInit(out result);
- if (lhs._length <= 1)
- {
- SetUInt64(out result, (ulong)lhs.ToUInt32() * value);
- return;
- }
-
if (value <= 1)
{
if (value == 0)
@@ -737,23 +717,16 @@ public static void Multiply(scoped ref BigInteger lhs, uint value, out BigIntege
}
else
{
- SetValue(out result, ref lhs);
+ SetValue(out result, in lhs);
}
return;
}
- int lhsLength = lhs._length;
- int index = 0;
- uint carry = 0;
+ int lhsLength = (int)lhs._length;
- while (index < lhsLength)
- {
- ulong product = ((ulong)lhs._blocks[index] * value) + carry;
- result._blocks[index] = (uint)product;
- carry = (uint)(product >> 32);
-
- index++;
- }
+ // Mul1 is safe in place, so this also serves the result-aliases-lhs callers.
+ Span resultBlocks = result._blocks;
+ nuint carry = BigIntegerCalculator.Mul1(resultBlocks[..lhsLength], lhs._blocks[..lhsLength], value);
int resultLength = lhsLength;
@@ -771,42 +744,42 @@ public static void Multiply(scoped ref BigInteger lhs, uint value, out BigIntege
return;
}
- result._blocks[index] = carry;
+ resultBlocks[resultLength] = carry;
resultLength += 1;
}
result._length = resultLength;
}
- public static void Multiply(scoped ref BigInteger lhs, scoped ref BigInteger rhs, out BigInteger result)
+ public static void Multiply(scoped ref readonly BigInteger lhs, scoped ref readonly BigInteger rhs, out BigInteger result)
{
Unsafe.SkipInit(out result);
if (lhs._length <= 1)
{
- Multiply(ref rhs, lhs.ToUInt32(), out result);
+ Multiply(in rhs, lhs.ToBlock(), out result);
return;
}
if (rhs._length <= 1)
{
- Multiply(ref lhs, rhs.ToUInt32(), out result);
+ Multiply(in lhs, rhs.ToBlock(), out result);
return;
}
ref readonly BigInteger large = ref lhs;
- int largeLength = lhs._length;
+ int largeLength = (int)lhs._length;
ref readonly BigInteger small = ref rhs;
- int smallLength = rhs._length;
+ int smallLength = (int)rhs._length;
if (largeLength < smallLength)
{
large = ref rhs;
- largeLength = rhs._length;
+ largeLength = (int)rhs._length;
small = ref lhs;
- smallLength = lhs._length;
+ smallLength = (int)lhs._length;
}
int maxResultLength = smallLength + largeLength;
@@ -822,42 +795,17 @@ public static void Multiply(scoped ref BigInteger lhs, scoped ref BigInteger rhs
return;
}
- // Zero out result internal blocks.
+ // The shared kernel accumulates into a zero-initialized result buffer.
result._length = maxResultLength;
- result.Clear(maxResultLength);
-
- int smallIndex = 0;
- int resultStartIndex = 0;
-
- while (smallIndex < smallLength)
- {
- // Multiply each block of large BigNum.
- if (small._blocks[smallIndex] != 0)
- {
- int largeIndex = 0;
- int resultIndex = resultStartIndex;
-
- ulong carry = 0;
- do
- {
- ulong product = result._blocks[resultIndex] + ((ulong)small._blocks[smallIndex] * large._blocks[largeIndex]) + carry;
- carry = product >> 32;
- result._blocks[resultIndex] = (uint)product;
+ ReadOnlySpan largeBlocks = ((ReadOnlySpan)large._blocks).Slice(0, largeLength);
+ ReadOnlySpan smallBlocks = ((ReadOnlySpan)small._blocks).Slice(0, smallLength);
+ Span resultBlocks = ((Span)result._blocks).Slice(0, maxResultLength);
- resultIndex++;
- largeIndex++;
- }
- while (largeIndex < largeLength);
+ resultBlocks.Clear();
+ BigIntegerCalculator.MultiplyNaive(largeBlocks, smallBlocks, resultBlocks);
- result._blocks[resultIndex] = (uint)carry;
- }
-
- smallIndex++;
- resultStartIndex++;
- }
-
- if ((maxResultLength > 0) && (result._blocks[maxResultLength - 1] == 0))
+ if ((maxResultLength > 0) && (resultBlocks[maxResultLength - 1] == 0))
{
result._length--;
}
@@ -867,7 +815,7 @@ public static void Pow2(int exponent, out BigInteger result)
{
Unsafe.SkipInit(out result);
- int blocksToShift = DivRem32(exponent, out int remainingBitsToShift);
+ int blocksToShift = DivRemBlock(exponent, out int remainingBitsToShift);
result._length = blocksToShift + 1;
Debug.Assert(unchecked((uint)result._length) <= MaxBlockCount);
@@ -886,48 +834,46 @@ public static void Pow2(int exponent, out BigInteger result)
{
result.Clear(blocksToShift);
}
- result._blocks[blocksToShift] = 1U << remainingBitsToShift;
+ result._blocks[blocksToShift] = (nuint)1 << remainingBitsToShift;
+ }
+
+ // Reinterprets the power stored at the given block index of Pow10BigNumTable as a BigInteger.
+ // The table packs each power as `nint _length` followed by its blocks, all in native-word
+ // sized elements, so a block index maps to a byte offset of index * sizeof(nuint). Viewing the
+ // bytes keeps the per-bitness element type (ulong on 64-bit, uint on 32-bit) out of the cast.
+ private static ref readonly BigInteger Pow10BigNum(int index)
+ {
+ ReadOnlySpan table = MemoryMarshal.AsBytes(Pow10BigNumTable);
+ int offset = index * sizeof(nuint);
+ Debug.Assert((offset + sizeof(BigInteger)) <= table.Length);
+ return ref Unsafe.As(ref MemoryMarshal.GetReference(table.Slice(offset)));
}
public static void Pow10(uint exponent, out BigInteger result)
{
- // We leverage two arrays - Pow10UInt32Table and Pow10BigNumTable to speed up the Pow10 calculation.
+ // We leverage two arrays - Pow10Table and Pow10BigNumTable to speed up the Pow10 calculation.
//
- // Pow10UInt32Table stores the results of 10^0 to 10^7.
- // Pow10BigNumTable stores the results of 10^8, 10^16, 10^32, 10^64, 10^128, 10^256, and 10^512
+ // The low Pow10SmallExpBits of the exponent select 10^smallExp, which fits in a single block
+ // and is applied with one lookup+multiply. Each remaining bit selects a power from
+ // Pow10BigNumTable, whose entries are 10^(2^Pow10SmallExpBits), 10^(2^(Pow10SmallExpBits+1)), ...
+ // So the result is 10^smallExp * product(selected big powers).
//
- // For example, let's say exp = 0b111111. We can split the exp to two parts, one is small exp,
- // which 10^smallExp can be represented as uint, another part is 10^bigExp, which must be represented as BigNum.
- // So the result should be 10^smallExp * 10^bigExp.
+ // The split is a clean power-of-two boundary: bits below Pow10SmallExpBits are small, the rest
+ // are big. That requires 10^(2^Pow10SmallExpBits - 1) to still fit in a single block, which is
+ // why the boundary tracks the block width (see Pow10SmallExpBits).
//
- // Calculating 10^smallExp is simple, we just lookup the 10^smallExp from Pow10UInt32Table.
- // But here's a bad news: although uint can represent 10^9, exp 9's binary representation is 1001.
- // That means 10^(1011), 10^(1101), 10^(1111) all cannot be stored as uint, we cannot easily say something like:
- // "Any bits <= 3 is small exp, any bits > 3 is big exp". So instead of involving 10^8, 10^9 to Pow10UInt32Table,
- // consider 10^8 and 10^9 as a bigNum, so they fall into Pow10BigNumTable. Now we can have a simple rule:
- // "Any bits <= 3 is small exp, any bits > 3 is big exp".
- //
- // For 0b111111, we first calculate 10^(smallExp), which is 10^(7), now we can shift right 3 bits, prepare to calculate the bigExp part,
- // the exp now becomes 0b000111.
- //
- // Apparently the lowest bit of bigExp should represent 10^8 because we have already shifted 3 bits for smallExp, so Pow10BigNumTable[0] = 10^8.
- // Now let's shift exp right 1 bit, the lowest bit should represent 10^(8 * 2) = 10^16, and so on...
- //
- // That's why we just need the values of Pow10BigNumTable be power of 2.
+ // For example on 32-bit (Pow10SmallExpBits == 3) exp = 0b111111 splits into 10^7 (0b111) and the
+ // remaining 0b111 selecting 10^8 * 10^16 * 10^32.
//
// More details of this implementation can be found at: https://github.com/dotnet/coreclr/pull/12894#discussion_r128890596
- // Validate that `Pow10BigNumTable` has exactly enough trailing elements to fill a BigInteger (which contains MaxBlockCount + 1 elements)
- // We validate here, since this is the only current consumer of the array
- Debug.Assert((Pow10BigNumTableIndices[^1] + MaxBlockCount + 2) == Pow10BigNumTable.Length);
-
- SetUInt32(out BigInteger temp1, Pow10UInt32Table[(int)(exponent & 0x7)]);
+ SetBlock(out BigInteger temp1, (nuint)Pow10Table[(int)(exponent & Pow10SmallExpMask)]);
ref BigInteger lhs = ref temp1;
SetZero(out BigInteger temp2);
ref BigInteger product = ref temp2;
- exponent >>= 3;
+ exponent >>= Pow10SmallExpBits;
int index = 0;
while (exponent != 0)
@@ -935,16 +881,9 @@ public static void Pow10(uint exponent, out BigInteger result)
// If the current bit is set, multiply it with the corresponding power of 10
if ((exponent & 1) != 0)
{
- // Multiply into the next temporary
- // It is safe to reinterpret the memory at the indexed position of Pow10BigNumTable as a BigInteger,
- // because Pow10BigNumTable is laid out such that each indexed position contains a valid BigInteger
- // representation with the correct structure and alignment.
-
- int pow10BigNumTableIndex = Pow10BigNumTableIndices[index];
- Debug.Assert((pow10BigNumTableIndex + ((sizeof(BigInteger) + sizeof(uint) - 1) / sizeof(uint))) < Pow10BigNumTable.Length);
-
- ref BigInteger rhs = ref Unsafe.As(ref Unsafe.AsRef(in Pow10BigNumTable[pow10BigNumTableIndex]));
- Multiply(ref lhs, ref rhs, out product);
+ // Multiply into the next temporary.
+ ref readonly BigInteger rhs = ref Pow10BigNum(Pow10BigNumTableIndices[index]);
+ Multiply(ref lhs, in rhs, out product);
// Swap to the next temporary
ref BigInteger temp = ref product;
@@ -960,101 +899,12 @@ public static void Pow10(uint exponent, out BigInteger result)
SetValue(out result, ref lhs);
}
- private static uint AddDivisor(ref BigInteger lhs, int lhsStartIndex, ref BigInteger rhs)
+ public void Add(nuint value)
{
- int lhsLength = lhs._length;
- int rhsLength = rhs._length;
-
- Debug.Assert(lhsLength >= 0);
- Debug.Assert(rhsLength >= 0);
- Debug.Assert(lhsLength >= rhsLength);
-
- // Repairs the dividend, if the last subtract was too much
-
- ulong carry = 0UL;
-
- for (int i = 0; i < rhsLength; i++)
- {
- ref uint lhsValue = ref lhs._blocks[lhsStartIndex + i];
-
- ulong digit = lhsValue + carry + rhs._blocks[i];
- lhsValue = unchecked((uint)digit);
- carry = digit >> 32;
- }
-
- return (uint)carry;
- }
-
- private static bool DivideGuessTooBig(ulong q, ulong valHi, uint valLo, uint divHi, uint divLo)
- {
- Debug.Assert(q <= 0xFFFFFFFF);
-
- // We multiply the two most significant limbs of the divisor
- // with the current guess for the quotient. If those are bigger
- // than the three most significant limbs of the current dividend
- // we return true, which means the current guess is still too big.
-
- ulong chkHi = divHi * q;
- ulong chkLo = divLo * q;
-
- chkHi += chkLo >> 32;
- chkLo &= uint.MaxValue;
-
- if (chkHi < valHi)
- return false;
-
- if (chkHi > valHi)
- return true;
-
- if (chkLo < valLo)
- return false;
-
- if (chkLo > valLo)
- return true;
-
- return false;
- }
-
- private static uint SubtractDivisor(ref BigInteger lhs, int lhsStartIndex, ref BigInteger rhs, ulong q)
- {
- int lhsLength = lhs._length - lhsStartIndex;
- int rhsLength = rhs._length;
-
- Debug.Assert(lhsLength >= 0);
- Debug.Assert(rhsLength >= 0);
- Debug.Assert(lhsLength >= rhsLength);
- Debug.Assert(q <= uint.MaxValue);
-
- // Combines a subtract and a multiply operation, which is naturally
- // more efficient than multiplying and then subtracting...
-
- ulong carry = 0;
-
- for (int i = 0; i < rhsLength; i++)
- {
- carry += rhs._blocks[i] * q;
- uint digit = unchecked((uint)carry);
- carry >>= 32;
-
- ref uint lhsValue = ref lhs._blocks[lhsStartIndex + i];
-
- if (lhsValue < digit)
- {
- carry++;
- }
-
- lhsValue = unchecked(lhsValue - digit);
- }
-
- return (uint)carry;
- }
-
- public void Add(uint value)
- {
- int length = _length;
+ int length = (int)_length;
if (length == 0)
{
- SetUInt32(out this, value);
+ SetBlock(out this, value);
return;
}
@@ -1091,15 +941,72 @@ public void Add(uint value)
_length = length + 1;
}
- public readonly uint GetBlock(int index)
+ public readonly nuint GetBlock(int index)
{
Debug.Assert((uint)index < _length);
return _blocks[index];
}
+ public readonly ulong GetBits64(int bitIndex)
+ {
+ // Extracts the 64-bit window starting at bitIndex, i.e. bits [bitIndex, bitIndex + 64),
+ // funnel-shifting across native blocks. Callers guarantee the window is in range.
+ int block = bitIndex >> BlockShift;
+ int shift = bitIndex & (BitsPerBlock - 1);
+
+#if TARGET_64BIT
+ Debug.Assert((uint)block < (uint)_length);
+ ulong low = (ulong)_blocks[block];
+
+ if (shift == 0)
+ {
+ return low;
+ }
+
+ Debug.Assert((uint)(block + 1) < (uint)_length);
+ ulong high = (ulong)_blocks[block + 1];
+ return (low >> shift) | (high << (BitsPerBlock - shift));
+#else
+ Debug.Assert((uint)(block + 1) < (uint)_length);
+ ulong low = _blocks[block] | ((ulong)_blocks[block + 1] << BitsPerBlock);
+
+ if (shift == 0)
+ {
+ return low;
+ }
+
+ Debug.Assert((uint)(block + 2) < (uint)_length);
+ ulong high = _blocks[block + 2];
+ return (low >> shift) | (high << (64 - shift));
+#endif
+ }
+
+ public readonly bool HasZeroTail(int bitCount)
+ {
+ // Returns true when the lowest bitCount bits are all zero.
+ Debug.Assert((uint)bitCount <= (uint)(_length * BitsPerBlock));
+
+ int fullBlocks = bitCount >> BlockShift;
+
+ if (((ReadOnlySpan)_blocks)[..fullBlocks].ContainsAnyExcept((nuint)0))
+ {
+ return false;
+ }
+
+ int remainingBits = bitCount & (BitsPerBlock - 1);
+
+ if (remainingBits != 0)
+ {
+ nuint mask = ((nuint)1 << remainingBits) - 1;
+ return (_blocks[fullBlocks] & mask) == 0;
+ }
+
+ return true;
+ }
+
public readonly int GetLength()
{
- return _length;
+ return (int)_length;
}
public readonly bool IsZero()
@@ -1107,7 +1014,7 @@ public readonly bool IsZero()
return _length == 0;
}
- public void Multiply(uint value)
+ public void Multiply(nuint value)
{
Multiply(ref this, value, out this);
}
@@ -1116,7 +1023,7 @@ public void Multiply(scoped ref BigInteger value)
{
if (value._length <= 1)
{
- Multiply(ref this, value.ToUInt32(), out this);
+ Multiply(ref this, value.ToBlock(), out this);
}
else
{
@@ -1132,19 +1039,33 @@ public void Multiply10()
return;
}
- int index = 0;
- int length = _length;
- ulong carry = 0;
+ int length = (int)_length;
- do
- {
- ulong block = _blocks[index];
- ulong product = (block << 3) + (block << 1) + carry;
- carry = product >> 32;
- _blocks[index] = (uint)product;
+ // Multiply-by-10 is called once per output digit in Dragon4's hot loop. Routing it
+ // through the shared Mul1 kernel costs a non-inlined call plus span setup per digit
+ // (the unrolled UInt128 body is too large to inline), which measurably regresses
+ // ToString of long fixed-precision formats. Keeping the single-block multiply inline
+ // here avoids that; nint.Size constant-folds to the native-width loop.
+ nuint carry = 0;
- index++;
- } while (index < length);
+ if (nint.Size == 8)
+ {
+ for (int i = 0; i < length; i++)
+ {
+ UInt128 product = (UInt128)(ulong)_blocks[i] * 10 + carry;
+ _blocks[i] = (nuint)product.Lower;
+ carry = (nuint)product.Upper;
+ }
+ }
+ else
+ {
+ for (int i = 0; i < length; i++)
+ {
+ ulong product = (ulong)(uint)_blocks[i] * 10 + carry;
+ _blocks[i] = (nuint)(uint)product;
+ carry = (nuint)(uint)(product >> 32);
+ }
+ }
if (carry != 0)
{
@@ -1160,16 +1081,16 @@ public void Multiply10()
return;
}
- _blocks[index] = (uint)carry;
+ _blocks[length] = carry;
_length = length + 1;
}
}
public void MultiplyPow10(uint exponent)
{
- if (exponent <= 9)
+ if (exponent < (uint)Pow10Table.Length)
{
- Multiply(Pow10UInt32Table[(int)exponent]);
+ Multiply((nuint)Pow10Table[(int)exponent]);
}
else if (!IsZero())
{
@@ -1197,6 +1118,18 @@ public static void SetUInt64(out BigInteger result, ulong value)
{
Unsafe.SkipInit(out result);
+#if TARGET_64BIT
+ // A ulong fits into a single 64-bit block.
+ if (value == 0)
+ {
+ SetZero(out result);
+ }
+ else
+ {
+ result._blocks[0] = (nuint)value;
+ result._length = 1;
+ }
+#else
if (value <= uint.MaxValue)
{
SetUInt32(out result, (uint)value);
@@ -1208,15 +1141,16 @@ public static void SetUInt64(out BigInteger result, ulong value)
result._length = 2;
}
+#endif
}
- public static void SetValue(out BigInteger result, scoped ref BigInteger value)
+ public static void SetValue(out BigInteger result, scoped ref readonly BigInteger value)
{
Unsafe.SkipInit(out result);
- int rhsLength = value._length;
+ int rhsLength = (int)value._length;
result._length = rhsLength;
- Buffer.Memmove(ref result._blocks[0], ref value._blocks[0], (nuint)rhsLength);
+ value._blocks[..rhsLength].CopyTo(result._blocks);
}
public static void SetZero(out BigInteger result)
@@ -1230,14 +1164,14 @@ public void ShiftLeft(int shift)
Debug.Assert(shift >= 0);
// Process blocks high to low so that we can safely process in place
- int length = _length;
+ int length = (int)_length;
if ((length == 0) || (shift == 0))
{
return;
}
- int blocksToShift = DivRem32(shift, out int remainingBitsToShift);
+ int blocksToShift = DivRemBlock(shift, out int remainingBitsToShift);
// Copy blocks from high to low
int readIndex = length - 1;
@@ -1258,12 +1192,8 @@ public void ShiftLeft(int shift)
return;
}
- while (readIndex >= 0)
- {
- _blocks[writeIndex] = _blocks[readIndex];
- readIndex--;
- writeIndex--;
- }
+ // CopyTo is overlap-safe, so it handles the in-place upward block shift.
+ ((Span)_blocks).Slice(0, length).CopyTo(((Span)_blocks).Slice(blocksToShift));
_length += blocksToShift;
@@ -1291,11 +1221,11 @@ public void ShiftLeft(int shift)
_length = writeIndex + 1;
// Output the initial blocks
- int lowBitsShift = 32 - remainingBitsToShift;
+ int lowBitsShift = BitsPerBlock - remainingBitsToShift;
- uint highBits = 0;
- uint block = _blocks[readIndex];
- uint lowBits = block >> lowBitsShift;
+ nuint highBits = 0;
+ nuint block = _blocks[readIndex];
+ nuint lowBits = block >> lowBitsShift;
while (readIndex > 0)
{
@@ -1317,25 +1247,24 @@ public void ShiftLeft(int shift)
Clear(blocksToShift);
// Check if the terminating block has no set bits
- if (_blocks[_length - 1] == 0)
+ if (_blocks[(int)_length - 1] == 0)
{
_length--;
}
}
}
- public readonly uint ToUInt32()
+ public readonly ulong ToUInt64()
{
+#if TARGET_64BIT
+ // A single 64-bit block already holds the full value.
if (_length > 0)
{
return _blocks[0];
}
return 0;
- }
-
- public readonly ulong ToUInt64()
- {
+#else
if (_length > 1)
{
return ((ulong)_blocks[1] << 32) + _blocks[0];
@@ -1347,10 +1276,24 @@ public readonly ulong ToUInt64()
}
return 0;
+#endif
}
public UInt128 ToUInt128()
{
+#if TARGET_64BIT
+ if (_length > 1)
+ {
+ return new UInt128(_blocks[1], _blocks[0]);
+ }
+
+ if (_length > 0)
+ {
+ return _blocks[0];
+ }
+
+ return 0;
+#else
if (_length > 3)
{
return new UInt128(((ulong)_blocks[3] << 32) + _blocks[2], ((ulong)(_blocks[1]) << 32) + _blocks[0]);
@@ -1372,21 +1315,47 @@ public UInt128 ToUInt128()
}
return 0;
+#endif
+ }
+
+ private readonly nuint ToBlock()
+ {
+ if (_length > 0)
+ {
+ return _blocks[0];
+ }
+
+ return 0;
+ }
+
+ private static void SetBlock(out BigInteger result, nuint value)
+ {
+ Unsafe.SkipInit(out result);
+
+ if (value == 0)
+ {
+ result._length = 0;
+ }
+ else
+ {
+ result._blocks[0] = value;
+ result._length = 1;
+ }
}
- private void Clear(int length) => ((Span)_blocks).Slice(0, length).Clear();
+ private void Clear(int length) => ((Span)_blocks).Slice(0, length).Clear();
- private static int DivRem32(int value, out int remainder)
+ private static int DivRemBlock(int value, out int remainder)
{
Debug.Assert(value >= 0);
- remainder = value & 31;
- return value >>> 5;
+ remainder = value & (BitsPerBlock - 1);
+ return value >>> BlockShift;
}
[InlineArray(MaxBlockCount)]
private struct BlocksBuffer
{
- public uint e0;
+ public nuint e0;
}
}
}
diff --git a/src/libraries/System.Private.CoreLib/src/System/Number.Dragon4.cs b/src/libraries/System.Private.CoreLib/src/System/Number.Dragon4.cs
index c097abe3ab58b5..c829e84423cc54 100644
--- a/src/libraries/System.Private.CoreLib/src/System/Number.Dragon4.cs
+++ b/src/libraries/System.Private.CoreLib/src/System/Number.Dragon4.cs
@@ -276,23 +276,23 @@ private static uint Dragon4(ulong mantissa, int exponent, uint mantissaHighBitId
// Output the exponent of the first digit we will print
decimalExponent = --digitExponent;
- // In preparation for calling BigInteger.HeuristicDivie(), we need to scale up our values such that the highest block of the denominator is greater than or equal to 8.
+ // In preparation for calling BigInteger.HeuristicDivide(), we need to scale up our values such that the highest block of the denominator is greater than or equal to 8.
// We also need to guarantee that the numerator can never have a length greater than the denominator after each loop iteration.
- // This requires the highest block of the denominator to be less than or equal to 429496729 which is the highest number that can be multiplied by 10 without overflowing to a new block.
+ // This requires the highest block of the denominator to be less than or equal to the highest number that can be multiplied by 10 without overflowing to a new block (nuint.MaxValue / 10).
Debug.Assert(scale.GetLength() > 0);
- uint hiBlock = scale.GetBlock(scale.GetLength() - 1);
+ nuint hiBlock = scale.GetBlock(scale.GetLength() - 1);
- if ((hiBlock < 8) || (hiBlock > 429496729))
+ if ((hiBlock < 8) || (hiBlock > (nuint.MaxValue / 10)))
{
- // Perform a bit shift on all values to get the highest block of the denominator into the range [8,429496729].
- // We are more likely to make accurate quotient estimations in BigInteger.HeuristicDivide() with higher denominator values so we shift the denominator to place the highest bit at index 27 of the highest block.
- // This is safe because (2^28 - 1) = 268435455 which is less than 429496729.
- // This means that all values with a highest bit at index 27 are within range.
+ // Perform a bit shift on all values to get the highest block of the denominator into the range [8, nuint.MaxValue / 10].
+ // We are more likely to make accurate quotient estimations in BigInteger.HeuristicDivide() with higher denominator values so we shift the denominator to place the highest bit at index (BitsPerBlock - 5) of the highest block.
+ // This is safe because (2^(BitsPerBlock - 4) - 1) is less than (nuint.MaxValue / 10).
+ // This means that all values with a highest bit at index (BitsPerBlock - 5) are within range.
Debug.Assert(hiBlock != 0);
int hiBlockLog2 = BitOperations.Log2(hiBlock);
- Debug.Assert((hiBlockLog2 < 3) || (hiBlockLog2 > 27));
- int shift = (32 + 27 - hiBlockLog2) % 32;
+ Debug.Assert((hiBlockLog2 < 3) || (hiBlockLog2 > (BigInteger.BitsPerBlock - 5)));
+ int shift = (BigInteger.BitsPerBlock + (BigInteger.BitsPerBlock - 5) - hiBlockLog2) % BigInteger.BitsPerBlock;
scale.ShiftLeft(shift);
scaledValue.ShiftLeft(shift);
@@ -307,7 +307,7 @@ private static uint Dragon4(ulong mantissa, int exponent, uint mantissaHighBitId
// These values are used to inspect why the print loop terminated so we can properly round the final digit.
bool low; // did the value get within marginLow distance from zero
bool high; // did the value get within marginHigh distance from one
- uint outputDigit; // current digit being output
+ uint outputDigit = 0; // current digit being output
if (cutoffNumber == -1)
{
@@ -370,24 +370,94 @@ private static uint Dragon4(ulong mantissa, int exponent, uint mantissaHighBitId
low = false;
high = false;
- while (true)
+ // This path never inspects per-digit margins, so when the per-digit divide is genuinely
+ // O(L > 1) we pull a whole native block's worth of digits at a time: `floor(scaledValue *
+ // 10^(K-1) / scale)` yields the next K digits in one big divide. K is the largest power of
+ // ten that fits a block (19 on 64-bit, 9 on 32-bit). The digit at cutoffExponent is always
+ // the rounding digit, so it is left to a per-digit tail whose leftover scaledValue drives
+ // the rounding logic. A single-block scale (e.g. small integers) already extracts a digit
+ // in O(1), so batching would only add block-divide overhead -- it stays on the scalar loop.
+ if (scale.GetLength() > 1)
{
- // divide out the scale to extract the digit
- outputDigit = BigInteger.HeuristicDivide(ref scaledValue, ref scale);
- Debug.Assert(outputDigit < 10);
+ int maxBatchDigits = (nint.Size == 8) ? 19 : 9;
+ bool roundingDigitExtracted = false;
- if (scaledValue.IsZero() || (digitExponent <= cutoffExponent))
+ while (digitExponent > cutoffExponent)
{
- break;
+ int batchDigits = Math.Min(digitExponent - cutoffExponent, maxBatchDigits);
+
+ if (batchDigits > 1)
+ {
+ scaledValue.MultiplyPow10((uint)(batchDigits - 1));
+ }
+
+ BigInteger.DivRem(ref scaledValue, ref scale, out BigInteger batch, out BigInteger remainder);
+ Debug.Assert(batch.GetLength() <= 1);
+
+ // Expand the block quotient into `batchDigits` zero-padded decimal digits.
+ nuint blockDigits = (batch.GetLength() == 0) ? 0 : batch.GetBlock(0);
+
+ for (int i = batchDigits - 1; i >= 0; i--)
+ {
+ (blockDigits, nuint digit) = Math.DivRem(blockDigits, (nuint)10);
+ buffer[curDigit + i] = (byte)('0' + digit);
+ }
+
+ if (remainder.IsZero())
+ {
+ // The value was captured exactly within this block. The trailing zeros are not
+ // significant, so the last non-zero digit becomes the rounding digit and the
+ // leftover scaledValue is exactly zero (an exact round-down below).
+ int last = curDigit + batchDigits - 1;
+
+ while ((last > curDigit) && (buffer[last] == '0'))
+ {
+ last--;
+ }
+
+ outputDigit = (uint)(buffer[last] - '0');
+ curDigit = last;
+ BigInteger.SetZero(out scaledValue);
+ roundingDigitExtracted = true;
+ break;
+ }
+
+ curDigit += batchDigits;
+
+ // multiply larger by the output base to line up the next block
+ remainder.Multiply10();
+ scaledValue = remainder;
+ digitExponent -= batchDigits;
}
- // store the output digit
- buffer[curDigit] = (byte)('0' + outputDigit);
- curDigit++;
+ if (!roundingDigitExtracted)
+ {
+ // extract the rounding digit at cutoffExponent
+ outputDigit = BigInteger.HeuristicDivide(ref scaledValue, ref scale);
+ Debug.Assert(outputDigit < 10);
+ }
+ }
+ else
+ {
+ while (true)
+ {
+ // divide out the scale to extract the digit
+ outputDigit = BigInteger.HeuristicDivide(ref scaledValue, ref scale);
+ Debug.Assert(outputDigit < 10);
- // multiply larger by the output base
- scaledValue.Multiply10();
- digitExponent--;
+ if (scaledValue.IsZero() || (digitExponent <= cutoffExponent))
+ {
+ break;
+ }
+
+ // store the output digit
+ buffer[curDigit] = (byte)('0' + outputDigit);
+ curDigit++;
+
+ // multiply larger by the output base
+ scaledValue.Multiply10();
+ digitExponent--;
+ }
}
}
else
diff --git a/src/libraries/System.Private.CoreLib/src/System/Number.NumberToFloatingPointBits.cs b/src/libraries/System.Private.CoreLib/src/System/Number.NumberToFloatingPointBits.cs
index c4dca930cb5d0d..86fe66de0d2da3 100644
--- a/src/libraries/System.Private.CoreLib/src/System/Number.NumberToFloatingPointBits.cs
+++ b/src/libraries/System.Private.CoreLib/src/System/Number.NumberToFloatingPointBits.cs
@@ -706,8 +706,22 @@ internal static void AccumulateDecimalDigitsIntoBigInteger(scoped ref NumberBuff
while (remaining != 0)
{
- uint count = Math.Min(remaining, 9);
- uint value = DigitsToUInt32(src, (int)(count));
+ // Batch as many digits as fill a single block -- 9 on 32-bit (10^9 fits a uint),
+ // 19 on 64-bit (10^19 fits a nuint) -- so wide builds halve the multiply/add
+ // iterations instead of staying 32-bit-granular. nint.Size constant-folds here.
+ uint count;
+ nuint value;
+
+ if (nint.Size == 8)
+ {
+ count = Math.Min(remaining, 19);
+ value = (nuint)DigitsToUInt64(src, (int)count);
+ }
+ else
+ {
+ count = Math.Min(remaining, 9);
+ value = DigitsToUInt32(src, (int)count);
+ }
result.MultiplyPow10(count);
result.Add(value);
@@ -847,45 +861,12 @@ private static ulong ConvertBigIntegerToFloatingPointBits(ref BigInteger
return AssembleFloatingPointBits(value.ToUInt64(), baseExponent, !hasNonZeroFractionalPart);
}
- (int topBlockIndex, int topBlockBits) = Math.DivRem(integerBitsOfPrecision, 32);
- int middleBlockIndex = topBlockIndex - 1;
- int bottomBlockIndex = middleBlockIndex - 1;
+ // The mantissa is the top 64 bits of the value; everything below that window is the tail.
+ int windowBitIndex = integerBitsOfPrecision - 64;
- ulong mantissa;
- int exponent = baseExponent + ((int)(bottomBlockIndex) * 32);
- bool hasZeroTail = !hasNonZeroFractionalPart;
-
- // When the top 64-bits perfectly span two blocks, we can get those blocks directly
- if (topBlockBits == 0)
- {
- mantissa = ((ulong)(value.GetBlock(middleBlockIndex)) << 32) + value.GetBlock(bottomBlockIndex);
- }
- else
- {
- // Otherwise, we need to read three blocks and combine them into a 64-bit mantissa
-
- int bottomBlockShift = (int)(topBlockBits);
- int topBlockShift = 64 - bottomBlockShift;
- int middleBlockShift = topBlockShift - 32;
-
- exponent += (int)(topBlockBits);
-
- uint bottomBlock = value.GetBlock(bottomBlockIndex);
- uint bottomBits = bottomBlock >> bottomBlockShift;
-
- ulong middleBits = (ulong)(value.GetBlock(middleBlockIndex)) << middleBlockShift;
- ulong topBits = (ulong)(value.GetBlock(topBlockIndex)) << topBlockShift;
-
- mantissa = topBits + middleBits + bottomBits;
-
- uint unusedBottomBlockBitsMask = (1u << (int)(topBlockBits)) - 1;
- hasZeroTail &= (bottomBlock & unusedBottomBlockBitsMask) == 0;
- }
-
- for (int i = 0; i < bottomBlockIndex; i++)
- {
- hasZeroTail &= (value.GetBlock(i) == 0);
- }
+ ulong mantissa = value.GetBits64(windowBitIndex);
+ int exponent = baseExponent + windowBitIndex;
+ bool hasZeroTail = !hasNonZeroFractionalPart && value.HasZeroTail(windowBitIndex);
return AssembleFloatingPointBits(mantissa, exponent, hasZeroTail);
}
diff --git a/src/libraries/System.Private.CoreLib/src/System/Number.Rounding.cs b/src/libraries/System.Private.CoreLib/src/System/Number.Rounding.cs
index c6efa31b6c5f29..5e258928718c9a 100644
--- a/src/libraries/System.Private.CoreLib/src/System/Number.Rounding.cs
+++ b/src/libraries/System.Private.CoreLib/src/System/Number.Rounding.cs
@@ -404,7 +404,7 @@ public static unsafe TNumber RoundToDecimalDigits(TNumber value, int di
do
{
BigInteger.DivRem(ref quotient, ref ten, out BigInteger newQuotient, out BigInteger digit);
- uint digitValue = digit.IsZero() ? 0 : digit.GetBlock(0);
+ uint digitValue = digit.IsZero() ? 0 : (uint)digit.GetBlock(0);
buffer[digitCount++] = (byte)('0' + digitValue);
BigInteger.SetValue(out quotient, ref newQuotient);
}
diff --git a/src/libraries/System.Private.CoreLib/src/System/Numerics/BigIntegerCalculator.Shared.cs b/src/libraries/System.Private.CoreLib/src/System/Numerics/BigIntegerCalculator.Shared.cs
new file mode 100644
index 00000000000000..f5bec8d2c6f5f8
--- /dev/null
+++ b/src/libraries/System.Private.CoreLib/src/System/Numerics/BigIntegerCalculator.Shared.cs
@@ -0,0 +1,747 @@
+// Licensed to the .NET Foundation under one or more agreements.
+// The .NET Foundation licenses this file to you under the MIT license.
+
+using System.Diagnostics;
+using System.Runtime.CompilerServices;
+using System.Runtime.Intrinsics.X86;
+
+namespace System.Numerics
+{
+ // Magnitude (unsigned, native-limb) arithmetic shared between the public
+ // System.Numerics.BigInteger and the internal System.Number.BigInteger used by
+ // floating-point parsing/formatting/rounding. Both are nuint-backed, so these kernels
+ // operate purely on Span/ReadOnlySpan and are free of any sign, allocation,
+ // or type-specific policy. This file is compiled into System.Private.CoreLib and linked
+ // into System.Runtime.Numerics.
+ internal static partial class BigIntegerCalculator
+ {
+ /// Number of bits per native-width limb: 32 on 32-bit, 64 on 64-bit.
+ internal static int BitsPerLimb
+ {
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ get => nint.Size * 8;
+ }
+
+ public static int Compare(ReadOnlySpan left, ReadOnlySpan right)
+ {
+ Debug.Assert(left.Length <= right.Length || left.Slice(right.Length).ContainsAnyExcept(0u));
+ Debug.Assert(left.Length >= right.Length || right.Slice(left.Length).ContainsAnyExcept(0u));
+
+ if (left.Length != right.Length)
+ {
+ return left.Length < right.Length ? -1 : 1;
+ }
+
+ int iv = left.Length;
+ while (--iv >= 0 && left[iv] == right[iv]) ;
+
+ if (iv < 0)
+ {
+ return 0;
+ }
+
+ 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
+ // of a given value may be less than the array's length
+
+ return value.LastIndexOfAnyExcept((nuint)0) + 1;
+ }
+
+ ///
+ /// Performs widening addition of two limbs plus a carry-in, returning the sum and carry-out.
+ /// On 64-bit: uses 128-bit arithmetic. On 32-bit: uses 64-bit arithmetic.
+ ///
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ internal static nuint AddWithCarry(nuint a, nuint b, nuint carryIn, out nuint carryOut)
+ {
+ if (nint.Size == 8)
+ {
+ nuint sum1 = a + b;
+ nuint c1 = (sum1 < a) ? 1 : (nuint)0;
+ nuint sum2 = sum1 + carryIn;
+ nuint c2 = (sum2 < sum1) ? 1 : (nuint)0;
+ carryOut = c1 + c2;
+ return sum2;
+ }
+ else
+ {
+ ulong sum = (ulong)a + b + carryIn;
+ carryOut = (uint)(sum >> 32);
+ return (uint)sum;
+ }
+ }
+
+ ///
+ /// Performs widening subtraction of two limbs with a borrow-in, returning the difference and borrow-out.
+ /// borrowOut is 0 (no borrow) or 1 (borrow occurred).
+ ///
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ internal static nuint SubWithBorrow(nuint a, nuint b, nuint borrowIn, out nuint borrowOut)
+ {
+ if (nint.Size == 8)
+ {
+ // Use unsigned underflow detection
+ nuint diff1 = a - b;
+ nuint b1 = (diff1 > a) ? 1 : (nuint)0;
+ nuint diff2 = diff1 - borrowIn;
+ nuint b2 = (diff2 > diff1) ? 1 : (nuint)0;
+ borrowOut = b1 + b2;
+ return diff2;
+ }
+ else
+ {
+ long diff = (long)a - (long)b - (long)borrowIn;
+ borrowOut = (uint)(-(int)(diff >> 32)); // 0 or 1
+ return (uint)diff;
+ }
+ }
+
+ ///
+ /// Widening divide: (hi:lo) / divisor -> (quotient, remainder).
+ ///
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ internal static nuint DivRem(nuint hi, nuint lo, nuint divisor, out nuint remainder)
+ {
+ if (nint.Size == 8)
+ {
+ // Compute (hi * 2^64 + lo) / divisor.
+ // hi < divisor is guaranteed by callers, so quotient fits in 64 bits.
+ Debug.Assert(hi < (ulong)divisor || divisor == 0);
+
+ if (hi == 0)
+ {
+ (ulong q, ulong r) = Math.DivRem(lo, (ulong)divisor);
+ remainder = (nuint)r;
+ return (nuint)q;
+ }
+
+ // When divisor fits in 32 bits, split lo into two 32-bit halves
+ // and chain two native 64-bit divisions (avoids UInt128 overhead):
+ // (hi * 2^32 + lo_hi) / divisor -> (q_hi, r1) [fits: hi < divisor < 2^32]
+ // (r1 * 2^32 + lo_lo) / divisor -> (q_lo, r2) [fits: r1 < divisor < 2^32]
+ if ((ulong)divisor <= uint.MaxValue)
+ {
+ ulong lo_hi = (ulong)lo >> 32;
+ ulong lo_lo = (ulong)lo & 0xFFFFFFFF;
+
+ (ulong q_hi, ulong r1) = Math.DivRem(((ulong)hi << 32) | lo_hi, divisor);
+ (ulong q_lo, ulong r2) = Math.DivRem((r1 << 32) | lo_lo, divisor);
+
+ remainder = (nuint)r2;
+ return (nuint)((q_hi << 32) | q_lo);
+ }
+
+ {
+#pragma warning disable SYSLIB5004 // X86Base.DivRem is experimental
+ if (X86Base.X64.IsSupported)
+ {
+ (ulong q, ulong r) = X86Base.X64.DivRem(lo, hi, divisor);
+ remainder = (nuint)r;
+ return (nuint)q;
+ }
+#pragma warning restore SYSLIB5004
+
+ UInt128 value = ((UInt128)(ulong)hi << 64) | (ulong)lo;
+ UInt128 digit = value / (ulong)divisor;
+ remainder = (nuint)(ulong)(value - digit * (ulong)divisor);
+ return (nuint)(ulong)digit;
+ }
+ }
+ else
+ {
+ ulong value = ((ulong)hi << 32) | lo;
+ ulong digit = value / divisor;
+ remainder = (uint)(value - digit * divisor);
+ return (uint)digit;
+ }
+ }
+
+ ///
+ /// Multiply by scalar: result[0..left.Length] = left * multiplier.
+ /// Returns the carry out. Unrolled by 4 on 64-bit.
+ /// Unlike MulAdd1, this writes to result rather than accumulating.
+ ///
+ internal static nuint Mul1(Span result, ReadOnlySpan left, nuint multiplier)
+ {
+ Debug.Assert(result.Length >= left.Length);
+
+ int length = left.Length;
+ int i = 0;
+ nuint carry = 0;
+
+ if (nint.Size == 8)
+ {
+ for (; i + 3 < length; i += 4)
+ {
+ UInt128 p0 = (UInt128)(ulong)left[i] * (ulong)multiplier + (ulong)carry;
+ result[i] = (nuint)(ulong)p0;
+
+ UInt128 p1 = (UInt128)(ulong)left[i + 1] * (ulong)multiplier + (ulong)(p0 >> 64);
+ result[i + 1] = (nuint)(ulong)p1;
+
+ UInt128 p2 = (UInt128)(ulong)left[i + 2] * (ulong)multiplier + (ulong)(p1 >> 64);
+ result[i + 2] = (nuint)(ulong)p2;
+
+ UInt128 p3 = (UInt128)(ulong)left[i + 3] * (ulong)multiplier + (ulong)(p2 >> 64);
+ result[i + 3] = (nuint)(ulong)p3;
+
+ carry = (nuint)(ulong)(p3 >> 64);
+ }
+
+ for (; i < length; i++)
+ {
+ UInt128 product = (UInt128)(ulong)left[i] * (ulong)multiplier + (ulong)carry;
+ result[i] = (nuint)(ulong)product;
+ carry = (nuint)(ulong)(product >> 64);
+ }
+ }
+ else
+ {
+ for (; i < length; i++)
+ {
+ ulong product = (ulong)left[i] * multiplier + carry;
+ result[i] = (uint)product;
+ carry = (uint)(product >> 32);
+ }
+ }
+
+ return carry;
+ }
+
+ ///
+ /// Fused multiply-accumulate by scalar: result[0..left.Length] += left * multiplier.
+ /// Returns the carry out. Unrolled by 4 on 64-bit to overlap multiply latencies.
+ ///
+ internal static nuint MulAdd1(Span result, ReadOnlySpan left, nuint multiplier)
+ {
+ Debug.Assert(result.Length >= left.Length);
+
+ int length = left.Length;
+ int i = 0;
+ nuint carry = 0;
+
+ if (nint.Size == 8)
+ {
+ // Unroll by 4: mulx has 3-5 cycle latency but 1 cycle throughput,
+ // so issuing 4 multiplies allows the CPU to pipeline them while
+ // carry chains complete sequentially behind.
+ for (; i + 3 < length; i += 4)
+ {
+ UInt128 p0 = (UInt128)(ulong)left[i] * (ulong)multiplier + (ulong)result[i] + (ulong)carry;
+ result[i] = (nuint)(ulong)p0;
+
+ UInt128 p1 = (UInt128)(ulong)left[i + 1] * (ulong)multiplier + (ulong)result[i + 1] + (ulong)(p0 >> 64);
+ result[i + 1] = (nuint)(ulong)p1;
+
+ UInt128 p2 = (UInt128)(ulong)left[i + 2] * (ulong)multiplier + (ulong)result[i + 2] + (ulong)(p1 >> 64);
+ result[i + 2] = (nuint)(ulong)p2;
+
+ UInt128 p3 = (UInt128)(ulong)left[i + 3] * (ulong)multiplier + (ulong)result[i + 3] + (ulong)(p2 >> 64);
+ result[i + 3] = (nuint)(ulong)p3;
+
+ carry = (nuint)(ulong)(p3 >> 64);
+ }
+
+ for (; i < length; i++)
+ {
+ UInt128 product = (UInt128)(ulong)left[i] * (ulong)multiplier + (ulong)result[i] + (ulong)carry;
+ result[i] = (nuint)(ulong)product;
+ carry = (nuint)(ulong)(product >> 64);
+ }
+ }
+ else
+ {
+ for (; i < length; i++)
+ {
+ ulong product = (ulong)left[i] * multiplier
+ + result[i] + carry;
+ result[i] = (uint)product;
+ carry = (uint)(product >> 32);
+ }
+ }
+
+ return carry;
+ }
+
+ ///
+ /// Fused subtract-multiply by scalar: result[0..right.Length] -= right * multiplier.
+ /// Returns the borrow out. Unrolled by 4 on 64-bit.
+ ///
+ internal static nuint SubMul1(Span result, ReadOnlySpan right, nuint multiplier)
+ {
+ Debug.Assert(result.Length >= right.Length);
+
+ int length = right.Length;
+ int i = 0;
+ nuint carry = 0;
+
+ if (nint.Size == 8)
+ {
+ for (; i + 3 < length; i += 4)
+ {
+ UInt128 prod0 = (UInt128)(ulong)right[i] * (ulong)multiplier + (ulong)carry;
+ nuint lo0 = (nuint)(ulong)prod0;
+ nuint hi0 = (nuint)(ulong)(prod0 >> 64);
+ nuint orig0 = result[i];
+ result[i] = orig0 - lo0;
+ hi0 += (orig0 < lo0) ? (nuint)1 : 0;
+
+ UInt128 prod1 = (UInt128)(ulong)right[i + 1] * (ulong)multiplier + (ulong)hi0;
+ nuint lo1 = (nuint)(ulong)prod1;
+ nuint hi1 = (nuint)(ulong)(prod1 >> 64);
+ nuint orig1 = result[i + 1];
+ result[i + 1] = orig1 - lo1;
+ hi1 += (orig1 < lo1) ? (nuint)1 : 0;
+
+ UInt128 prod2 = (UInt128)(ulong)right[i + 2] * (ulong)multiplier + (ulong)hi1;
+ nuint lo2 = (nuint)(ulong)prod2;
+ nuint hi2 = (nuint)(ulong)(prod2 >> 64);
+ nuint orig2 = result[i + 2];
+ result[i + 2] = orig2 - lo2;
+ hi2 += (orig2 < lo2) ? (nuint)1 : 0;
+
+ UInt128 prod3 = (UInt128)(ulong)right[i + 3] * (ulong)multiplier + (ulong)hi2;
+ nuint lo3 = (nuint)(ulong)prod3;
+ nuint hi3 = (nuint)(ulong)(prod3 >> 64);
+ nuint orig3 = result[i + 3];
+ result[i + 3] = orig3 - lo3;
+ hi3 += (orig3 < lo3) ? (nuint)1 : 0;
+
+ carry = hi3;
+ }
+
+ for (; i < length; i++)
+ {
+ UInt128 product = (UInt128)(ulong)right[i] * (ulong)multiplier + (ulong)carry;
+ nuint lo = (nuint)(ulong)product;
+ nuint hi = (nuint)(ulong)(product >> 64);
+ nuint orig = result[i];
+ result[i] = orig - lo;
+ hi += (orig < lo) ? (nuint)1 : 0;
+ carry = hi;
+ }
+ }
+ else
+ {
+ for (; i < length; i++)
+ {
+ ulong product = (ulong)right[i] * multiplier + carry;
+ uint lo = (uint)product;
+ uint hi = (uint)(product >> 32);
+
+ uint orig = (uint)result[i];
+ result[i] = orig - lo;
+ hi += (orig < lo) ? 1u : 0;
+
+ carry = hi;
+ }
+ }
+
+ return carry;
+ }
+
+ private const int CopyToThreshold = 8;
+
+ private static void CopyTail(ReadOnlySpan source, Span dest, int start)
+ {
+ source.Slice(start).CopyTo(dest.Slice(start));
+ }
+
+ public static void Add(ReadOnlySpan left, nuint right, Span bits)
+ {
+ Debug.Assert(left.Length >= 1);
+ Debug.Assert(bits.Length == left.Length + 1);
+
+ Add(left, bits, startIndex: 0, initialCarry: right);
+ }
+
+ public static void Add(ReadOnlySpan left, ReadOnlySpan right, Span bits)
+ {
+ Debug.Assert(right.Length >= 1);
+ Debug.Assert(left.Length >= right.Length);
+ Debug.Assert(bits.Length == left.Length + 1);
+
+ // Establish cross-span length relationships so the JIT can
+ // elide bounds checks for left[i] and bits[i] in the loop.
+ _ = left[right.Length - 1];
+ _ = bits[right.Length];
+
+ nuint carry = 0;
+
+ for (int i = 0; i < right.Length; i++)
+ {
+ bits[i] = AddWithCarry(left[i], right[i], carry, out carry);
+ }
+
+ Add(left, bits, startIndex: right.Length, initialCarry: carry);
+ }
+
+ public static void AddSelf(Span left, ReadOnlySpan right)
+ {
+ Debug.Assert(left.Length >= right.Length);
+
+ int i = 0;
+ nuint carry = 0;
+
+ if (right.Length != 0)
+ {
+ _ = left[right.Length - 1];
+ }
+
+ for (; i < right.Length; i++)
+ {
+ left[i] = AddWithCarry(left[i], right[i], carry, out carry);
+ }
+
+ for (; carry != 0 && i < left.Length; i++)
+ {
+ nuint sum = left[i] + carry;
+ carry = (sum < carry) ? 1 : (nuint)0;
+ left[i] = sum;
+ }
+
+ Debug.Assert(carry == 0);
+ }
+
+ public static void Subtract(ReadOnlySpan left, nuint right, Span bits)
+ {
+ Debug.Assert(left.Length >= 1);
+ Debug.Assert(left[0] >= right || left.Length >= 2);
+ Debug.Assert(bits.Length == left.Length);
+
+ Subtract(left, bits, startIndex: 0, initialBorrow: right);
+ }
+
+ public static void Subtract(ReadOnlySpan left, ReadOnlySpan right, Span bits)
+ {
+ Debug.Assert(right.Length >= 1);
+ Debug.Assert(left.Length >= right.Length);
+ Debug.Assert(CompareActual(left, right) >= 0);
+ Debug.Assert(bits.Length == left.Length);
+
+ _ = left[right.Length - 1];
+ _ = bits[right.Length - 1];
+
+ nuint borrow = 0;
+
+ for (int i = 0; i < right.Length; i++)
+ {
+ bits[i] = SubWithBorrow(left[i], right[i], borrow, out borrow);
+ }
+
+ Subtract(left, bits, startIndex: right.Length, initialBorrow: borrow);
+ }
+
+ public static void SubtractSelf(Span left, ReadOnlySpan right)
+ {
+ Debug.Assert(left.Length >= right.Length);
+
+ int i = 0;
+ nuint borrow = 0;
+
+ if (right.Length != 0)
+ {
+ _ = left[right.Length - 1];
+ }
+
+ for (; i < right.Length; i++)
+ {
+ left[i] = SubWithBorrow(left[i], right[i], borrow, out borrow);
+ }
+
+ for (; borrow != 0 && i < left.Length; i++)
+ {
+ nuint val = left[i];
+ left[i] = val - borrow;
+ borrow = val == 0 ? 1 : (nuint)0;
+ }
+
+ Debug.Assert(borrow == 0);
+ }
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ private static void Add(ReadOnlySpan left, Span bits, int startIndex, nuint initialCarry)
+ {
+ // Executes the addition for one big and one single-limb integer.
+
+ int i = startIndex;
+ nuint carry = initialCarry;
+
+ _ = bits[left.Length];
+
+ if (left.Length <= CopyToThreshold)
+ {
+ for (; i < left.Length; i++)
+ {
+ nuint sum = left[i] + carry;
+ carry = (sum < carry) ? 1 : (nuint)0;
+ bits[i] = sum;
+ }
+
+ bits[left.Length] = carry;
+ }
+ else
+ {
+ for (; i < left.Length;)
+ {
+ nuint sum = left[i] + carry;
+ carry = (sum < carry) ? 1 : (nuint)0;
+ bits[i] = sum;
+ i++;
+
+ // Once carry is set to 0 it can not be 1 anymore.
+ // So the tail of the loop is just the movement of argument values to result span.
+ if (carry == 0)
+ {
+ break;
+ }
+ }
+
+ bits[left.Length] = carry;
+
+ if (i < left.Length)
+ {
+ CopyTail(left, bits, i);
+ }
+ }
+ }
+
+ [MethodImpl(MethodImplOptions.AggressiveInlining)]
+ private static void Subtract(ReadOnlySpan left, Span bits, int startIndex, nuint initialBorrow)
+ {
+ // Executes the subtraction for one big and one single-limb integer.
+
+ int i = startIndex;
+ nuint borrow = initialBorrow;
+
+ if (left.Length != 0)
+ {
+ _ = bits[left.Length - 1];
+ }
+
+ if (left.Length <= CopyToThreshold)
+ {
+ for (; i < left.Length; i++)
+ {
+ nuint val = left[i];
+ nuint diff = val - borrow;
+ borrow = (diff > val) ? 1 : (nuint)0;
+ bits[i] = diff;
+ }
+ }
+ else
+ {
+ for (; i < left.Length;)
+ {
+ nuint val = left[i];
+ nuint diff = val - borrow;
+ borrow = (diff > val) ? 1 : (nuint)0;
+ bits[i] = diff;
+ i++;
+
+ // Once borrow is set to 0 it can not be 1 anymore.
+ // So the tail of the loop is just the movement of argument values to result span.
+ if (borrow == 0)
+ {
+ break;
+ }
+ }
+
+ if (i < left.Length)
+ {
+ CopyTail(left, bits, i);
+ }
+ }
+ }
+
+ ///
+ /// Multiply using the "grammar-school" method: bits += left * right[i] per limb of right.
+ /// Callers guarantee left.Length >= right.Length and bits large enough for the product.
+ ///
+ public static void MultiplyNaive(ReadOnlySpan left, ReadOnlySpan right, Span bits)
+ {
+ Debug.Assert(left.Length >= right.Length);
+ Debug.Assert(right.IsEmpty || bits.Length >= left.Length + right.Length);
+
+ // Multiplies the bits using the "grammar-school" method.
+ // Envisioning the "rhombus" of a pen-and-paper calculation
+ // should help getting the idea of these two loops...
+ // The inner multiplication operations are safe, because
+ // z_i+j + a_j * b_i + c <= 2(2^n - 1) + (2^n - 1)^2 =
+ // = 2^(2n) - 1, where n = BitsPerLimb.
+
+ for (int i = 0; i < right.Length; i++)
+ {
+ nuint carry = MulAdd1(bits.Slice(i), left, right[i]);
+ bits[i + left.Length] = carry;
+ }
+ }
+
+ internal 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
+ // block of the divisor. Thus, guessing digits of the quotient
+ // will be more precise. Additionally we'll get r = a % b.
+
+ nuint divHi = right[^1];
+ nuint divLo = right.Length > 1 ? right[^2] : 0;
+
+ // We measure the leading zeros of the divisor
+ int shift = (int)nuint.LeadingZeroCount(divHi);
+ int backShift = BitsPerLimb - shift;
+
+ // And, we make sure the most significant bit is set
+ if (shift > 0)
+ {
+ nuint divNx = right.Length > 2 ? right[^3] : 0;
+
+ divHi = (divHi << shift) | (divLo >> backShift);
+ divLo = (divLo << shift) | (divNx >> backShift);
+ }
+
+ // Then, we divide all of the bits as we would do it using
+ // pen and paper: guessing the next digit, subtracting, ...
+ for (int i = left.Length; i >= right.Length; i--)
+ {
+ int n = i - right.Length;
+ nuint t = (uint)i < (uint)left.Length ? left[i] : 0;
+
+ nuint valHi1 = t;
+ nuint valHi0 = left[i - 1];
+ nuint valLo = i > 1 ? left[i - 2] : 0;
+
+ // We shifted the divisor, we shift the dividend too
+ if (shift > 0)
+ {
+ nuint valNx = i > 2 ? left[i - 3] : 0;
+
+ valHi1 = (valHi1 << shift) | (valHi0 >> backShift);
+ valHi0 = (valHi0 << shift) | (valLo >> backShift);
+ valLo = (valLo << shift) | (valNx >> backShift);
+ }
+
+ // First guess for the current digit of the quotient,
+ // which naturally must have only native-width bits...
+ nuint digit = (valHi1 >= divHi) ? nuint.MaxValue : DivRem(valHi1, valHi0, divHi, out _);
+
+ // Our first guess may be a little bit too big
+ while (DivideGuessTooBig(digit, valHi1, valHi0, valLo, divHi, divLo))
+ {
+ --digit;
+ }
+
+ if (digit > 0)
+ {
+ // Now it's time to subtract our current quotient
+ nuint carry = SubtractDivisor(left.Slice(n), right, digit);
+ if (carry != t)
+ {
+ Debug.Assert(carry == t + 1);
+
+ // Our guess was still exactly one too high
+ carry = AddDivisor(left.Slice(n), right);
+ --digit;
+
+ Debug.Assert(carry == 1);
+ }
+ }
+
+ // We have the digit!
+ if ((uint)n < (uint)quotient.Length)
+ {
+ quotient[n] = digit;
+ }
+
+ if ((uint)i < (uint)left.Length)
+ {
+ left[i] = 0;
+ }
+ }
+ }
+
+ private static nuint AddDivisor(Span left, ReadOnlySpan right)
+ {
+ Debug.Assert(left.Length >= right.Length);
+
+ // Repairs the dividend, if the last subtract was too much
+
+ nuint carry = 0;
+
+ for (int i = 0; i < right.Length; i++)
+ {
+ ref nuint leftElement = ref left[i];
+ leftElement = AddWithCarry(leftElement, right[i], carry, out carry);
+ }
+
+ return carry;
+ }
+
+ private static nuint SubtractDivisor(Span left, ReadOnlySpan right, nuint q)
+ {
+ Debug.Assert(left.Length >= right.Length);
+
+ return SubMul1(left, right, q);
+ }
+
+ private static bool DivideGuessTooBig(nuint q, nuint valHi1, nuint valHi0,
+ nuint valLo, nuint divHi, nuint divLo)
+ {
+ // We multiply the two most significant limbs of the divisor
+ // with the current guess for the quotient. If those are bigger
+ // than the three most significant limbs of the current dividend
+ // we return true, which means the current guess is still too big.
+
+ nuint chkHiHi = nuint.BigMul(divHi, q, out nuint chkHiLo);
+ nuint chkLoHi = nuint.BigMul(divLo, q, out nuint chkLoLo);
+
+ chkHiLo += chkLoHi;
+ if (chkHiLo < chkLoHi)
+ {
+ chkHiHi++;
+ }
+
+ return (chkHiHi > valHi1)
+ || ((chkHiHi == valHi1) && ((chkHiLo > valHi0) || ((chkHiLo == valHi0) && (chkLoLo > valLo))));
+ }
+ }
+}
diff --git a/src/libraries/System.Private.CoreLib/src/System/UInt32.cs b/src/libraries/System.Private.CoreLib/src/System/UInt32.cs
index 59c5226f61599d..988838c154e5b4 100644
--- a/src/libraries/System.Private.CoreLib/src/System/UInt32.cs
+++ b/src/libraries/System.Private.CoreLib/src/System/UInt32.cs
@@ -304,7 +304,7 @@ public static uint Log10(uint value)
}
// Lookup table for power-of-10 boundaries corrections
- private static ReadOnlySpan PowersOf10 =>
+ internal static ReadOnlySpan PowersOf10 =>
[
1,
10,
diff --git a/src/libraries/System.Runtime.Numerics/src/System.Runtime.Numerics.csproj b/src/libraries/System.Runtime.Numerics/src/System.Runtime.Numerics.csproj
index 38a2f401fc7f36..96ba9ea061676a 100644
--- a/src/libraries/System.Runtime.Numerics/src/System.Runtime.Numerics.csproj
+++ b/src/libraries/System.Runtime.Numerics/src/System.Runtime.Numerics.csproj
@@ -12,7 +12,6 @@
-
@@ -34,6 +33,8 @@
+
- /// Specifies the minimum number of elements required to trigger a copy operation using an optimized path.
- ///
- ///
- /// This threshold is determined based on benchmarking and may be adjusted in the future to balance the overhead of copying versus the benefits of reduced loop iterations.
- ///
- private const int CopyToThreshold = 8;
-
- private static void CopyTail(ReadOnlySpan source, Span dest, int start)
- {
- source.Slice(start).CopyTo(dest.Slice(start));
- }
-
- public static void Add(ReadOnlySpan left, nuint right, Span bits)
- {
- Debug.Assert(left.Length >= 1);
- Debug.Assert(bits.Length == left.Length + 1);
-
- Add(left, bits, startIndex: 0, initialCarry: right);
- }
-
- public static void Add(ReadOnlySpan left, ReadOnlySpan right, Span bits)
- {
- Debug.Assert(right.Length >= 1);
- Debug.Assert(left.Length >= right.Length);
- Debug.Assert(bits.Length == left.Length + 1);
-
- // Establish cross-span length relationships so the JIT can
- // elide bounds checks for left[i] and bits[i] in the loop.
- _ = left[right.Length - 1];
- _ = bits[right.Length];
-
- nuint carry = 0;
-
- for (int i = 0; i < right.Length; i++)
- {
- bits[i] = AddWithCarry(left[i], right[i], carry, out carry);
- }
-
- Add(left, bits, startIndex: right.Length, initialCarry: carry);
- }
-
- public static void AddSelf(Span left, ReadOnlySpan right)
- {
- Debug.Assert(left.Length >= right.Length);
-
- int i = 0;
- nuint carry = 0;
-
- if (right.Length != 0)
- {
- _ = left[right.Length - 1];
- }
-
- for (; i < right.Length; i++)
- {
- left[i] = AddWithCarry(left[i], right[i], carry, out carry);
- }
-
- for (; carry != 0 && i < left.Length; i++)
- {
- nuint sum = left[i] + carry;
- carry = (sum < carry) ? 1 : (nuint)0;
- left[i] = sum;
- }
-
- Debug.Assert(carry == 0);
- }
-
- public static void Subtract(ReadOnlySpan left, nuint right, Span bits)
- {
- Debug.Assert(left.Length >= 1);
- Debug.Assert(left[0] >= right || left.Length >= 2);
- Debug.Assert(bits.Length == left.Length);
-
- Subtract(left, bits, startIndex: 0, initialBorrow: right);
- }
-
- public static void Subtract(ReadOnlySpan left, ReadOnlySpan right, Span bits)
- {
- Debug.Assert(right.Length >= 1);
- Debug.Assert(left.Length >= right.Length);
- Debug.Assert(CompareActual(left, right) >= 0);
- Debug.Assert(bits.Length == left.Length);
-
- _ = left[right.Length - 1];
- _ = bits[right.Length - 1];
-
- nuint borrow = 0;
-
- for (int i = 0; i < right.Length; i++)
- {
- bits[i] = SubWithBorrow(left[i], right[i], borrow, out borrow);
- }
-
- Subtract(left, bits, startIndex: right.Length, initialBorrow: borrow);
- }
-
- public static void SubtractSelf(Span left, ReadOnlySpan right)
- {
- Debug.Assert(left.Length >= right.Length);
-
- int i = 0;
- nuint borrow = 0;
-
- if (right.Length != 0)
- {
- _ = left[right.Length - 1];
- }
-
- for (; i < right.Length; i++)
- {
- left[i] = SubWithBorrow(left[i], right[i], borrow, out borrow);
- }
-
- for (; borrow != 0 && i < left.Length; i++)
- {
- nuint val = left[i];
- left[i] = val - borrow;
- borrow = val == 0 ? 1 : (nuint)0;
- }
-
- Debug.Assert(borrow == 0);
- }
-
- [MethodImpl(MethodImplOptions.AggressiveInlining)]
- private static void Add(ReadOnlySpan left, Span bits, int startIndex, nuint initialCarry)
- {
- // Executes the addition for one big and one single-limb integer.
-
- int i = startIndex;
- nuint carry = initialCarry;
-
- _ = bits[left.Length];
-
- if (left.Length <= CopyToThreshold)
- {
- for (; i < left.Length; i++)
- {
- nuint sum = left[i] + carry;
- carry = (sum < carry) ? 1 : (nuint)0;
- bits[i] = sum;
- }
-
- bits[left.Length] = carry;
- }
- else
- {
- for (; i < left.Length;)
- {
- nuint sum = left[i] + carry;
- carry = (sum < carry) ? 1 : (nuint)0;
- bits[i] = sum;
- i++;
-
- // Once carry is set to 0 it can not be 1 anymore.
- // So the tail of the loop is just the movement of argument values to result span.
- if (carry == 0)
- {
- break;
- }
- }
-
- bits[left.Length] = carry;
-
- if (i < left.Length)
- {
- CopyTail(left, bits, i);
- }
- }
- }
-
- [MethodImpl(MethodImplOptions.AggressiveInlining)]
- private static void Subtract(ReadOnlySpan left, Span bits, int startIndex, nuint initialBorrow)
- {
- // Executes the subtraction for one big and one single-limb integer.
-
- int i = startIndex;
- nuint borrow = initialBorrow;
-
- if (left.Length != 0)
- {
- _ = bits[left.Length - 1];
- }
-
- if (left.Length <= CopyToThreshold)
- {
- for (; i < left.Length; i++)
- {
- nuint val = left[i];
- nuint diff = val - borrow;
- borrow = (diff > val) ? 1 : (nuint)0;
- bits[i] = diff;
- }
- }
- else
- {
- for (; i < left.Length;)
- {
- nuint val = left[i];
- nuint diff = val - borrow;
- borrow = (diff > val) ? 1 : (nuint)0;
- bits[i] = diff;
- i++;
-
- // Once borrow is set to 0 it can not be 1 anymore.
- // So the tail of the loop is just the movement of argument values to result span.
- if (borrow == 0)
- {
- break;
- }
- }
-
- if (i < left.Length)
- {
- CopyTail(left, bits, i);
- }
- }
- }
- }
-}
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 119c241460e192..29aa27dff5f932 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
@@ -174,142 +174,6 @@ private static void DivRem(Span left, ReadOnlySpan right, 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
- // block of the divisor. Thus, guessing digits of the quotient
- // will be more precise. Additionally we'll get r = a % b.
-
- nuint divHi = right[^1];
- nuint divLo = right.Length > 1 ? right[^2] : 0;
-
- // We measure the leading zeros of the divisor
- int shift = (int)nuint.LeadingZeroCount(divHi);
- int backShift = BitsPerLimb - shift;
-
- // And, we make sure the most significant bit is set
- if (shift > 0)
- {
- nuint divNx = right.Length > 2 ? right[^3] : 0;
-
- divHi = (divHi << shift) | (divLo >> backShift);
- divLo = (divLo << shift) | (divNx >> backShift);
- }
-
- // Then, we divide all of the bits as we would do it using
- // pen and paper: guessing the next digit, subtracting, ...
- for (int i = left.Length; i >= right.Length; i--)
- {
- int n = i - right.Length;
- nuint t = (uint)i < (uint)left.Length ? left[i] : 0;
-
- nuint valHi1 = t;
- nuint valHi0 = left[i - 1];
- nuint valLo = i > 1 ? left[i - 2] : 0;
-
- // We shifted the divisor, we shift the dividend too
- if (shift > 0)
- {
- nuint valNx = i > 2 ? left[i - 3] : 0;
-
- valHi1 = (valHi1 << shift) | (valHi0 >> backShift);
- valHi0 = (valHi0 << shift) | (valLo >> backShift);
- valLo = (valLo << shift) | (valNx >> backShift);
- }
-
- // First guess for the current digit of the quotient,
- // which naturally must have only native-width bits...
- nuint digit = (valHi1 >= divHi) ? nuint.MaxValue : DivRem(valHi1, valHi0, divHi, out _);
-
- // Our first guess may be a little bit to big
- while (DivideGuessTooBig(digit, valHi1, valHi0, valLo, divHi, divLo))
- {
- --digit;
- }
-
- if (digit > 0)
- {
- // Now it's time to subtract our current quotient
- nuint carry = SubtractDivisor(left.Slice(n), right, digit);
- if (carry != t)
- {
- Debug.Assert(carry == t + 1);
-
- // Our guess was still exactly one too high
- carry = AddDivisor(left.Slice(n), right);
- --digit;
-
- Debug.Assert(carry == 1);
- }
- }
-
- // We have the digit!
- if ((uint)n < (uint)quotient.Length)
- {
- quotient[n] = digit;
- }
-
- if ((uint)i < (uint)left.Length)
- {
- left[i] = 0;
- }
- }
- }
-
- private static nuint AddDivisor(Span left, ReadOnlySpan right)
- {
- Debug.Assert(left.Length >= right.Length);
-
- // Repairs the dividend, if the last subtract was too much
-
- nuint carry = 0;
-
- for (int i = 0; i < right.Length; i++)
- {
- ref nuint leftElement = ref left[i];
- leftElement = AddWithCarry(leftElement, right[i], carry, out carry);
- }
-
- return carry;
- }
-
- private static nuint SubtractDivisor(Span left, ReadOnlySpan right, nuint q)
- {
- Debug.Assert(left.Length >= right.Length);
-
- return SubMul1(left, right, q);
- }
-
- private static bool DivideGuessTooBig(nuint q, nuint valHi1, nuint valHi0,
- nuint valLo, nuint divHi, nuint divLo)
- {
- // We multiply the two most significant limbs of the divisor
- // with the current guess for the quotient. If those are bigger
- // than the three most significant limbs of the current dividend
- // we return true, which means the current guess is still too big.
-
- nuint chkHiHi = nuint.BigMul(divHi, q, out nuint chkHiLo);
- nuint chkLoHi = nuint.BigMul(divLo, q, out nuint chkLoLo);
-
- chkHiLo += chkLoHi;
- if (chkHiLo < chkLoHi)
- {
- chkHiHi++;
- }
-
- return (chkHiHi > valHi1)
- || ((chkHiHi == valHi1) && ((chkHiLo > valHi0) || ((chkHiLo == valHi0) && (chkLoLo > valLo))));
- }
-
private static void DivideBurnikelZiegler(ReadOnlySpan left, ReadOnlySpan right, Span quotient, Span remainder)
{
Debug.Assert(left.Length >= 1);
diff --git a/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.SquMul.cs b/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.SquMul.cs
index ee501d657048ab..42e5bc80da1530 100644
--- a/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.SquMul.cs
+++ b/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.SquMul.cs
@@ -238,7 +238,7 @@ public static void Multiply(ReadOnlySpan left, ReadOnlySpan right,
if (right.Length < MultiplyKaratsubaThreshold)
{
- Naive(left, right, bits);
+ MultiplyNaive(left, right, bits);
}
else if ((left.Length + 1) >> 1 is int n && right.Length <= n)
{
@@ -514,24 +514,6 @@ static void RightSmall(ReadOnlySpan left, ReadOnlySpan right, Span
carryBuffer.Dispose();
}
-
- static void Naive(ReadOnlySpan left, ReadOnlySpan right, Span bits)
- {
- Debug.Assert(right.Length < MultiplyKaratsubaThreshold);
-
- // Multiplies the bits using the "grammar-school" method.
- // Envisioning the "rhombus" of a pen-and-paper calculation
- // should help getting the idea of these two loops...
- // The inner multiplication operations are safe, because
- // z_i+j + a_j * b_i + c <= 2(2^n - 1) + (2^n - 1)^2 =
- // = 2^(2n) - 1, where n = BitsPerLimb.
-
- for (int i = 0; i < right.Length; i++)
- {
- nuint carry = MulAdd1(bits.Slice(i), left, right[i]);
- bits[i + left.Length] = carry;
- }
- }
}
[StructLayout(LayoutKind.Auto)]
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 024d30f148ee50..1e9d4862e13a7f 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
@@ -2,84 +2,11 @@
// The .NET Foundation licenses this file to you under the MIT license.
using System.Diagnostics;
-using System.Runtime.CompilerServices;
-using System.Runtime.Intrinsics.X86;
-using X86Base = System.Runtime.Intrinsics.X86.X86Base;
namespace System.Numerics
{
internal static partial class BigIntegerCalculator
{
- /// Number of bits per native-width limb: 32 on 32-bit, 64 on 64-bit.
- internal static int BitsPerLimb
- {
- [MethodImpl(MethodImplOptions.AggressiveInlining)]
- get => nint.Size * 8;
- }
-
- public static int Compare(ReadOnlySpan left, ReadOnlySpan right)
- {
- Debug.Assert(left.Length <= right.Length || left.Slice(right.Length).ContainsAnyExcept(0u));
- Debug.Assert(left.Length >= right.Length || right.Slice(left.Length).ContainsAnyExcept(0u));
-
- if (left.Length != right.Length)
- {
- return left.Length < right.Length ? -1 : 1;
- }
-
- int iv = left.Length;
- while (--iv >= 0 && left[iv] == right[iv]) ;
-
- if (iv < 0)
- {
- return 0;
- }
-
- 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
- // of a given value may be less then the array's length
-
- int length = value.Length;
-
- while (length > 0 && value[length - 1] == 0)
- {
- --length;
- }
-
- return length;
- }
-
private static int Reduce(Span bits, ReadOnlySpan modulus)
{
// Executes a modulo operation using the divide operation.
@@ -100,298 +27,5 @@ public static void InitializeForDebug(Span bits)
// Reproduce the case where the return value of `stackalloc nuint` is not initialized to zero.
bits.Fill(0xCD);
}
-
- ///
- /// Performs widening addition of two limbs plus a carry-in, returning the sum and carry-out.
- /// On 64-bit: uses 128-bit arithmetic. On 32-bit: uses 64-bit arithmetic.
- ///
- [MethodImpl(MethodImplOptions.AggressiveInlining)]
- internal static nuint AddWithCarry(nuint a, nuint b, nuint carryIn, out nuint carryOut)
- {
- if (nint.Size == 8)
- {
- nuint sum1 = a + b;
- nuint c1 = (sum1 < a) ? 1 : (nuint)0;
- nuint sum2 = sum1 + carryIn;
- nuint c2 = (sum2 < sum1) ? 1 : (nuint)0;
- carryOut = c1 + c2;
- return sum2;
- }
- else
- {
- ulong sum = (ulong)a + b + carryIn;
- carryOut = (uint)(sum >> 32);
- return (uint)sum;
- }
- }
-
- ///
- /// Performs widening subtraction of two limbs with a borrow-in, returning the difference and borrow-out.
- /// borrowOut is 0 (no borrow) or 1 (borrow occurred).
- ///
- [MethodImpl(MethodImplOptions.AggressiveInlining)]
- internal static nuint SubWithBorrow(nuint a, nuint b, nuint borrowIn, out nuint borrowOut)
- {
- if (nint.Size == 8)
- {
- // Use unsigned underflow detection
- nuint diff1 = a - b;
- nuint b1 = (diff1 > a) ? 1 : (nuint)0;
- nuint diff2 = diff1 - borrowIn;
- nuint b2 = (diff2 > diff1) ? 1 : (nuint)0;
- borrowOut = b1 + b2;
- return diff2;
- }
- else
- {
- long diff = (long)a - (long)b - (long)borrowIn;
- borrowOut = (uint)(-(int)(diff >> 32)); // 0 or 1
- return (uint)diff;
- }
- }
-
- ///
- /// Widening divide: (hi:lo) / divisor -> (quotient, remainder).
- ///
- [MethodImpl(MethodImplOptions.AggressiveInlining)]
- internal static nuint DivRem(nuint hi, nuint lo, nuint divisor, out nuint remainder)
- {
- if (nint.Size == 8)
- {
- // Compute (hi * 2^64 + lo) / divisor.
- // hi < divisor is guaranteed by callers, so quotient fits in 64 bits.
- Debug.Assert(hi < (ulong)divisor || divisor == 0);
-
- if (hi == 0)
- {
- (ulong q, ulong r) = Math.DivRem(lo, (ulong)divisor);
- remainder = (nuint)r;
- return (nuint)q;
- }
-
- // When divisor fits in 32 bits, split lo into two 32-bit halves
- // and chain two native 64-bit divisions (avoids UInt128 overhead):
- // (hi * 2^32 + lo_hi) / divisor -> (q_hi, r1) [fits: hi < divisor < 2^32]
- // (r1 * 2^32 + lo_lo) / divisor -> (q_lo, r2) [fits: r1 < divisor < 2^32]
- if ((ulong)divisor <= uint.MaxValue)
- {
- ulong lo_hi = (ulong)lo >> 32;
- ulong lo_lo = (ulong)lo & 0xFFFFFFFF;
-
- (ulong q_hi, ulong r1) = Math.DivRem(((ulong)hi << 32) | lo_hi, divisor);
- (ulong q_lo, ulong r2) = Math.DivRem((r1 << 32) | lo_lo, divisor);
-
- remainder = (nuint)r2;
- return (nuint)((q_hi << 32) | q_lo);
- }
-
- {
-#pragma warning disable SYSLIB5004 // X86Base.DivRem is experimental
- if (X86Base.X64.IsSupported)
- {
- (ulong q, ulong r) = X86Base.X64.DivRem(lo, hi, divisor);
- remainder = (nuint)r;
- return (nuint)q;
- }
-#pragma warning restore SYSLIB5004
-
- UInt128 value = ((UInt128)(ulong)hi << 64) | (ulong)lo;
- UInt128 digit = value / (ulong)divisor;
- remainder = (nuint)(ulong)(value - digit * (ulong)divisor);
- return (nuint)(ulong)digit;
- }
- }
- else
- {
- ulong value = ((ulong)hi << 32) | lo;
- ulong digit = value / divisor;
- remainder = (uint)(value - digit * divisor);
- return (uint)digit;
- }
- }
-
- ///
- /// Multiply by scalar: result[0..left.Length] = left * multiplier.
- /// Returns the carry out. Unrolled by 4 on 64-bit.
- /// Unlike MulAdd1, this writes to result rather than accumulating.
- ///
- internal static nuint Mul1(Span result, ReadOnlySpan left, nuint multiplier)
- {
- Debug.Assert(result.Length >= left.Length);
-
- int length = left.Length;
- int i = 0;
- nuint carry = 0;
-
- if (nint.Size == 8)
- {
- for (; i + 3 < length; i += 4)
- {
- UInt128 p0 = (UInt128)(ulong)left[i] * (ulong)multiplier + (ulong)carry;
- result[i] = (nuint)(ulong)p0;
-
- UInt128 p1 = (UInt128)(ulong)left[i + 1] * (ulong)multiplier + (ulong)(p0 >> 64);
- result[i + 1] = (nuint)(ulong)p1;
-
- UInt128 p2 = (UInt128)(ulong)left[i + 2] * (ulong)multiplier + (ulong)(p1 >> 64);
- result[i + 2] = (nuint)(ulong)p2;
-
- UInt128 p3 = (UInt128)(ulong)left[i + 3] * (ulong)multiplier + (ulong)(p2 >> 64);
- result[i + 3] = (nuint)(ulong)p3;
-
- carry = (nuint)(ulong)(p3 >> 64);
- }
-
- for (; i < length; i++)
- {
- UInt128 product = (UInt128)(ulong)left[i] * (ulong)multiplier + (ulong)carry;
- result[i] = (nuint)(ulong)product;
- carry = (nuint)(ulong)(product >> 64);
- }
- }
- else
- {
- for (; i < length; i++)
- {
- ulong product = (ulong)left[i] * multiplier + carry;
- result[i] = (uint)product;
- carry = (uint)(product >> 32);
- }
- }
-
- return carry;
- }
-
- ///
- /// Fused multiply-accumulate by scalar: result[0..left.Length] += left * multiplier.
- /// Returns the carry out. Unrolled by 4 on 64-bit to overlap multiply latencies.
- ///
- internal static nuint MulAdd1(Span result, ReadOnlySpan left, nuint multiplier)
- {
- Debug.Assert(result.Length >= left.Length);
-
- int length = left.Length;
- int i = 0;
- nuint carry = 0;
-
- if (nint.Size == 8)
- {
- // Unroll by 4: mulx has 3-5 cycle latency but 1 cycle throughput,
- // so issuing 4 multiplies allows the CPU to pipeline them while
- // carry chains complete sequentially behind.
- for (; i + 3 < length; i += 4)
- {
- UInt128 p0 = (UInt128)(ulong)left[i] * (ulong)multiplier + (ulong)result[i] + (ulong)carry;
- result[i] = (nuint)(ulong)p0;
-
- UInt128 p1 = (UInt128)(ulong)left[i + 1] * (ulong)multiplier + (ulong)result[i + 1] + (ulong)(p0 >> 64);
- result[i + 1] = (nuint)(ulong)p1;
-
- UInt128 p2 = (UInt128)(ulong)left[i + 2] * (ulong)multiplier + (ulong)result[i + 2] + (ulong)(p1 >> 64);
- result[i + 2] = (nuint)(ulong)p2;
-
- UInt128 p3 = (UInt128)(ulong)left[i + 3] * (ulong)multiplier + (ulong)result[i + 3] + (ulong)(p2 >> 64);
- result[i + 3] = (nuint)(ulong)p3;
-
- carry = (nuint)(ulong)(p3 >> 64);
- }
-
- for (; i < length; i++)
- {
- UInt128 product = (UInt128)(ulong)left[i] * (ulong)multiplier + (ulong)result[i] + (ulong)carry;
- result[i] = (nuint)(ulong)product;
- carry = (nuint)(ulong)(product >> 64);
- }
- }
- else
- {
- for (; i < length; i++)
- {
- ulong product = (ulong)left[i] * multiplier
- + result[i] + carry;
- result[i] = (uint)product;
- carry = (uint)(product >> 32);
- }
- }
-
- return carry;
- }
-
- ///
- /// Fused subtract-multiply by scalar: result[0..right.Length] -= right * multiplier.
- /// Returns the borrow out. Unrolled by 4 on 64-bit.
- ///
- internal static nuint SubMul1(Span result, ReadOnlySpan right, nuint multiplier)
- {
- Debug.Assert(result.Length >= right.Length);
-
- int length = right.Length;
- int i = 0;
- nuint carry = 0;
-
- if (nint.Size == 8)
- {
- for (; i + 3 < length; i += 4)
- {
- UInt128 prod0 = (UInt128)(ulong)right[i] * (ulong)multiplier + (ulong)carry;
- nuint lo0 = (nuint)(ulong)prod0;
- nuint hi0 = (nuint)(ulong)(prod0 >> 64);
- nuint orig0 = result[i];
- result[i] = orig0 - lo0;
- hi0 += (orig0 < lo0) ? (nuint)1 : 0;
-
- UInt128 prod1 = (UInt128)(ulong)right[i + 1] * (ulong)multiplier + (ulong)hi0;
- nuint lo1 = (nuint)(ulong)prod1;
- nuint hi1 = (nuint)(ulong)(prod1 >> 64);
- nuint orig1 = result[i + 1];
- result[i + 1] = orig1 - lo1;
- hi1 += (orig1 < lo1) ? (nuint)1 : 0;
-
- UInt128 prod2 = (UInt128)(ulong)right[i + 2] * (ulong)multiplier + (ulong)hi1;
- nuint lo2 = (nuint)(ulong)prod2;
- nuint hi2 = (nuint)(ulong)(prod2 >> 64);
- nuint orig2 = result[i + 2];
- result[i + 2] = orig2 - lo2;
- hi2 += (orig2 < lo2) ? (nuint)1 : 0;
-
- UInt128 prod3 = (UInt128)(ulong)right[i + 3] * (ulong)multiplier + (ulong)hi2;
- nuint lo3 = (nuint)(ulong)prod3;
- nuint hi3 = (nuint)(ulong)(prod3 >> 64);
- nuint orig3 = result[i + 3];
- result[i + 3] = orig3 - lo3;
- hi3 += (orig3 < lo3) ? (nuint)1 : 0;
-
- carry = hi3;
- }
-
- for (; i < length; i++)
- {
- UInt128 product = (UInt128)(ulong)right[i] * (ulong)multiplier + (ulong)carry;
- nuint lo = (nuint)(ulong)product;
- nuint hi = (nuint)(ulong)(product >> 64);
- nuint orig = result[i];
- result[i] = orig - lo;
- hi += (orig < lo) ? (nuint)1 : 0;
- carry = hi;
- }
- }
- else
- {
- for (; i < length; i++)
- {
- ulong product = (ulong)right[i] * multiplier + carry;
- uint lo = (uint)product;
- uint hi = (uint)(product >> 32);
-
- uint orig = (uint)result[i];
- result[i] = orig - lo;
- hi += (orig < lo) ? 1u : 0;
-
- carry = hi;
- }
- }
-
- return carry;
- }
}
}