Number.BigInteger.cs 46 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350
  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.
  4. using System.Diagnostics;
  5. using System.Runtime.InteropServices;
  6. using Internal.Runtime.CompilerServices;
  7. namespace System
  8. {
  9. internal static partial class Number
  10. {
  11. [StructLayout(LayoutKind.Sequential, Pack = 1)]
  12. internal unsafe ref struct BigInteger
  13. {
  14. // The longest binary mantissa requires: explicit mantissa bits + abs(min exponent)
  15. // * Half: 10 + 14 = 24
  16. // * Single: 23 + 126 = 149
  17. // * Double: 52 + 1022 = 1074
  18. // * Quad: 112 + 16382 = 16494
  19. private const int BitsForLongestBinaryMantissa = 1074;
  20. // The longest digit sequence requires: ceil(log2(pow(10, max significant digits + 1 rounding digit)))
  21. // * Half: ceil(log2(pow(10, 21 + 1))) = 74
  22. // * Single: ceil(log2(pow(10, 112 + 1))) = 376
  23. // * Double: ceil(log2(pow(10, 767 + 1))) = 2552
  24. // * Quad: ceil(log2(pow(10, 11563 + 1))) = 38415
  25. private const int BitsForLongestDigitSequence = 2552;
  26. // We require BitsPerBlock additional bits for shift space used during the pre-division preparation
  27. private const int MaxBits = BitsForLongestBinaryMantissa + BitsForLongestDigitSequence + BitsPerBlock;
  28. private const int BitsPerBlock = sizeof(int) * 8;
  29. private const int MaxBlockCount = (MaxBits + (BitsPerBlock - 1)) / BitsPerBlock;
  30. private static readonly uint[] s_Pow10UInt32Table = new uint[]
  31. {
  32. 1, // 10^0
  33. 10, // 10^1
  34. 100, // 10^2
  35. 1000, // 10^3
  36. 10000, // 10^4
  37. 100000, // 10^5
  38. 1000000, // 10^6
  39. 10000000, // 10^7
  40. };
  41. private static readonly int[] s_Pow10BigNumTableIndices = new int[]
  42. {
  43. 0, // 10^8
  44. 2, // 10^16
  45. 5, // 10^32
  46. 10, // 10^64
  47. 18, // 10^128
  48. 33, // 10^256
  49. 61, // 10^512
  50. 116, // 10^1024
  51. };
  52. private static readonly uint[] s_Pow10BigNumTable = new uint[]
  53. {
  54. // 10^8
  55. 1, // _length
  56. 100000000, // _blocks
  57. // 10^16
  58. 2, // _length
  59. 0x6FC10000, // _blocks
  60. 0x002386F2,
  61. // 10^32
  62. 4, // _length
  63. 0x00000000, // _blocks
  64. 0x85ACEF81,
  65. 0x2D6D415B,
  66. 0x000004EE,
  67. // 10^64
  68. 7, // _length
  69. 0x00000000, // _blocks
  70. 0x00000000,
  71. 0xBF6A1F01,
  72. 0x6E38ED64,
  73. 0xDAA797ED,
  74. 0xE93FF9F4,
  75. 0x00184F03,
  76. // 10^128
  77. 14, // _length
  78. 0x00000000, // _blocks
  79. 0x00000000,
  80. 0x00000000,
  81. 0x00000000,
  82. 0x2E953E01,
  83. 0x03DF9909,
  84. 0x0F1538FD,
  85. 0x2374E42F,
  86. 0xD3CFF5EC,
  87. 0xC404DC08,
  88. 0xBCCDB0DA,
  89. 0xA6337F19,
  90. 0xE91F2603,
  91. 0x0000024E,
  92. // 10^256
  93. 27, // _length
  94. 0x00000000, // _blocks
  95. 0x00000000,
  96. 0x00000000,
  97. 0x00000000,
  98. 0x00000000,
  99. 0x00000000,
  100. 0x00000000,
  101. 0x00000000,
  102. 0x982E7C01,
  103. 0xBED3875B,
  104. 0xD8D99F72,
  105. 0x12152F87,
  106. 0x6BDE50C6,
  107. 0xCF4A6E70,
  108. 0xD595D80F,
  109. 0x26B2716E,
  110. 0xADC666B0,
  111. 0x1D153624,
  112. 0x3C42D35A,
  113. 0x63FF540E,
  114. 0xCC5573C0,
  115. 0x65F9EF17,
  116. 0x55BC28F2,
  117. 0x80DCC7F7,
  118. 0xF46EEDDC,
  119. 0x5FDCEFCE,
  120. 0x000553F7,
  121. // 10^512
  122. 54, // _length
  123. 0x00000000, // _blocks
  124. 0x00000000,
  125. 0x00000000,
  126. 0x00000000,
  127. 0x00000000,
  128. 0x00000000,
  129. 0x00000000,
  130. 0x00000000,
  131. 0x00000000,
  132. 0x00000000,
  133. 0x00000000,
  134. 0x00000000,
  135. 0x00000000,
  136. 0x00000000,
  137. 0x00000000,
  138. 0x00000000,
  139. 0xFC6CF801,
  140. 0x77F27267,
  141. 0x8F9546DC,
  142. 0x5D96976F,
  143. 0xB83A8A97,
  144. 0xC31E1AD9,
  145. 0x46C40513,
  146. 0x94E65747,
  147. 0xC88976C1,
  148. 0x4475B579,
  149. 0x28F8733B,
  150. 0xAA1DA1BF,
  151. 0x703ED321,
  152. 0x1E25CFEA,
  153. 0xB21A2F22,
  154. 0xBC51FB2E,
  155. 0x96E14F5D,
  156. 0xBFA3EDAC,
  157. 0x329C57AE,
  158. 0xE7FC7153,
  159. 0xC3FC0695,
  160. 0x85A91924,
  161. 0xF95F635E,
  162. 0xB2908EE0,
  163. 0x93ABADE4,
  164. 0x1366732A,
  165. 0x9449775C,
  166. 0x69BE5B0E,
  167. 0x7343AFAC,
  168. 0xB099BC81,
  169. 0x45A71D46,
  170. 0xA2699748,
  171. 0x8CB07303,
  172. 0x8A0B1F13,
  173. 0x8CAB8A97,
  174. 0xC1D238D9,
  175. 0x633415D4,
  176. 0x0000001C,
  177. // 10^1024
  178. 107, // _length
  179. 0x00000000, // _blocks
  180. 0x00000000,
  181. 0x00000000,
  182. 0x00000000,
  183. 0x00000000,
  184. 0x00000000,
  185. 0x00000000,
  186. 0x00000000,
  187. 0x00000000,
  188. 0x00000000,
  189. 0x00000000,
  190. 0x00000000,
  191. 0x00000000,
  192. 0x00000000,
  193. 0x00000000,
  194. 0x00000000,
  195. 0x00000000,
  196. 0x00000000,
  197. 0x00000000,
  198. 0x00000000,
  199. 0x00000000,
  200. 0x00000000,
  201. 0x00000000,
  202. 0x00000000,
  203. 0x00000000,
  204. 0x00000000,
  205. 0x00000000,
  206. 0x00000000,
  207. 0x00000000,
  208. 0x00000000,
  209. 0x00000000,
  210. 0x00000000,
  211. 0x2919F001,
  212. 0xF55B2B72,
  213. 0x6E7C215B,
  214. 0x1EC29F86,
  215. 0x991C4E87,
  216. 0x15C51A88,
  217. 0x140AC535,
  218. 0x4C7D1E1A,
  219. 0xCC2CD819,
  220. 0x0ED1440E,
  221. 0x896634EE,
  222. 0x7DE16CFB,
  223. 0x1E43F61F,
  224. 0x9FCE837D,
  225. 0x231D2B9C,
  226. 0x233E55C7,
  227. 0x65DC60D7,
  228. 0xF451218B,
  229. 0x1C5CD134,
  230. 0xC9635986,
  231. 0x922BBB9F,
  232. 0xA7E89431,
  233. 0x9F9F2A07,
  234. 0x62BE695A,
  235. 0x8E1042C4,
  236. 0x045B7A74,
  237. 0x1ABE1DE3,
  238. 0x8AD822A5,
  239. 0xBA34C411,
  240. 0xD814B505,
  241. 0xBF3FDEB3,
  242. 0x8FC51A16,
  243. 0xB1B896BC,
  244. 0xF56DEEEC,
  245. 0x31FB6BFD,
  246. 0xB6F4654B,
  247. 0x101A3616,
  248. 0x6B7595FB,
  249. 0xDC1A47FE,
  250. 0x80D98089,
  251. 0x80BDA5A5,
  252. 0x9A202882,
  253. 0x31EB0F66,
  254. 0xFC8F1F90,
  255. 0x976A3310,
  256. 0xE26A7B7E,
  257. 0xDF68368A,
  258. 0x3CE3A0B8,
  259. 0x8E4262CE,
  260. 0x75A351A2,
  261. 0x6CB0B6C9,
  262. 0x44597583,
  263. 0x31B5653F,
  264. 0xC356E38A,
  265. 0x35FAABA6,
  266. 0x0190FBA0,
  267. 0x9FC4ED52,
  268. 0x88BC491B,
  269. 0x1640114A,
  270. 0x005B8041,
  271. 0xF4F3235E,
  272. 0x1E8D4649,
  273. 0x36A8DE06,
  274. 0x73C55349,
  275. 0xA7E6BD2A,
  276. 0xC1A6970C,
  277. 0x47187094,
  278. 0xD2DB49EF,
  279. 0x926C3F5B,
  280. 0xAE6209D4,
  281. 0x2D433949,
  282. 0x34F4A3C6,
  283. 0xD4305D94,
  284. 0xD9D61A05,
  285. 0x00000325,
  286. // 9 Trailing blocks to ensure MaxBlockCount
  287. 0x00000000,
  288. 0x00000000,
  289. 0x00000000,
  290. 0x00000000,
  291. 0x00000000,
  292. 0x00000000,
  293. 0x00000000,
  294. 0x00000000,
  295. 0x00000000,
  296. };
  297. private static readonly uint[] s_MultiplyDeBruijnBitPosition = new uint[]
  298. {
  299. 0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30,
  300. 8, 12, 20, 28, 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31
  301. };
  302. private int _length;
  303. private fixed uint _blocks[MaxBlockCount];
  304. public BigInteger(uint value)
  305. {
  306. _blocks[0] = value;
  307. _length = (value == 0) ? 0 : 1;
  308. }
  309. public BigInteger(ulong value)
  310. {
  311. var lower = (uint)(value);
  312. var upper = (uint)(value >> 32);
  313. _blocks[0] = lower;
  314. _blocks[1] = upper;
  315. _length = (upper == 0) ? 1 : 2;
  316. }
  317. public static void Add(ref BigInteger lhs, uint value, ref BigInteger result)
  318. {
  319. if (lhs.IsZero())
  320. {
  321. result.SetUInt32(value);
  322. return;
  323. }
  324. if (value == 0)
  325. {
  326. result.SetValue(ref lhs);
  327. return;
  328. }
  329. int lhsLength = lhs._length;
  330. int index = 0;
  331. uint carry = value;
  332. while (index < lhsLength)
  333. {
  334. ulong sum = (ulong)(lhs._blocks[index]) + carry;
  335. lhs._blocks[index] = (uint)(sum);
  336. carry = (uint)(sum >> 32);
  337. index++;
  338. }
  339. if (carry != 0)
  340. {
  341. Debug.Assert(unchecked((uint)(lhsLength)) + 1 <= MaxBlockCount);
  342. result._blocks[index] = carry;
  343. result._length = (lhsLength + 1);
  344. }
  345. }
  346. public static void Add(ref BigInteger lhs, ref BigInteger rhs, out BigInteger result)
  347. {
  348. // determine which operand has the smaller length
  349. ref BigInteger large = ref (lhs._length < rhs._length) ? ref rhs : ref lhs;
  350. ref BigInteger small = ref (lhs._length < rhs._length) ? ref lhs : ref rhs;
  351. int largeLength = large._length;
  352. int smallLength = small._length;
  353. // The output will be at least as long as the largest input
  354. result = new BigInteger(0);
  355. result._length = largeLength;
  356. // Add each block and add carry the overflow to the next block
  357. ulong carry = 0;
  358. int largeIndex = 0;
  359. int smallIndex = 0;
  360. int resultIndex = 0;
  361. while (smallIndex < smallLength)
  362. {
  363. ulong sum = carry + large._blocks[largeIndex] + small._blocks[smallIndex];
  364. carry = sum >> 32;
  365. result._blocks[resultIndex] = (uint)(sum);
  366. largeIndex++;
  367. smallIndex++;
  368. resultIndex++;
  369. }
  370. // Add the carry to any blocks that only exist in the large operand
  371. while (largeIndex < largeLength)
  372. {
  373. ulong sum = carry + large._blocks[largeIndex];
  374. carry = sum >> 32;
  375. result._blocks[resultIndex] = (uint)(sum);
  376. largeIndex++;
  377. resultIndex++;
  378. }
  379. // If there's still a carry, append a new block
  380. if (carry != 0)
  381. {
  382. Debug.Assert(carry == 1);
  383. Debug.Assert((resultIndex == largeLength) && (largeLength < MaxBlockCount));
  384. result._blocks[resultIndex] = 1;
  385. result._length += 1;
  386. }
  387. }
  388. public static int Compare(ref BigInteger lhs, ref BigInteger rhs)
  389. {
  390. Debug.Assert(unchecked((uint)(lhs._length)) <= MaxBlockCount);
  391. Debug.Assert(unchecked((uint)(rhs._length)) <= MaxBlockCount);
  392. int lhsLength = lhs._length;
  393. int rhsLength = rhs._length;
  394. int lengthDelta = (lhsLength - rhsLength);
  395. if (lengthDelta != 0)
  396. {
  397. return lengthDelta;
  398. }
  399. if (lhsLength == 0)
  400. {
  401. Debug.Assert(rhsLength == 0);
  402. return 0;
  403. }
  404. for (int index = (lhsLength - 1); index >= 0; index--)
  405. {
  406. long delta = (long)(lhs._blocks[index]) - rhs._blocks[index];
  407. if (delta != 0)
  408. {
  409. return delta > 0 ? 1 : -1;
  410. }
  411. }
  412. return 0;
  413. }
  414. public static uint CountSignificantBits(uint value)
  415. {
  416. return (value != 0) ? (1 + LogBase2(value)) : 0;
  417. }
  418. public static uint CountSignificantBits(ulong value)
  419. {
  420. uint upper = (uint)(value >> 32);
  421. if (upper != 0)
  422. {
  423. return 32 + CountSignificantBits(upper);
  424. }
  425. return CountSignificantBits((uint)(value));
  426. }
  427. public static uint CountSignificantBits(ref BigInteger value)
  428. {
  429. if (value.IsZero())
  430. {
  431. return 0;
  432. }
  433. // We don't track any unused blocks, so we only need to do a BSR on the
  434. // last index and add that to the number of bits we skipped.
  435. uint lastIndex = (uint)(value._length - 1);
  436. return (lastIndex * BitsPerBlock) + CountSignificantBits(value._blocks[lastIndex]);
  437. }
  438. public static void DivRem(ref BigInteger lhs, ref BigInteger rhs, out BigInteger quo, out BigInteger rem)
  439. {
  440. // This is modified from the CoreFX BigIntegerCalculator.DivRem.cs implementation:
  441. // https://github.com/dotnet/corefx/blob/0bb106232745aedfc0d0c5a84ff2b244bf190317/src/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.DivRem.cs
  442. Debug.Assert(!rhs.IsZero());
  443. quo = new BigInteger(0);
  444. rem = new BigInteger(0);
  445. if (lhs.IsZero())
  446. {
  447. return;
  448. }
  449. int lhsLength = lhs._length;
  450. int rhsLength = rhs._length;
  451. if ((lhsLength == 1) && (rhsLength == 1))
  452. {
  453. uint quotient = Math.DivRem(lhs._blocks[0], rhs._blocks[0], out uint remainder);
  454. quo = new BigInteger(quotient);
  455. rem = new BigInteger(remainder);
  456. return;
  457. }
  458. if (rhsLength == 1)
  459. {
  460. // We can make the computation much simpler if the rhs is only one block
  461. int quoLength = lhsLength;
  462. ulong rhsValue = rhs._blocks[0];
  463. ulong carry = 0;
  464. for (int i = quoLength - 1; i >= 0; i--)
  465. {
  466. ulong value = (carry << 32) | lhs._blocks[i];
  467. ulong digit = Math.DivRem(value, rhsValue, out carry);
  468. if ((digit == 0) && (i == (quoLength - 1)))
  469. {
  470. quoLength--;
  471. }
  472. else
  473. {
  474. quo._blocks[i] = (uint)(digit);
  475. }
  476. }
  477. quo._length = quoLength;
  478. rem.SetUInt32((uint)(carry));
  479. return;
  480. }
  481. else if (rhsLength > lhsLength)
  482. {
  483. // Handle the case where we have no quotient
  484. rem.SetValue(ref lhs);
  485. return;
  486. }
  487. else
  488. {
  489. int quoLength = lhsLength - rhsLength + 1;
  490. rem.SetValue(ref lhs);
  491. int remLength = lhsLength;
  492. // Executes the "grammar-school" algorithm for computing q = a / b.
  493. // Before calculating q_i, we get more bits into the highest bit
  494. // block of the divisor. Thus, guessing digits of the quotient
  495. // will be more precise. Additionally we'll get r = a % b.
  496. uint divHi = rhs._blocks[rhsLength - 1];
  497. uint divLo = rhs._blocks[rhsLength - 2];
  498. // We measure the leading zeros of the divisor
  499. int shiftLeft = (int)(LeadingZeroCount(divHi));
  500. int shiftRight = 32 - shiftLeft;
  501. // And, we make sure the most significant bit is set
  502. if (shiftLeft > 0)
  503. {
  504. divHi = (divHi << shiftLeft) | (divLo >> shiftRight);
  505. divLo = (divLo << shiftLeft);
  506. if (rhsLength > 2)
  507. {
  508. divLo |= (rhs._blocks[rhsLength - 3] >> shiftRight);
  509. }
  510. }
  511. // Then, we divide all of the bits as we would do it using
  512. // pen and paper: guessing the next digit, subtracting, ...
  513. for (int i = lhsLength; i >= rhsLength; i--)
  514. {
  515. int n = i - rhsLength;
  516. uint t = i < lhsLength ? rem._blocks[i] : 0;
  517. ulong valHi = ((ulong)(t) << 32) | rem._blocks[i - 1];
  518. uint valLo = i > 1 ? rem._blocks[i - 2] : 0;
  519. // We shifted the divisor, we shift the dividend too
  520. if (shiftLeft > 0)
  521. {
  522. valHi = (valHi << shiftLeft) | (valLo >> shiftRight);
  523. valLo = (valLo << shiftLeft);
  524. if (i > 2)
  525. {
  526. valLo |= (rem._blocks[i - 3] >> shiftRight);
  527. }
  528. }
  529. // First guess for the current digit of the quotient,
  530. // which naturally must have only 32 bits...
  531. ulong digit = valHi / divHi;
  532. if (digit > uint.MaxValue)
  533. {
  534. digit = uint.MaxValue;
  535. }
  536. // Our first guess may be a little bit to big
  537. while (DivideGuessTooBig(digit, valHi, valLo, divHi, divLo))
  538. {
  539. digit--;
  540. }
  541. if (digit > 0)
  542. {
  543. // Now it's time to subtract our current quotient
  544. uint carry = SubtractDivisor(ref rem, n, ref rhs, digit);
  545. if (carry != t)
  546. {
  547. Debug.Assert(carry == t + 1);
  548. // Our guess was still exactly one too high
  549. carry = AddDivisor(ref rem, n, ref rhs);
  550. digit--;
  551. Debug.Assert(carry == 1);
  552. }
  553. }
  554. // We have the digit!
  555. if (quoLength != 0)
  556. {
  557. if ((digit == 0) && (n == (quoLength - 1)))
  558. {
  559. quoLength--;
  560. }
  561. else
  562. {
  563. quo._blocks[n] = (uint)(digit);
  564. }
  565. }
  566. if (i < remLength)
  567. {
  568. remLength--;
  569. }
  570. }
  571. quo._length = quoLength;
  572. // We need to check for the case where remainder is zero
  573. for (int i = remLength - 1; i >= 0; i--)
  574. {
  575. if (rem._blocks[i] == 0)
  576. {
  577. remLength--;
  578. }
  579. }
  580. rem._length = remLength;
  581. }
  582. }
  583. public static uint HeuristicDivide(ref BigInteger dividend, ref BigInteger divisor)
  584. {
  585. int divisorLength = divisor._length;
  586. if (dividend._length < divisorLength)
  587. {
  588. return 0;
  589. }
  590. // This is an estimated quotient. Its error should be less than 2.
  591. // Reference inequality:
  592. // a/b - floor(floor(a)/(floor(b) + 1)) < 2
  593. int lastIndex = (divisorLength - 1);
  594. uint quotient = dividend._blocks[lastIndex] / (divisor._blocks[lastIndex] + 1);
  595. if (quotient != 0)
  596. {
  597. // Now we use our estimated quotient to update each block of dividend.
  598. // dividend = dividend - divisor * quotient
  599. int index = 0;
  600. ulong borrow = 0;
  601. ulong carry = 0;
  602. do
  603. {
  604. ulong product = ((ulong)(divisor._blocks[index]) * quotient) + carry;
  605. carry = product >> 32;
  606. ulong difference = (ulong)(dividend._blocks[index]) - (uint)(product) - borrow;
  607. borrow = (difference >> 32) & 1;
  608. dividend._blocks[index] = (uint)(difference);
  609. index++;
  610. }
  611. while (index < divisorLength);
  612. // Remove all leading zero blocks from dividend
  613. while ((divisorLength > 0) && (dividend._blocks[divisorLength - 1] == 0))
  614. {
  615. divisorLength--;
  616. }
  617. dividend._length = divisorLength;
  618. }
  619. // If the dividend is still larger than the divisor, we overshot our estimate quotient. To correct,
  620. // we increment the quotient and subtract one more divisor from the dividend (Because we guaranteed the error range).
  621. if (Compare(ref dividend, ref divisor) >= 0)
  622. {
  623. quotient++;
  624. // dividend = dividend - divisor
  625. int index = 0;
  626. ulong borrow = 0;
  627. do
  628. {
  629. ulong difference = (ulong)(dividend._blocks[index]) - divisor._blocks[index] - borrow;
  630. borrow = (difference >> 32) & 1;
  631. dividend._blocks[index] = (uint)(difference);
  632. index++;
  633. }
  634. while (index < divisorLength);
  635. // Remove all leading zero blocks from dividend
  636. while ((divisorLength > 0) && (dividend._blocks[divisorLength - 1] == 0))
  637. {
  638. divisorLength--;
  639. }
  640. dividend._length = divisorLength;
  641. }
  642. return quotient;
  643. }
  644. public static uint LeadingZeroCount(uint value)
  645. {
  646. return 32 - CountSignificantBits(value);
  647. }
  648. public static uint LeadingZeroCount(ulong value)
  649. {
  650. return 64 - CountSignificantBits(value);
  651. }
  652. public static uint LogBase2(uint value)
  653. {
  654. Debug.Assert(value != 0);
  655. // This comes from the Stanford Bit Widdling Hacks by Sean Eron Anderson:
  656. // http://graphics.stanford.edu/~seander/bithacks.html#IntegerLogDeBruijn
  657. value |= (value >> 1); // first round down to one less than a power of 2
  658. value |= (value >> 2);
  659. value |= (value >> 4);
  660. value |= (value >> 8);
  661. value |= (value >> 16);
  662. uint index = (value * 0x07C4ACDD) >> 27;
  663. return s_MultiplyDeBruijnBitPosition[(int)(index)];
  664. }
  665. public static uint LogBase2(ulong value)
  666. {
  667. Debug.Assert(value != 0);
  668. uint upper = (uint)(value >> 32);
  669. if (upper != 0)
  670. {
  671. return 32 + LogBase2(upper);
  672. }
  673. return LogBase2((uint)(value));
  674. }
  675. public static void Multiply(ref BigInteger lhs, uint value, ref BigInteger result)
  676. {
  677. if (lhs.IsZero() || (value == 1))
  678. {
  679. result.SetValue(ref lhs);
  680. return;
  681. }
  682. if (value == 0)
  683. {
  684. result.SetZero();
  685. return;
  686. }
  687. int lhsLength = lhs._length;
  688. int index = 0;
  689. uint carry = 0;
  690. while (index < lhsLength)
  691. {
  692. ulong product = ((ulong)(lhs._blocks[index]) * value) + carry;
  693. result._blocks[index] = (uint)(product);
  694. carry = (uint)(product >> 32);
  695. index++;
  696. }
  697. if (carry != 0)
  698. {
  699. Debug.Assert(unchecked((uint)(lhsLength)) + 1 <= MaxBlockCount);
  700. result._blocks[index] = carry;
  701. result._length = (lhsLength + 1);
  702. }
  703. }
  704. public static void Multiply(ref BigInteger lhs, ref BigInteger rhs, ref BigInteger result)
  705. {
  706. if (lhs.IsZero() || rhs.IsOne())
  707. {
  708. result.SetValue(ref lhs);
  709. return;
  710. }
  711. if (rhs.IsZero())
  712. {
  713. result.SetZero();
  714. return;
  715. }
  716. ref readonly BigInteger large = ref lhs;
  717. int largeLength = lhs._length;
  718. ref readonly BigInteger small = ref rhs;
  719. int smallLength = rhs._length;
  720. if (largeLength < smallLength)
  721. {
  722. large = ref rhs;
  723. largeLength = rhs._length;
  724. small = ref lhs;
  725. smallLength = lhs._length;
  726. }
  727. int maxResultLength = smallLength + largeLength;
  728. Debug.Assert(unchecked((uint)(maxResultLength)) <= MaxBlockCount);
  729. // Zero out result internal blocks.
  730. Buffer.ZeroMemory((byte*)(result.GetBlocksPointer()), (maxResultLength * sizeof(uint)));
  731. int smallIndex = 0;
  732. int resultStartIndex = 0;
  733. while (smallIndex < smallLength)
  734. {
  735. // Multiply each block of large BigNum.
  736. if (small._blocks[smallIndex] != 0)
  737. {
  738. int largeIndex = 0;
  739. int resultIndex = resultStartIndex;
  740. ulong carry = 0;
  741. do
  742. {
  743. ulong product = result._blocks[resultIndex] + ((ulong)(small._blocks[smallIndex]) * large._blocks[largeIndex]) + carry;
  744. carry = product >> 32;
  745. result._blocks[resultIndex] = (uint)(product);
  746. resultIndex++;
  747. largeIndex++;
  748. }
  749. while (largeIndex < largeLength);
  750. result._blocks[resultIndex] = (uint)(carry);
  751. }
  752. smallIndex++;
  753. resultStartIndex++;
  754. }
  755. if ((maxResultLength > 0) && (result._blocks[maxResultLength - 1] == 0))
  756. {
  757. result._length = (maxResultLength - 1);
  758. }
  759. else
  760. {
  761. result._length = maxResultLength;
  762. }
  763. }
  764. public static void Pow2(uint exponent, out BigInteger result)
  765. {
  766. result = new BigInteger(0);
  767. ShiftLeft(1, exponent, ref result);
  768. }
  769. public static void Pow10(uint exponent, out BigInteger result)
  770. {
  771. // We leverage two arrays - s_Pow10UInt32Table and s_Pow10BigNumTable to speed up the Pow10 calculation.
  772. //
  773. // s_Pow10UInt32Table stores the results of 10^0 to 10^7.
  774. // s_Pow10BigNumTable stores the results of 10^8, 10^16, 10^32, 10^64, 10^128, 10^256, and 10^512
  775. //
  776. // For example, let's say exp = 0b111111. We can split the exp to two parts, one is small exp,
  777. // which 10^smallExp can be represented as uint, another part is 10^bigExp, which must be represented as BigNum.
  778. // So the result should be 10^smallExp * 10^bigExp.
  779. //
  780. // Calculating 10^smallExp is simple, we just lookup the 10^smallExp from s_Pow10UInt32Table.
  781. // But here's a bad news: although uint can represent 10^9, exp 9's binary representation is 1001.
  782. // That means 10^(1011), 10^(1101), 10^(1111) all cannot be stored as uint, we cannot easily say something like:
  783. // "Any bits <= 3 is small exp, any bits > 3 is big exp". So instead of involving 10^8, 10^9 to s_Pow10UInt32Table,
  784. // consider 10^8 and 10^9 as a bigNum, so they fall into s_Pow10BigNumTable. Now we can have a simple rule:
  785. // "Any bits <= 3 is small exp, any bits > 3 is big exp".
  786. //
  787. // For 0b111111, we first calculate 10^(smallExp), which is 10^(7), now we can shift right 3 bits, prepare to calculate the bigExp part,
  788. // the exp now becomes 0b000111.
  789. //
  790. // 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.
  791. // Now let's shift exp right 1 bit, the lowest bit should represent 10^(8 * 2) = 10^16, and so on...
  792. //
  793. // That's why we just need the values of s_Pow10BigNumTable be power of 2.
  794. //
  795. // More details of this implementation can be found at: https://github.com/dotnet/coreclr/pull/12894#discussion_r128890596
  796. // Validate that `s_Pow10BigNumTable` has exactly enough trailing elements to fill a BigInteger (which contains MaxBlockCount + 1 elements)
  797. // We validate here, since this is the only current consumer of the array
  798. Debug.Assert((s_Pow10BigNumTableIndices[s_Pow10BigNumTableIndices.Length - 1] + MaxBlockCount + 2) == s_Pow10BigNumTable.Length);
  799. BigInteger temp1 = new BigInteger(s_Pow10UInt32Table[exponent & 0x7]);
  800. ref BigInteger lhs = ref temp1;
  801. BigInteger temp2 = new BigInteger(0);
  802. ref BigInteger product = ref temp2;
  803. exponent >>= 3;
  804. uint index = 0;
  805. while (exponent != 0)
  806. {
  807. // If the current bit is set, multiply it with the corresponding power of 10
  808. if ((exponent & 1) != 0)
  809. {
  810. // Multiply into the next temporary
  811. fixed (uint* pBigNumEntry = &s_Pow10BigNumTable[s_Pow10BigNumTableIndices[index]])
  812. {
  813. ref BigInteger rhs = ref *(BigInteger*)(pBigNumEntry);
  814. Multiply(ref lhs, ref rhs, ref product);
  815. }
  816. // Swap to the next temporary
  817. ref BigInteger temp = ref product;
  818. product = ref lhs;
  819. lhs = ref temp;
  820. }
  821. // Advance to the next bit
  822. ++index;
  823. exponent >>= 1;
  824. }
  825. result = new BigInteger(0);
  826. result.SetValue(ref lhs);
  827. }
  828. public static void ShiftLeft(ulong input, uint shift, ref BigInteger output)
  829. {
  830. if (shift == 0)
  831. {
  832. return;
  833. }
  834. uint blocksToShift = Math.DivRem(shift, 32, out uint remainingBitsToShift);
  835. if (blocksToShift > 0)
  836. {
  837. // If blocks shifted, we should fill the corresponding blocks with zero.
  838. output.ExtendBlocks(0, blocksToShift);
  839. }
  840. if (remainingBitsToShift == 0)
  841. {
  842. // We shift 32 * n (n >= 1) bits. No remaining bits.
  843. output.ExtendBlock((uint)(input));
  844. uint highBits = (uint)(input >> 32);
  845. if (highBits != 0)
  846. {
  847. output.ExtendBlock(highBits);
  848. }
  849. }
  850. else
  851. {
  852. // Extract the high position bits which would be shifted out of range.
  853. uint highPositionBits = (uint)(input) >> (int)(64 - remainingBitsToShift);
  854. // Shift the input. The result should be stored to current block.
  855. ulong shiftedInput = input << (int)(remainingBitsToShift);
  856. output.ExtendBlock((uint)(shiftedInput));
  857. uint highBits = (uint)(input >> 32);
  858. if (highBits != 0)
  859. {
  860. output.ExtendBlock(highBits);
  861. }
  862. if (highPositionBits != 0)
  863. {
  864. // If the high position bits is not 0, we should store them to next block.
  865. output.ExtendBlock(highPositionBits);
  866. }
  867. }
  868. }
  869. private static unsafe uint AddDivisor(ref BigInteger lhs, int lhsStartIndex, ref BigInteger rhs)
  870. {
  871. int lhsLength = lhs._length;
  872. int rhsLength = rhs._length;
  873. Debug.Assert(lhsLength >= 0);
  874. Debug.Assert(rhsLength >= 0);
  875. Debug.Assert(lhsLength >= rhsLength);
  876. // Repairs the dividend, if the last subtract was too much
  877. ulong carry = 0UL;
  878. for (int i = 0; i < rhsLength; i++)
  879. {
  880. ref uint lhsValue = ref lhs._blocks[lhsStartIndex + i];
  881. ulong digit = (lhsValue + carry) + rhs._blocks[i];
  882. lhsValue = unchecked((uint)digit);
  883. carry = digit >> 32;
  884. }
  885. return (uint)(carry);
  886. }
  887. private static bool DivideGuessTooBig(ulong q, ulong valHi, uint valLo, uint divHi, uint divLo)
  888. {
  889. Debug.Assert(q <= 0xFFFFFFFF);
  890. // We multiply the two most significant limbs of the divisor
  891. // with the current guess for the quotient. If those are bigger
  892. // than the three most significant limbs of the current dividend
  893. // we return true, which means the current guess is still too big.
  894. ulong chkHi = divHi * q;
  895. ulong chkLo = divLo * q;
  896. chkHi = chkHi + (chkLo >> 32);
  897. chkLo = chkLo & uint.MaxValue;
  898. if (chkHi < valHi)
  899. return false;
  900. if (chkHi > valHi)
  901. return true;
  902. if (chkLo < valLo)
  903. return false;
  904. if (chkLo > valLo)
  905. return true;
  906. return false;
  907. }
  908. private static unsafe uint SubtractDivisor(ref BigInteger lhs, int lhsStartIndex, ref BigInteger rhs, ulong q)
  909. {
  910. int lhsLength = lhs._length - lhsStartIndex;
  911. int rhsLength = rhs._length;
  912. Debug.Assert((lhsLength) >= 0);
  913. Debug.Assert(rhsLength >= 0);
  914. Debug.Assert(lhsLength >= rhsLength);
  915. Debug.Assert(q <= uint.MaxValue);
  916. // Combines a subtract and a multiply operation, which is naturally
  917. // more efficient than multiplying and then subtracting...
  918. ulong carry = 0;
  919. for (int i = 0; i < rhsLength; i++)
  920. {
  921. carry += rhs._blocks[i] * q;
  922. uint digit = unchecked((uint)carry);
  923. carry = carry >> 32;
  924. ref uint lhsValue = ref lhs._blocks[lhsStartIndex + i];
  925. if (lhsValue < digit)
  926. {
  927. carry++;
  928. }
  929. lhsValue = unchecked(lhsValue - digit);
  930. }
  931. return (uint)(carry);
  932. }
  933. public void Add(uint value)
  934. {
  935. Add(ref this, value, ref this);
  936. }
  937. public void ExtendBlock(uint blockValue)
  938. {
  939. _blocks[_length] = blockValue;
  940. _length++;
  941. }
  942. public void ExtendBlocks(uint blockValue, uint blockCount)
  943. {
  944. Debug.Assert(blockCount > 0);
  945. if (blockCount == 1)
  946. {
  947. ExtendBlock(blockValue);
  948. return;
  949. }
  950. Buffer.ZeroMemory((byte*)(GetBlocksPointer() + _length), ((blockCount - 1) * sizeof(uint)));
  951. _length += (int)(blockCount);
  952. _blocks[_length - 1] = blockValue;
  953. }
  954. public uint GetBlock(uint index)
  955. {
  956. Debug.Assert(index < _length);
  957. return _blocks[index];
  958. }
  959. public int GetLength()
  960. {
  961. return _length;
  962. }
  963. public bool IsOne()
  964. {
  965. return (_length == 1)
  966. && (_blocks[0] == 1);
  967. }
  968. public bool IsZero()
  969. {
  970. return _length == 0;
  971. }
  972. public void Multiply(uint value)
  973. {
  974. Multiply(ref this, value, ref this);
  975. }
  976. public void Multiply(ref BigInteger value)
  977. {
  978. var result = new BigInteger(0);
  979. Multiply(ref this, ref value, ref result);
  980. Buffer.Memcpy((byte*)(GetBlocksPointer()), (byte*)(result.GetBlocksPointer()), (result._length) * sizeof(uint));
  981. _length = result._length;
  982. }
  983. public void Multiply10()
  984. {
  985. if (IsZero())
  986. {
  987. return;
  988. }
  989. int index = 0;
  990. int length = _length;
  991. ulong carry = 0;
  992. while (index < length)
  993. {
  994. var block = (ulong)(_blocks[index]);
  995. ulong product = (block << 3) + (block << 1) + carry;
  996. carry = product >> 32;
  997. _blocks[index] = (uint)(product);
  998. index++;
  999. }
  1000. if (carry != 0)
  1001. {
  1002. Debug.Assert(unchecked((uint)(_length)) + 1 <= MaxBlockCount);
  1003. _blocks[index] = (uint)(carry);
  1004. _length += 1;
  1005. }
  1006. }
  1007. public void MultiplyPow10(uint exponent)
  1008. {
  1009. Pow10(exponent, out BigInteger poweredValue);
  1010. if (poweredValue._length == 1)
  1011. {
  1012. Multiply(poweredValue._blocks[0]);
  1013. }
  1014. else
  1015. {
  1016. Multiply(ref poweredValue);
  1017. }
  1018. }
  1019. public void SetUInt32(uint value)
  1020. {
  1021. _blocks[0] = value;
  1022. _length = 1;
  1023. }
  1024. public void SetUInt64(ulong value)
  1025. {
  1026. var lower = (uint)(value);
  1027. var upper = (uint)(value >> 32);
  1028. _blocks[0] = lower;
  1029. _blocks[1] = upper;
  1030. _length = (upper == 0) ? 1 : 2;
  1031. }
  1032. public void SetValue(ref BigInteger rhs)
  1033. {
  1034. int rhsLength = rhs._length;
  1035. Buffer.Memcpy((byte*)(GetBlocksPointer()), (byte*)(rhs.GetBlocksPointer()), (rhsLength * sizeof(uint)));
  1036. _length = rhsLength;
  1037. }
  1038. public void SetZero()
  1039. {
  1040. _length = 0;
  1041. }
  1042. public void ShiftLeft(uint shift)
  1043. {
  1044. // Process blocks high to low so that we can safely process in place
  1045. var length = _length;
  1046. if ((length == 0) || (shift == 0))
  1047. {
  1048. return;
  1049. }
  1050. uint blocksToShift = Math.DivRem(shift, 32, out uint remainingBitsToShift);
  1051. // Copy blocks from high to low
  1052. int readIndex = (length - 1);
  1053. int writeIndex = readIndex + (int)(blocksToShift);
  1054. // Check if the shift is block aligned
  1055. if (remainingBitsToShift == 0)
  1056. {
  1057. Debug.Assert(writeIndex < MaxBlockCount);
  1058. while (readIndex >= 0)
  1059. {
  1060. _blocks[writeIndex] = _blocks[readIndex];
  1061. readIndex--;
  1062. writeIndex--;
  1063. }
  1064. _length += (int)(blocksToShift);
  1065. // Zero the remaining low blocks
  1066. Buffer.ZeroMemory((byte*)(GetBlocksPointer()), (blocksToShift * sizeof(uint)));
  1067. }
  1068. else
  1069. {
  1070. // We need an extra block for the partial shift
  1071. writeIndex++;
  1072. Debug.Assert(writeIndex < MaxBlockCount);
  1073. // Set the length to hold the shifted blocks
  1074. _length = writeIndex + 1;
  1075. // Output the initial blocks
  1076. uint lowBitsShift = (32 - remainingBitsToShift);
  1077. uint highBits = 0;
  1078. uint block = _blocks[readIndex];
  1079. uint lowBits = block >> (int)(lowBitsShift);
  1080. while (readIndex > 0)
  1081. {
  1082. _blocks[writeIndex] = highBits | lowBits;
  1083. highBits = block << (int)(remainingBitsToShift);
  1084. --readIndex;
  1085. --writeIndex;
  1086. block = _blocks[readIndex];
  1087. lowBits = block >> (int)lowBitsShift;
  1088. }
  1089. // Output the final blocks
  1090. _blocks[writeIndex] = highBits | lowBits;
  1091. _blocks[writeIndex - 1] = block << (int)(remainingBitsToShift);
  1092. // Zero the remaining low blocks
  1093. Buffer.ZeroMemory((byte*)(GetBlocksPointer()), (blocksToShift * sizeof(uint)));
  1094. // Check if the terminating block has no set bits
  1095. if (_blocks[_length - 1] == 0)
  1096. {
  1097. _length--;
  1098. }
  1099. }
  1100. }
  1101. public ulong ToUInt64()
  1102. {
  1103. if (_length > 1)
  1104. {
  1105. return ((ulong)(_blocks[1]) << 32) + _blocks[0];
  1106. }
  1107. if (_length > 0)
  1108. {
  1109. return _blocks[0];
  1110. }
  1111. return 0;
  1112. }
  1113. private uint* GetBlocksPointer()
  1114. {
  1115. // This is safe to do since we are a ref struct
  1116. return (uint*)(Unsafe.AsPointer(ref _blocks[0]));
  1117. }
  1118. }
  1119. }
  1120. }