1 // Licensed to the .NET Foundation under one or more agreements.
2 // The .NET Foundation licenses this file to you under the MIT license.
3 // See the LICENSE file in the project root for more information.
5 using System.Diagnostics;
6 using System.Runtime.InteropServices;
7 using Internal.Runtime.CompilerServices;
11 internal static partial class Number
13 [StructLayout(LayoutKind.Sequential, Pack = 1)]
14 internal unsafe ref struct BigInteger
16 private const int MaxBlockCount = 35;
18 private static readonly uint[] s_Pow10UInt32Table = new uint[]
30 private static readonly int[] s_s_Pow10BigNumTableIndices = new int[]
40 private static readonly uint[] s_Pow10BigNumTable = new uint[]
48 0x6FC10000, // _blocks
53 0x00000000, // _blocks
60 0x00000000, // _blocks
70 0x00000000, // _blocks
87 0x00000000, // _blocks
115 // Trailing blocks to ensure MaxBlockCount
126 private static readonly uint[] s_MultiplyDeBruijnBitPosition = new uint[]
128 0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30,
129 8, 12, 20, 28, 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31
133 private BlocksBuffer _blocks;
135 public BigInteger(uint value)
138 _length = (value == 0) ? 0 : 1;
141 public BigInteger(ulong value)
143 var lower = (uint)(value);
144 var upper = (uint)(value >> 32);
149 _length = (upper == 0) ? 1 : 2;
152 public static uint BitScanReverse(uint mask)
154 // This comes from the Stanford Bit Widdling Hacks by Sean Eron Anderson:
155 // http://graphics.stanford.edu/~seander/bithacks.html#IntegerLogDeBruijn
157 mask |= (mask >> 1); // first round down to one less than a power of 2
161 mask |= (mask >> 16);
163 uint index = (mask * 0x07C4ACDD) >> 27;
164 return s_MultiplyDeBruijnBitPosition[(int)(index)];
167 public static int Compare(ref BigInteger lhs, ref BigInteger rhs)
169 Debug.Assert(unchecked((uint)(lhs._length)) <= MaxBlockCount);
170 Debug.Assert(unchecked((uint)(rhs._length)) <= MaxBlockCount);
172 int lhsLength = lhs._length;
173 int rhsLength = rhs._length;
175 int lengthDelta = (lhsLength - rhsLength);
177 if (lengthDelta != 0)
184 Debug.Assert(rhsLength == 0);
188 for (int index = (lhsLength - 1); index >= 0; index--)
190 long delta = (long)(lhs._blocks[index]) - rhs._blocks[index];
194 return delta > 0 ? 1 : -1;
201 public static uint HeuristicDivide(ref BigInteger dividend, ref BigInteger divisor)
203 int divisorLength = divisor._length;
205 if (dividend._length < divisorLength)
210 // This is an estimated quotient. Its error should be less than 2.
211 // Reference inequality:
212 // a/b - floor(floor(a)/(floor(b) + 1)) < 2
213 int lastIndex = (divisorLength - 1);
214 uint quotient = dividend._blocks[lastIndex] / (divisor._blocks[lastIndex] + 1);
218 // Now we use our estimated quotient to update each block of dividend.
219 // dividend = dividend - divisor * quotient
227 ulong product = ((ulong)(divisor._blocks[index]) * quotient) + carry;
228 carry = product >> 32;
230 ulong difference = (ulong)(dividend._blocks[index]) - (uint)(product) - borrow;
231 borrow = (difference >> 32) & 1;
233 dividend._blocks[index] = (uint)(difference);
237 while (index < divisorLength);
239 // Remove all leading zero blocks from dividend
240 while ((divisorLength > 0) && (dividend._blocks[divisorLength - 1] == 0))
245 dividend._length = divisorLength;
248 // If the dividend is still larger than the divisor, we overshot our estimate quotient. To correct,
249 // we increment the quotient and subtract one more divisor from the dividend (Because we guaranteed the error range).
250 if (Compare(ref dividend, ref divisor) >= 0)
254 // dividend = dividend - divisor
260 ulong difference = (ulong)(dividend._blocks[index]) - divisor._blocks[index] - borrow;
261 borrow = (difference >> 32) & 1;
263 dividend._blocks[index] = (uint)(difference);
267 while (index < divisorLength);
269 // Remove all leading zero blocks from dividend
270 while ((divisorLength > 0) && (dividend._blocks[divisorLength - 1] == 0))
275 dividend._length = divisorLength;
281 public static uint LogBase2(uint value)
283 Debug.Assert(value != 0);
284 return BitScanReverse(value);
287 public static uint LogBase2(ulong value)
289 Debug.Assert(value != 0);
291 uint upper = (uint)(value >> 32);
295 return 32 + LogBase2(upper);
298 return LogBase2((uint)(value));
301 public static void Multiply(ref BigInteger lhs, uint value, ref BigInteger result)
303 if (lhs.IsZero() || (value == 1))
305 result.SetValue(ref lhs);
315 int lhsLength = lhs._length;
320 while (index < lhsLength)
322 ulong product = ((ulong)(lhs._blocks[index]) * value) + carry;
323 carry = product >> 32;
324 result._blocks[index] = (uint)(product);
331 Debug.Assert(unchecked((uint)(lhsLength)) + 1 <= MaxBlockCount);
332 result._blocks[index] = (uint)(carry);
333 result._length += (lhsLength + 1);
337 public static void Multiply(ref BigInteger lhs, ref BigInteger rhs, ref BigInteger result)
339 if (lhs.IsZero() || rhs.IsOne())
341 result.SetValue(ref lhs);
351 ref readonly BigInteger large = ref lhs;
352 int largeLength = lhs._length;
354 ref readonly BigInteger small = ref rhs;
355 int smallLength = rhs._length;
357 if (largeLength < smallLength)
360 largeLength = rhs._length;
363 smallLength = lhs._length;
366 int maxResultLength = smallLength + largeLength;
367 Debug.Assert(unchecked((uint)(maxResultLength)) <= MaxBlockCount);
369 // Zero out result internal blocks.
370 Buffer.ZeroMemory((byte*)(result._blocks.GetPointer()), (maxResultLength * sizeof(uint)));
373 int resultStartIndex = 0;
375 while (smallIndex < smallLength)
377 // Multiply each block of large BigNum.
378 if (small._blocks[smallIndex] != 0)
381 int resultIndex = resultStartIndex;
387 ulong product = result._blocks[resultIndex] + ((ulong)(small._blocks[smallIndex]) * large._blocks[largeIndex]) + carry;
388 carry = product >> 32;
389 result._blocks[resultIndex] = (uint)(product);
394 while (largeIndex < largeLength);
396 result._blocks[resultIndex] = (uint)(carry);
403 if ((maxResultLength > 0) && (result._blocks[maxResultLength - 1] == 0))
405 result._length = (maxResultLength - 1);
409 result._length = maxResultLength;
413 public static void Pow10(uint exponent, ref BigInteger result)
415 // We leverage two arrays - s_Pow10UInt32Table and s_Pow10BigNumTable to speed up the Pow10 calculation.
417 // s_Pow10UInt32Table stores the results of 10^0 to 10^7.
418 // s_Pow10BigNumTable stores the results of 10^8, 10^16, 10^32, 10^64, 10^128 and 10^256.
420 // For example, let's say exp = 0b111111. We can split the exp to two parts, one is small exp,
421 // which 10^smallExp can be represented as uint, another part is 10^bigExp, which must be represented as BigNum.
422 // So the result should be 10^smallExp * 10^bigExp.
424 // Calculating 10^smallExp is simple, we just lookup the 10^smallExp from s_Pow10UInt32Table.
425 // But here's a bad news: although uint can represent 10^9, exp 9's binary representation is 1001.
426 // That means 10^(1011), 10^(1101), 10^(1111) all cannot be stored as uint, we cannot easily say something like:
427 // "Any bits <= 3 is small exp, any bits > 3 is big exp". So instead of involving 10^8, 10^9 to s_Pow10UInt32Table,
428 // consider 10^8 and 10^9 as a bigNum, so they fall into s_Pow10BigNumTable. Now we can have a simple rule:
429 // "Any bits <= 3 is small exp, any bits > 3 is big exp".
431 // For 0b111111, we first calculate 10^(smallExp), which is 10^(7), now we can shift right 3 bits, prepare to calculate the bigExp part,
432 // the exp now becomes 0b000111.
434 // Apparently the lowest bit of bigExp should represent 10^8 because we have already shifted 3 bits for smallExp, so s_Pow10BigNumTable[0] = 10^8.
435 // Now let's shift exp right 1 bit, the lowest bit should represent 10^(8 * 2) = 10^16, and so on...
437 // That's why we just need the values of s_Pow10BigNumTable be power of 2.
439 // More details of this implementation can be found at: https://github.com/dotnet/coreclr/pull/12894#discussion_r128890596
441 BigInteger temp1 = new BigInteger(s_Pow10UInt32Table[exponent & 0x7]);
442 ref BigInteger lhs = ref temp1;
444 BigInteger temp2 = new BigInteger(0);
445 ref BigInteger product = ref temp2;
450 while (exponent != 0)
452 // If the current bit is set, multiply it with the corresponding power of 10
453 if ((exponent & 1) != 0)
455 // Multiply into the next temporary
456 ref BigInteger rhs = ref *(BigInteger*)(Unsafe.AsPointer(ref s_Pow10BigNumTable[s_s_Pow10BigNumTableIndices[index]]));
457 Multiply(ref lhs, ref rhs, ref product);
459 // Swap to the next temporary
460 ref BigInteger temp = ref product;
465 // Advance to the next bit
470 result.SetValue(ref lhs);
473 public static void PrepareHeuristicDivide(ref BigInteger dividend, ref BigInteger divisor)
475 uint hiBlock = divisor._blocks[divisor._length - 1];
477 if ((hiBlock < 8) || (hiBlock > 429496729))
479 // Inspired by http://www.ryanjuckett.com/programming/printing-floating-point-numbers/
480 // Perform a bit shift on all values to get the highest block of the divisor into
481 // the range [8,429496729]. We are more likely to make accurate quotient estimations
482 // in heuristicDivide() with higher divisor values so
483 // we shift the divisor to place the highest bit at index 27 of the highest block.
484 // This is safe because (2^28 - 1) = 268435455 which is less than 429496729. This means
485 // that all values with a highest bit at index 27 are within range.
486 uint hiBlockLog2 = LogBase2(hiBlock);
487 uint shift = (59 - hiBlockLog2) % 32;
489 divisor.ShiftLeft(shift);
490 dividend.ShiftLeft(shift);
494 public static void ShiftLeft(ulong input, uint shift, ref BigInteger output)
501 uint blocksToShift = Math.DivRem(shift, 32, out uint remainingBitsToShift);
503 if (blocksToShift > 0)
505 // If blocks shifted, we should fill the corresponding blocks with zero.
506 output.ExtendBlocks(0, blocksToShift);
509 if (remainingBitsToShift == 0)
511 // We shift 32 * n (n >= 1) bits. No remaining bits.
512 output.ExtendBlock((uint)(input));
514 uint highBits = (uint)(input >> 32);
518 output.ExtendBlock(highBits);
523 // Extract the high position bits which would be shifted out of range.
524 uint highPositionBits = (uint)(input) >> (int)(64 - remainingBitsToShift);
526 // Shift the input. The result should be stored to current block.
527 ulong shiftedInput = input << (int)(remainingBitsToShift);
528 output.ExtendBlock((uint)(shiftedInput));
530 uint highBits = (uint)(input >> 32);
534 output.ExtendBlock(highBits);
537 if (highPositionBits != 0)
539 // If the high position bits is not 0, we should store them to next block.
540 output.ExtendBlock(highPositionBits);
545 public void ExtendBlock(uint blockValue)
547 _blocks[_length] = blockValue;
551 public void ExtendBlocks(uint blockValue, uint blockCount)
553 Debug.Assert(blockCount > 0);
557 ExtendBlock(blockValue);
562 Buffer.ZeroMemory((byte*)(_blocks.GetPointer() + _length), ((blockCount - 1) * sizeof(uint)));
563 _length += (int)(blockCount);
564 _blocks[_length - 1] = blockValue;
569 return (_length == 1)
570 && (_blocks[0] == 1);
578 public void Multiply(uint value)
580 Multiply(ref this, value, ref this);
583 public void Multiply(ref BigInteger value)
585 var result = new BigInteger(0);
586 Multiply(ref this, ref value, ref result);
588 Buffer.Memcpy((byte*)(_blocks.GetPointer()), (byte*)(result._blocks.GetPointer()), (result._length) * sizeof(uint));
589 _length = result._length;
592 public void Multiply10()
600 int length = _length;
603 while (index < length)
605 var block = (ulong)(_blocks[index]);
606 ulong product = (block << 3) + (block << 1) + carry;
607 carry = product >> 32;
608 _blocks[index] = (uint)(product);
615 Debug.Assert(unchecked((uint)(_length)) + 1 <= MaxBlockCount);
616 _blocks[index] = (uint)(carry);
621 public void SetUInt32(uint value)
627 public void SetUInt64(ulong value)
629 var lower = (uint)(value);
630 var upper = (uint)(value >> 32);
635 _length = (upper == 0) ? 1 : 2;
638 public void SetValue(ref BigInteger rhs)
640 int rhsLength = rhs._length;
641 Buffer.Memcpy((byte*)(_blocks.GetPointer()), (byte*)(rhs._blocks.GetPointer()), (rhsLength * sizeof(uint)));
645 public void SetZero()
650 public void ShiftLeft(uint shift)
652 // Process blocks high to low so that we can safely process in place
653 var length = _length;
655 if ((length == 0) || (shift == 0))
660 uint blocksToShift = Math.DivRem(shift, 32, out uint remainingBitsToShift);
662 // Copy blocks from high to low
663 int readIndex = (length - 1);
664 int writeIndex = readIndex + (int)(blocksToShift);
666 // Check if the shift is block aligned
667 if (remainingBitsToShift == 0)
669 Debug.Assert(writeIndex < MaxBlockCount);
671 while (readIndex >= 0)
673 _blocks[writeIndex] = _blocks[readIndex];
678 _length += (int)(blocksToShift);
680 // Zero the remaining low blocks
681 Buffer.ZeroMemory((byte*)(_blocks.GetPointer()), (blocksToShift * sizeof(uint)));
685 // We need an extra block for the partial shift
687 Debug.Assert(writeIndex < MaxBlockCount);
689 // Set the length to hold the shifted blocks
690 _length = writeIndex + 1;
692 // Output the initial blocks
693 uint lowBitsShift = (32 - remainingBitsToShift);
695 uint block = _blocks[readIndex];
696 uint lowBits = block >> (int)(lowBitsShift);
697 while (readIndex > 0)
699 _blocks[writeIndex] = highBits | lowBits;
700 highBits = block << (int)(remainingBitsToShift);
705 block = _blocks[readIndex];
706 lowBits = block >> (int)lowBitsShift;
709 // Output the final blocks
710 _blocks[writeIndex] = highBits | lowBits;
711 _blocks[writeIndex - 1] = block << (int)(remainingBitsToShift);
713 // Zero the remaining low blocks
714 Buffer.ZeroMemory((byte*)(_blocks.GetPointer()), (blocksToShift * sizeof(uint)));
716 // Check if the terminating block has no set bits
717 if (_blocks[_length - 1] == 0)
724 [StructLayout(LayoutKind.Sequential, Pack = 1)]
725 private struct BlocksBuffer
727 private fixed uint _blocks[MaxBlockCount];
729 public ref uint this[int index]
733 Debug.Assert(unchecked((uint)(index)) <= MaxBlockCount);
734 return ref Unsafe.Add(ref GetPinnableReference(), index);
738 public ref uint GetPinnableReference()
740 var pThis = Unsafe.AsPointer(ref this);
741 return ref Unsafe.AsRef<uint>(pThis);
744 public uint* GetPointer()
746 return (uint*)(Unsafe.AsPointer(ref this));