PageRenderTime 64ms CodeModel.GetById 26ms RepoModel.GetById 0ms app.codeStats 1ms

/SSISSSHFTP/SSISSharpSSHFTP/DilffieHellman/mono/BigInteger.cs

#
C# | 2257 lines | 1532 code | 477 blank | 248 comment | 421 complexity | 95186552756b601e33603a4d74a7684c MD5 | raw file

Large files files are truncated, but you can click here to view the full file

  1. //
  2. // BigInteger.cs - Big Integer implementation
  3. //
  4. // Authors:
  5. // Ben Maurer
  6. // Chew Keong TAN
  7. // Sebastien Pouliot (spouliot@motus.com)
  8. //
  9. // Copyright (c) 2003 Ben Maurer
  10. // All rights reserved
  11. //
  12. // Copyright (c) 2002 Chew Keong TAN
  13. // All rights reserved.
  14. using System;
  15. using System.Security.Cryptography;
  16. using Mono.Math.Prime.Generator;
  17. using Mono.Math.Prime;
  18. namespace Mono.Math {
  19. [CLSCompliant(false)]
  20. internal class BigInteger {
  21. #region Data Storage
  22. /// <summary>
  23. /// The Length of this BigInteger
  24. /// </summary>
  25. uint length = 1;
  26. /// <summary>
  27. /// The data for this BigInteger
  28. /// </summary>
  29. uint [] data;
  30. #endregion
  31. #region Constants
  32. /// <summary>
  33. /// Default length of a BigInteger in bytes
  34. /// </summary>
  35. const uint DEFAULT_LEN = 20;
  36. /// <summary>
  37. /// Table of primes below 2000.
  38. /// </summary>
  39. /// <remarks>
  40. /// <para>
  41. /// This table was generated using Mathematica 4.1 using the following function:
  42. /// </para>
  43. /// <para>
  44. /// <code>
  45. /// PrimeTable [x_] := Prime [Range [1, PrimePi [x]]]
  46. /// PrimeTable [6000]
  47. /// </code>
  48. /// </para>
  49. /// </remarks>
  50. public static readonly uint [] smallPrimes = {
  51. 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
  52. 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
  53. 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233,
  54. 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317,
  55. 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419,
  56. 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503,
  57. 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607,
  58. 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701,
  59. 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811,
  60. 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911,
  61. 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997,
  62. 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087,
  63. 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181,
  64. 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279,
  65. 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373,
  66. 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471,
  67. 1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559,
  68. 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637,
  69. 1657, 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747,
  70. 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861, 1867,
  71. 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973,
  72. 1979, 1987, 1993, 1997, 1999,
  73. 2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087, 2089,
  74. 2099, 2111, 2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207,
  75. 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, 2293, 2297,
  76. 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377, 2381, 2383, 2389,
  77. 2393, 2399, 2411, 2417, 2423, 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503,
  78. 2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, 2621,
  79. 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707,
  80. 2711, 2713, 2719, 2729, 2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791, 2797,
  81. 2801, 2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903,
  82. 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999,
  83. 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, 3083, 3089, 3109,
  84. 3119, 3121, 3137, 3163, 3167, 3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221,
  85. 3229, 3251, 3253, 3257, 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329,
  86. 3331, 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, 3433, 3449,
  87. 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533, 3539,
  88. 3541, 3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631,
  89. 3637, 3643, 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, 3733,
  90. 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, 3823, 3833, 3847, 3851,
  91. 3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923, 3929, 3931, 3943,
  92. 3947, 3967, 3989,
  93. 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, 4073, 4079, 4091,
  94. 4093, 4099, 4111, 4127, 4129, 4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211,
  95. 4217, 4219, 4229, 4231, 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289,
  96. 4297, 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, 4421, 4423,
  97. 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, 4507, 4513, 4517, 4519, 4523,
  98. 4547, 4549, 4561, 4567, 4583, 4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649,
  99. 4651, 4657, 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, 4759,
  100. 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, 4861, 4871, 4877, 4889,
  101. 4903, 4909, 4919, 4931, 4933, 4937, 4943, 4951, 4957, 4967, 4969, 4973, 4987,
  102. 4993, 4999,
  103. 5003, 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, 5099, 5101,
  104. 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179, 5189, 5197, 5209, 5227, 5231,
  105. 5233, 5237, 5261, 5273, 5279, 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351,
  106. 5381, 5387, 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, 5449,
  107. 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, 5527, 5531, 5557, 5563,
  108. 5569, 5573, 5581, 5591, 5623, 5639, 5641, 5647, 5651, 5653, 5657, 5659, 5669,
  109. 5683, 5689, 5693, 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791,
  110. 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857, 5861, 5867, 5869,
  111. 5879, 5881, 5897, 5903, 5923, 5927, 5939, 5953, 5981, 5987
  112. };
  113. public enum Sign : int {
  114. Negative = -1,
  115. Zero = 0,
  116. Positive = 1
  117. };
  118. #region Exception Messages
  119. const string WouldReturnNegVal = "Operation would return a negative value";
  120. #endregion
  121. #endregion
  122. #region Constructors
  123. public BigInteger ()
  124. {
  125. data = new uint [DEFAULT_LEN];
  126. }
  127. public BigInteger (Sign sign, uint len)
  128. {
  129. this.data = new uint [len];
  130. this.length = len;
  131. }
  132. public BigInteger (BigInteger bi)
  133. {
  134. this.data = (uint [])bi.data.Clone ();
  135. this.length = bi.length;
  136. }
  137. public BigInteger (BigInteger bi, uint len)
  138. {
  139. this.data = new uint [len];
  140. for (uint i = 0; i < bi.length; i++)
  141. this.data [i] = bi.data [i];
  142. this.length = bi.length;
  143. }
  144. #endregion
  145. public static BigInteger Parse(string number) {
  146. if (number == null)
  147. throw new ArgumentNullException(number);
  148. int i = 0, len = number.Length;
  149. char c;
  150. bool digits_seen = false;
  151. BigInteger val = new BigInteger(0);
  152. if (number[i] == '+') {
  153. i++;
  154. } else if(number[i] == '-') {
  155. throw new FormatException("Only positive integers are allowed.");
  156. }
  157. for(; i < len; i++) {
  158. c = number[i];
  159. if (c == '\0') {
  160. i = len;
  161. continue;
  162. }
  163. if (c >= '0' && c <= '9'){
  164. val = val * 10 + (c - '0');
  165. digits_seen = true;
  166. } else {
  167. if (Char.IsWhiteSpace(c)){
  168. for (i++; i < len; i++){
  169. if (!Char.IsWhiteSpace (number[i]))
  170. throw new FormatException();
  171. }
  172. break;
  173. } else
  174. throw new FormatException();
  175. }
  176. }
  177. if (!digits_seen)
  178. throw new FormatException();
  179. return val;
  180. }
  181. #region Conversions
  182. public BigInteger (byte [] inData)
  183. {
  184. length = (uint)inData.Length >> 2;
  185. int leftOver = inData.Length & 0x3;
  186. // length not multiples of 4
  187. if (leftOver != 0) length++;
  188. data = new uint [length];
  189. for (int i = inData.Length - 1, j = 0; i >= 3; i -= 4, j++) {
  190. data [j] = (uint)(
  191. (inData [i-3] << (3*8)) |
  192. (inData [i-2] << (2*8)) |
  193. (inData [i-1] << (1*8)) |
  194. (inData [i-0] << (0*8))
  195. );
  196. }
  197. switch (leftOver) {
  198. case 1: data [length-1] = (uint)inData [0]; break;
  199. case 2: data [length-1] = (uint)((inData [0] << 8) | inData [1]); break;
  200. case 3: data [length-1] = (uint)((inData [0] << 16) | (inData [1] << 8) | inData [2]); break;
  201. }
  202. this.Normalize ();
  203. }
  204. public BigInteger (uint [] inData)
  205. {
  206. length = (uint)inData.Length;
  207. data = new uint [length];
  208. for (int i = (int)length - 1, j = 0; i >= 0; i--, j++)
  209. data [j] = inData [i];
  210. this.Normalize ();
  211. }
  212. public BigInteger (uint ui)
  213. {
  214. data = new uint [] {ui};
  215. }
  216. public BigInteger (ulong ul)
  217. {
  218. data = new uint [2] { (uint)ul, (uint)(ul >> 32)};
  219. length = 2;
  220. this.Normalize ();
  221. }
  222. public static implicit operator BigInteger (uint value)
  223. {
  224. return (new BigInteger (value));
  225. }
  226. public static implicit operator BigInteger (int value)
  227. {
  228. if (value < 0) throw new ArgumentOutOfRangeException ("value");
  229. return (new BigInteger ((uint)value));
  230. }
  231. public static implicit operator BigInteger (ulong value)
  232. {
  233. return (new BigInteger (value));
  234. }
  235. #endregion
  236. #region Operators
  237. public static BigInteger operator + (BigInteger bi1, BigInteger bi2)
  238. {
  239. if (bi1 == 0)
  240. return new BigInteger (bi2);
  241. else if (bi2 == 0)
  242. return new BigInteger (bi1);
  243. else
  244. return Kernel.AddSameSign (bi1, bi2);
  245. }
  246. public static BigInteger operator - (BigInteger bi1, BigInteger bi2)
  247. {
  248. if (bi2 == 0)
  249. return new BigInteger (bi1);
  250. if (bi1 == 0)
  251. throw new ArithmeticException (WouldReturnNegVal);
  252. switch (Kernel.Compare (bi1, bi2)) {
  253. case Sign.Zero:
  254. return 0;
  255. case Sign.Positive:
  256. return Kernel.Subtract (bi1, bi2);
  257. case Sign.Negative:
  258. throw new ArithmeticException (WouldReturnNegVal);
  259. default:
  260. throw new Exception ();
  261. }
  262. }
  263. public static int operator % (BigInteger bi, int i)
  264. {
  265. if (i > 0)
  266. return (int)Kernel.DwordMod (bi, (uint)i);
  267. else
  268. return -(int)Kernel.DwordMod (bi, (uint)-i);
  269. }
  270. public static uint operator % (BigInteger bi, uint ui)
  271. {
  272. return Kernel.DwordMod (bi, (uint)ui);
  273. }
  274. public static BigInteger operator % (BigInteger bi1, BigInteger bi2)
  275. {
  276. return Kernel.multiByteDivide (bi1, bi2)[1];
  277. }
  278. public static BigInteger operator / (BigInteger bi, int i)
  279. {
  280. if (i > 0)
  281. return Kernel.DwordDiv (bi, (uint)i);
  282. throw new ArithmeticException (WouldReturnNegVal);
  283. }
  284. public static BigInteger operator / (BigInteger bi1, BigInteger bi2)
  285. {
  286. return Kernel.multiByteDivide (bi1, bi2)[0];
  287. }
  288. public static BigInteger operator * (BigInteger bi1, BigInteger bi2)
  289. {
  290. if (bi1 == 0 || bi2 == 0) return 0;
  291. //
  292. // Validate pointers
  293. //
  294. if (bi1.data.Length < bi1.length) throw new IndexOutOfRangeException ("bi1 out of range");
  295. if (bi2.data.Length < bi2.length) throw new IndexOutOfRangeException ("bi2 out of range");
  296. BigInteger ret = new BigInteger (Sign.Positive, bi1.length + bi2.length);
  297. Kernel.Multiply (bi1.data, 0, bi1.length, bi2.data, 0, bi2.length, ret.data, 0);
  298. ret.Normalize ();
  299. return ret;
  300. }
  301. public static BigInteger operator * (BigInteger bi, int i)
  302. {
  303. if (i < 0) throw new ArithmeticException (WouldReturnNegVal);
  304. if (i == 0) return 0;
  305. if (i == 1) return new BigInteger (bi);
  306. return Kernel.MultiplyByDword (bi, (uint)i);
  307. }
  308. public static BigInteger operator << (BigInteger bi1, int shiftVal)
  309. {
  310. return Kernel.LeftShift (bi1, shiftVal);
  311. }
  312. public static BigInteger operator >> (BigInteger bi1, int shiftVal)
  313. {
  314. return Kernel.RightShift (bi1, shiftVal);
  315. }
  316. #endregion
  317. #region Random
  318. private static RandomNumberGenerator rng;
  319. private static RandomNumberGenerator Rng {
  320. get {
  321. if (rng == null)
  322. rng = RandomNumberGenerator.Create ();
  323. return rng;
  324. }
  325. }
  326. /// <summary>
  327. /// Generates a new, random BigInteger of the specified length.
  328. /// </summary>
  329. /// <param name="bits">The number of bits for the new number.</param>
  330. /// <param name="rng">A random number generator to use to obtain the bits.</param>
  331. /// <returns>A random number of the specified length.</returns>
  332. public static BigInteger genRandom (int bits, RandomNumberGenerator rng)
  333. {
  334. int dwords = bits >> 5;
  335. int remBits = bits & 0x1F;
  336. if (remBits != 0)
  337. dwords++;
  338. BigInteger ret = new BigInteger (Sign.Positive, (uint)dwords + 1);
  339. byte [] random = new byte [dwords << 2];
  340. rng.GetBytes (random);
  341. Buffer.BlockCopy (random, 0, ret.data, 0, (int)dwords << 2);
  342. if (remBits != 0) {
  343. uint mask = (uint)(0x01 << (remBits-1));
  344. ret.data [dwords-1] |= mask;
  345. mask = (uint)(0xFFFFFFFF >> (32 - remBits));
  346. ret.data [dwords-1] &= mask;
  347. }
  348. else
  349. ret.data [dwords-1] |= 0x80000000;
  350. ret.Normalize ();
  351. return ret;
  352. }
  353. /// <summary>
  354. /// Generates a new, random BigInteger of the specified length using the default RNG crypto service provider.
  355. /// </summary>
  356. /// <param name="bits">The number of bits for the new number.</param>
  357. /// <returns>A random number of the specified length.</returns>
  358. public static BigInteger genRandom (int bits)
  359. {
  360. return genRandom (bits, Rng);
  361. }
  362. /// <summary>
  363. /// Randomizes the bits in "this" from the specified RNG.
  364. /// </summary>
  365. /// <param name="rng">A RNG.</param>
  366. public void randomize (RandomNumberGenerator rng)
  367. {
  368. int bits = this.bitCount ();
  369. int dwords = bits >> 5;
  370. int remBits = bits & 0x1F;
  371. if (remBits != 0)
  372. dwords++;
  373. byte [] random = new byte [dwords << 2];
  374. rng.GetBytes (random);
  375. Buffer.BlockCopy (random, 0, data, 0, (int)dwords << 2);
  376. if (remBits != 0) {
  377. uint mask = (uint)(0x01 << (remBits-1));
  378. data [dwords-1] |= mask;
  379. mask = (uint)(0xFFFFFFFF >> (32 - remBits));
  380. data [dwords-1] &= mask;
  381. }
  382. else
  383. data [dwords-1] |= 0x80000000;
  384. Normalize ();
  385. }
  386. /// <summary>
  387. /// Randomizes the bits in "this" from the default RNG.
  388. /// </summary>
  389. public void randomize ()
  390. {
  391. randomize (Rng);
  392. }
  393. #endregion
  394. #region Bitwise
  395. public int bitCount ()
  396. {
  397. this.Normalize ();
  398. uint value = data [length - 1];
  399. uint mask = 0x80000000;
  400. uint bits = 32;
  401. while (bits > 0 && (value & mask) == 0) {
  402. bits--;
  403. mask >>= 1;
  404. }
  405. bits += ((length - 1) << 5);
  406. return (int)bits;
  407. }
  408. /// <summary>
  409. /// Tests if the specified bit is 1.
  410. /// </summary>
  411. /// <param name="bitNum">The bit to test. The least significant bit is 0.</param>
  412. /// <returns>True if bitNum is set to 1, else false.</returns>
  413. public bool testBit (uint bitNum)
  414. {
  415. uint bytePos = bitNum >> 5; // divide by 32
  416. byte bitPos = (byte)(bitNum & 0x1F); // get the lowest 5 bits
  417. uint mask = (uint)1 << bitPos;
  418. return ((this.data [bytePos] & mask) != 0);
  419. }
  420. public bool testBit (int bitNum)
  421. {
  422. if (bitNum < 0) throw new IndexOutOfRangeException ("bitNum out of range");
  423. uint bytePos = (uint)bitNum >> 5; // divide by 32
  424. byte bitPos = (byte)(bitNum & 0x1F); // get the lowest 5 bits
  425. uint mask = (uint)1 << bitPos;
  426. return ((this.data [bytePos] | mask) == this.data [bytePos]);
  427. }
  428. public void setBit (uint bitNum)
  429. {
  430. setBit (bitNum, true);
  431. }
  432. public void clearBit (uint bitNum)
  433. {
  434. setBit (bitNum, false);
  435. }
  436. public void setBit (uint bitNum, bool val)
  437. {
  438. uint bytePos = bitNum >> 5; // divide by 32
  439. if (bytePos < this.length) {
  440. uint mask = (uint)1 << (int)(bitNum & 0x1F);
  441. if (val)
  442. this.data [bytePos] |= mask;
  443. else
  444. this.data [bytePos] &= ~mask;
  445. }
  446. }
  447. public int LowestSetBit ()
  448. {
  449. if (this == 0) return -1;
  450. int i = 0;
  451. while (!testBit (i)) i++;
  452. return i;
  453. }
  454. public byte [] getBytes ()
  455. {
  456. if (this == 0) return new byte [1];
  457. int numBits = bitCount ();
  458. int numBytes = numBits >> 3;
  459. if ((numBits & 0x7) != 0)
  460. numBytes++;
  461. byte [] result = new byte [numBytes];
  462. int numBytesInWord = numBytes & 0x3;
  463. if (numBytesInWord == 0) numBytesInWord = 4;
  464. int pos = 0;
  465. for (int i = (int)length - 1; i >= 0; i--) {
  466. uint val = data [i];
  467. for (int j = numBytesInWord - 1; j >= 0; j--) {
  468. result [pos+j] = (byte)(val & 0xFF);
  469. val >>= 8;
  470. }
  471. pos += numBytesInWord;
  472. numBytesInWord = 4;
  473. }
  474. return result;
  475. }
  476. #endregion
  477. #region Compare
  478. public static bool operator == (BigInteger bi1, uint ui)
  479. {
  480. if (bi1.length != 1) bi1.Normalize ();
  481. return bi1.length == 1 && bi1.data [0] == ui;
  482. }
  483. public static bool operator != (BigInteger bi1, uint ui)
  484. {
  485. if (bi1.length != 1) bi1.Normalize ();
  486. return !(bi1.length == 1 && bi1.data [0] == ui);
  487. }
  488. public static bool operator == (BigInteger bi1, BigInteger bi2)
  489. {
  490. // we need to compare with null
  491. if ((bi1 as object) == (bi2 as object))
  492. return true;
  493. if (null == bi1 || null == bi2)
  494. return false;
  495. return Kernel.Compare (bi1, bi2) == 0;
  496. }
  497. public static bool operator != (BigInteger bi1, BigInteger bi2)
  498. {
  499. // we need to compare with null
  500. if ((bi1 as object) == (bi2 as object))
  501. return false;
  502. if (null == bi1 || null == bi2)
  503. return true;
  504. return Kernel.Compare (bi1, bi2) != 0;
  505. }
  506. public static bool operator > (BigInteger bi1, BigInteger bi2)
  507. {
  508. return Kernel.Compare (bi1, bi2) > 0;
  509. }
  510. public static bool operator < (BigInteger bi1, BigInteger bi2)
  511. {
  512. return Kernel.Compare (bi1, bi2) < 0;
  513. }
  514. public static bool operator >= (BigInteger bi1, BigInteger bi2)
  515. {
  516. return Kernel.Compare (bi1, bi2) >= 0;
  517. }
  518. public static bool operator <= (BigInteger bi1, BigInteger bi2)
  519. {
  520. return Kernel.Compare (bi1, bi2) <= 0;
  521. }
  522. public Sign Compare (BigInteger bi)
  523. {
  524. return Kernel.Compare (this, bi);
  525. }
  526. #endregion
  527. #region Formatting
  528. public string ToString (uint radix)
  529. {
  530. return ToString (radix, "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ");
  531. }
  532. public string ToString (uint radix, string charSet)
  533. {
  534. if (charSet.Length < radix)
  535. throw new ArgumentException ("charSet length less than radix", "charSet");
  536. if (radix == 1)
  537. throw new ArgumentException ("There is no such thing as radix one notation", "radix");
  538. if (this == 0) return "0";
  539. if (this == 1) return "1";
  540. string result = "";
  541. BigInteger a = new BigInteger (this);
  542. while (a != 0) {
  543. uint rem = Kernel.SingleByteDivideInPlace (a, radix);
  544. result = charSet [ (int)rem] + result;
  545. }
  546. return result;
  547. }
  548. #endregion
  549. #region Misc
  550. /// <summary>
  551. /// Normalizes this by setting the length to the actual number of
  552. /// uints used in data and by setting the sign to Sign.Zero if the
  553. /// value of this is 0.
  554. /// </summary>
  555. private void Normalize ()
  556. {
  557. // Normalize length
  558. while (length > 0 && data [length-1] == 0) length--;
  559. // Check for zero
  560. if (length == 0)
  561. length++;
  562. }
  563. public void Clear ()
  564. {
  565. for (int i=0; i < length; i++)
  566. data [i] = 0x00;
  567. }
  568. #endregion
  569. #region Object Impl
  570. public override int GetHashCode ()
  571. {
  572. uint val = 0;
  573. for (uint i = 0; i < this.length; i++)
  574. val ^= this.data [i];
  575. return (int)val;
  576. }
  577. public override string ToString ()
  578. {
  579. return ToString (10);
  580. }
  581. public override bool Equals (object o)
  582. {
  583. if (o == null) return false;
  584. if (o is int) return (int)o >= 0 && this == (uint)o;
  585. return Kernel.Compare (this, (BigInteger)o) == 0;
  586. }
  587. #endregion
  588. #region Number Theory
  589. public BigInteger gcd (BigInteger bi)
  590. {
  591. return Kernel.gcd (this, bi);
  592. }
  593. public BigInteger modInverse (BigInteger mod)
  594. {
  595. return Kernel.modInverse (this, mod);
  596. }
  597. public BigInteger modPow (BigInteger exp, BigInteger n)
  598. {
  599. ModulusRing mr = new ModulusRing (n);
  600. return mr.Pow (this, exp);
  601. }
  602. #endregion
  603. #region Prime Testing
  604. public bool isProbablePrime ()
  605. {
  606. for (int p = 0; p < smallPrimes.Length; p++) {
  607. if (this == smallPrimes [p])
  608. return true;
  609. if (this % smallPrimes [p] == 0)
  610. return false;
  611. }
  612. return PrimalityTests.RabinMillerTest (this, Prime.ConfidenceFactor.Medium);
  613. }
  614. [Obsolete]
  615. public bool isProbablePrime (int notUsed)
  616. {
  617. for (int p = 0; p < smallPrimes.Length; p++) {
  618. if (this % smallPrimes [p] == 0)
  619. return false;
  620. }
  621. return
  622. PrimalityTests.SmallPrimeSppTest (this, Prime.ConfidenceFactor.Medium);
  623. }
  624. #endregion
  625. #region Prime Number Generation
  626. /// <summary>
  627. /// Generates the smallest prime >= bi
  628. /// </summary>
  629. /// <param name="bi">A BigInteger</param>
  630. /// <returns>The smallest prime >= bi. More mathematically, if bi is prime: bi, else Prime [PrimePi [bi] + 1].</returns>
  631. public static BigInteger NextHightestPrime (BigInteger bi)
  632. {
  633. NextPrimeFinder npf = new NextPrimeFinder ();
  634. return npf.GenerateNewPrime (0, bi);
  635. }
  636. public static BigInteger genPseudoPrime (int bits)
  637. {
  638. SequentialSearchPrimeGeneratorBase sspg = new SequentialSearchPrimeGeneratorBase ();
  639. return sspg.GenerateNewPrime (bits);
  640. }
  641. /// <summary>
  642. /// Increments this by two
  643. /// </summary>
  644. public void Incr2 ()
  645. {
  646. int i = 0;
  647. data [0] += 2;
  648. // If there was no carry, nothing to do
  649. if (data [0] < 2) {
  650. // Account for the first carry
  651. data [++i]++;
  652. // Keep adding until no carry
  653. while (data [i++] == 0x0)
  654. data [i]++;
  655. // See if we increased the data length
  656. if (length == (uint)i)
  657. length++;
  658. }
  659. }
  660. #endregion
  661. public sealed class ModulusRing {
  662. BigInteger mod, constant;
  663. public ModulusRing (BigInteger mod)
  664. {
  665. this.mod = mod;
  666. // calculate constant = b^ (2k) / m
  667. uint i = mod.length << 1;
  668. constant = new BigInteger (Sign.Positive, i + 1);
  669. constant.data [i] = 0x00000001;
  670. constant = constant / mod;
  671. }
  672. public void BarrettReduction (BigInteger x)
  673. {
  674. BigInteger n = mod;
  675. uint k = n.length,
  676. kPlusOne = k+1,
  677. kMinusOne = k-1;
  678. // x < mod, so nothing to do.
  679. if (x.length < k) return;
  680. BigInteger q3;
  681. //
  682. // Validate pointers
  683. //
  684. if (x.data.Length < x.length) throw new IndexOutOfRangeException ("x out of range");
  685. // q1 = x / b^ (k-1)
  686. // q2 = q1 * constant
  687. // q3 = q2 / b^ (k+1), Needs to be accessed with an offset of kPlusOne
  688. // TODO: We should the method in HAC p 604 to do this (14.45)
  689. q3 = new BigInteger (Sign.Positive, x.length - kMinusOne + constant.length);
  690. Kernel.Multiply (x.data, kMinusOne, x.length - kMinusOne, constant.data, 0, constant.length, q3.data, 0);
  691. // r1 = x mod b^ (k+1)
  692. // i.e. keep the lowest (k+1) words
  693. uint lengthToCopy = (x.length > kPlusOne) ? kPlusOne : x.length;
  694. x.length = lengthToCopy;
  695. x.Normalize ();
  696. // r2 = (q3 * n) mod b^ (k+1)
  697. // partial multiplication of q3 and n
  698. BigInteger r2 = new BigInteger (Sign.Positive, kPlusOne);
  699. Kernel.MultiplyMod2p32pmod (q3.data, (int)kPlusOne, (int)q3.length - (int)kPlusOne, n.data, 0, (int)n.length, r2.data, 0, (int)kPlusOne);
  700. r2.Normalize ();
  701. if (r2 < x) {
  702. Kernel.MinusEq (x, r2);
  703. } else {
  704. BigInteger val = new BigInteger (Sign.Positive, kPlusOne + 1);
  705. val.data [kPlusOne] = 0x00000001;
  706. Kernel.MinusEq (val, r2);
  707. Kernel.PlusEq (x, val);
  708. }
  709. while (x >= n)
  710. Kernel.MinusEq (x, n);
  711. }
  712. public BigInteger Multiply (BigInteger a, BigInteger b)
  713. {
  714. if (a == 0 || b == 0) return 0;
  715. if (a.length >= mod.length << 1)
  716. a %= mod;
  717. if (b.length >= mod.length << 1)
  718. b %= mod;
  719. if (a.length >= mod.length)
  720. BarrettReduction (a);
  721. if (b.length >= mod.length)
  722. BarrettReduction (b);
  723. BigInteger ret = new BigInteger (a * b);
  724. BarrettReduction (ret);
  725. return ret;
  726. }
  727. public BigInteger Difference (BigInteger a, BigInteger b)
  728. {
  729. Sign cmp = Kernel.Compare (a, b);
  730. BigInteger diff;
  731. switch (cmp) {
  732. case Sign.Zero:
  733. return 0;
  734. case Sign.Positive:
  735. diff = a - b; break;
  736. case Sign.Negative:
  737. diff = b - a; break;
  738. default:
  739. throw new Exception ();
  740. }
  741. if (diff >= mod) {
  742. if (diff.length >= mod.length << 1)
  743. diff %= mod;
  744. else
  745. BarrettReduction (diff);
  746. }
  747. if (cmp == Sign.Negative)
  748. diff = mod - diff;
  749. return diff;
  750. }
  751. public BigInteger Pow (BigInteger b, BigInteger exp)
  752. {
  753. if ((mod.data [0] & 1) == 1) return OddPow (b, exp);
  754. else return EvenPow (b, exp);
  755. }
  756. public BigInteger EvenPow (BigInteger b, BigInteger exp)
  757. {
  758. BigInteger resultNum = new BigInteger ((BigInteger)1, mod.length << 1);
  759. BigInteger tempNum = new BigInteger (b % mod, mod.length << 1); // ensures (tempNum * tempNum) < b^ (2k)
  760. uint totalBits = (uint)exp.bitCount ();
  761. uint [] wkspace = new uint [mod.length << 1];
  762. // perform squaring and multiply exponentiation
  763. for (uint pos = 0; pos < totalBits; pos++) {
  764. if (exp.testBit (pos)) {
  765. Array.Clear (wkspace, 0, wkspace.Length);
  766. Kernel.Multiply (resultNum.data, 0, resultNum.length, tempNum.data, 0, tempNum.length, wkspace, 0);
  767. resultNum.length += tempNum.length;
  768. uint [] t = wkspace;
  769. wkspace = resultNum.data;
  770. resultNum.data = t;
  771. BarrettReduction (resultNum);
  772. }
  773. Kernel.SquarePositive (tempNum, ref wkspace);
  774. BarrettReduction (tempNum);
  775. if (tempNum == 1) {
  776. return resultNum;
  777. }
  778. }
  779. return resultNum;
  780. }
  781. private BigInteger OddPow (BigInteger b, BigInteger exp)
  782. {
  783. BigInteger resultNum = new BigInteger (Montgomery.ToMont (1, mod), mod.length << 1);
  784. BigInteger tempNum = new BigInteger (Montgomery.ToMont (b, mod), mod.length << 1); // ensures (tempNum * tempNum) < b^ (2k)
  785. uint mPrime = Montgomery.Inverse (mod.data [0]);
  786. uint totalBits = (uint)exp.bitCount ();
  787. uint [] wkspace = new uint [mod.length << 1];
  788. // perform squaring and multiply exponentiation
  789. for (uint pos = 0; pos < totalBits; pos++) {
  790. if (exp.testBit (pos)) {
  791. Array.Clear (wkspace, 0, wkspace.Length);
  792. Kernel.Multiply (resultNum.data, 0, resultNum.length, tempNum.data, 0, tempNum.length, wkspace, 0);
  793. resultNum.length += tempNum.length;
  794. uint [] t = wkspace;
  795. wkspace = resultNum.data;
  796. resultNum.data = t;
  797. Montgomery.Reduce (resultNum, mod, mPrime);
  798. }
  799. Kernel.SquarePositive (tempNum, ref wkspace);
  800. Montgomery.Reduce (tempNum, mod, mPrime);
  801. }
  802. Montgomery.Reduce (resultNum, mod, mPrime);
  803. return resultNum;
  804. }
  805. #region Pow Small Base
  806. // TODO: Make tests for this, not really needed b/c prime stuff
  807. // checks it, but still would be nice
  808. public BigInteger Pow (uint b, BigInteger exp)
  809. {
  810. if (b != 2) {
  811. if ((mod.data [0] & 1) == 1) return OddPow (b, exp);
  812. else return EvenPow (b, exp);
  813. } else {
  814. if ((mod.data [0] & 1) == 1) return OddModTwoPow (exp);
  815. else return EvenModTwoPow (exp);
  816. }
  817. }
  818. private unsafe BigInteger OddPow (uint b, BigInteger exp)
  819. {
  820. exp.Normalize ();
  821. uint [] wkspace = new uint [mod.length << 1 + 1];
  822. BigInteger resultNum = Montgomery.ToMont ((BigInteger)b, this.mod);
  823. resultNum = new BigInteger (resultNum, mod.length << 1 +1);
  824. uint mPrime = Montgomery.Inverse (mod.data [0]);
  825. uint pos = (uint)exp.bitCount () - 2;
  826. //
  827. // We know that the first itr will make the val b
  828. //
  829. do {
  830. //
  831. // r = r ^ 2 % m
  832. //
  833. Kernel.SquarePositive (resultNum, ref wkspace);
  834. resultNum = Montgomery.Reduce (resultNum, mod, mPrime);
  835. if (exp.testBit (pos)) {
  836. //
  837. // r = r * b % m
  838. //
  839. // TODO: Is Unsafe really speeding things up?
  840. fixed (uint* u = resultNum.data) {
  841. uint i = 0;
  842. ulong mc = 0;
  843. do {
  844. mc += (ulong)u [i] * (ulong)b;
  845. u [i] = (uint)mc;
  846. mc >>= 32;
  847. } while (++i < resultNum.length);
  848. if (resultNum.length < mod.length) {
  849. if (mc != 0) {
  850. u [i] = (uint)mc;
  851. resultNum.length++;
  852. while (resultNum >= mod)
  853. Kernel.MinusEq (resultNum, mod);
  854. }
  855. } else if (mc != 0) {
  856. //
  857. // First, we estimate the quotient by dividing
  858. // the first part of each of the numbers. Then
  859. // we correct this, if necessary, with a subtraction.
  860. //
  861. uint cc = (uint)mc;
  862. // We would rather have this estimate overshoot,
  863. // so we add one to the divisor
  864. uint divEstimate = (uint) ((((ulong)cc << 32) | (ulong) u [i -1]) /
  865. (mod.data [mod.length-1] + 1));
  866. uint t;
  867. i = 0;
  868. mc = 0;
  869. do {
  870. mc += (ulong)mod.data [i] * (ulong)divEstimate;
  871. t = u [i];
  872. u [i] -= (uint)mc;
  873. mc >>= 32;
  874. if (u [i] > t) mc++;
  875. i++;
  876. } while (i < resultNum.length);
  877. cc -= (uint)mc;
  878. if (cc != 0) {
  879. uint sc = 0, j = 0;
  880. uint [] s = mod.data;
  881. do {
  882. uint a = s [j];
  883. if (((a += sc) < sc) | ((u [j] -= a) > ~a)) sc = 1;
  884. else sc = 0;
  885. j++;
  886. } while (j < resultNum.length);
  887. cc -= sc;
  888. }
  889. while (resultNum >= mod)
  890. Kernel.MinusEq (resultNum, mod);
  891. } else {
  892. while (resultNum >= mod)
  893. Kernel.MinusEq (resultNum, mod);
  894. }
  895. }
  896. }
  897. } while (pos-- > 0);
  898. resultNum = Montgomery.Reduce (resultNum, mod, mPrime);
  899. return resultNum;
  900. }
  901. private unsafe BigInteger EvenPow (uint b, BigInteger exp)
  902. {
  903. exp.Normalize ();
  904. uint [] wkspace = new uint [mod.length << 1 + 1];
  905. BigInteger resultNum = new BigInteger ((BigInteger)b, mod.length << 1 + 1);
  906. uint pos = (uint)exp.bitCount () - 2;
  907. //
  908. // We know that the first itr will make the val b
  909. //
  910. do {
  911. //
  912. // r = r ^ 2 % m
  913. //
  914. Kernel.SquarePositive (resultNum, ref wkspace);
  915. if (!(resultNum.length < mod.length))
  916. BarrettReduction (resultNum);
  917. if (exp.testBit (pos)) {
  918. //
  919. // r = r * b % m
  920. //
  921. // TODO: Is Unsafe really speeding things up?
  922. fixed (uint* u = resultNum.data) {
  923. uint i = 0;
  924. ulong mc = 0;
  925. do {
  926. mc += (ulong)u [i] * (ulong)b;
  927. u [i] = (uint)mc;
  928. mc >>= 32;
  929. } while (++i < resultNum.length);
  930. if (resultNum.length < mod.length) {
  931. if (mc != 0) {
  932. u [i] = (uint)mc;
  933. resultNum.length++;
  934. while (resultNum >= mod)
  935. Kernel.MinusEq (resultNum, mod);
  936. }
  937. } else if (mc != 0) {
  938. //
  939. // First, we estimate the quotient by dividing
  940. // the first part of each of the numbers. Then
  941. // we correct this, if necessary, with a subtraction.
  942. //
  943. uint cc = (uint)mc;
  944. // We would rather have this estimate overshoot,
  945. // so we add one to the divisor
  946. uint divEstimate = (uint) ((((ulong)cc << 32) | (ulong) u [i -1]) /
  947. (mod.data [mod.length-1] + 1));
  948. uint t;
  949. i = 0;
  950. mc = 0;
  951. do {
  952. mc += (ulong)mod.data [i] * (ulong)divEstimate;
  953. t = u [i];
  954. u [i] -= (uint)mc;
  955. mc >>= 32;
  956. if (u [i] > t) mc++;
  957. i++;
  958. } while (i < resultNum.length);
  959. cc -= (uint)mc;
  960. if (cc != 0) {
  961. uint sc = 0, j = 0;
  962. uint [] s = mod.data;
  963. do {
  964. uint a = s [j];
  965. if (((a += sc) < sc) | ((u [j] -= a) > ~a)) sc = 1;
  966. else sc = 0;
  967. j++;
  968. } while (j < resultNum.length);
  969. cc -= sc;
  970. }
  971. while (resultNum >= mod)
  972. Kernel.MinusEq (resultNum, mod);
  973. } else {
  974. while (resultNum >= mod)
  975. Kernel.MinusEq (resultNum, mod);
  976. }
  977. }
  978. }
  979. } while (pos-- > 0);
  980. return resultNum;
  981. }
  982. private unsafe BigInteger EvenModTwoPow (BigInteger exp)
  983. {
  984. exp.Normalize ();
  985. uint [] wkspace = new uint [mod.length << 1 + 1];
  986. BigInteger resultNum = new BigInteger (2, mod.length << 1 +1);
  987. uint value = exp.data [exp.length - 1];
  988. uint mask = 0x80000000;
  989. // Find the first bit of the exponent
  990. while ((value & mask) == 0)
  991. mask >>= 1;
  992. //
  993. // We know that the first itr will make the val 2,
  994. // so eat one bit of the exponent
  995. //
  996. mask >>= 1;
  997. uint wPos = exp.length - 1;
  998. do {
  999. value = exp.data [wPos];
  1000. do {
  1001. Kernel.SquarePositive (resultNum, ref wkspace);
  1002. if (resultNum.length >= mod.length)
  1003. BarrettReduction (resultNum);
  1004. if ((value & mask) != 0) {
  1005. //
  1006. // resultNum = (resultNum * 2) % mod
  1007. //
  1008. fixed (uint* u = resultNum.data) {
  1009. //
  1010. // Double
  1011. //
  1012. uint* uu = u;
  1013. uint* uuE = u + resultNum.length;
  1014. uint x, carry = 0;
  1015. while (uu < uuE) {
  1016. x = *uu;
  1017. *uu = (x << 1) | carry;
  1018. carry = x >> (32 - 1);
  1019. uu++;
  1020. }
  1021. // subtraction inlined because we know it is square
  1022. if (carry != 0 || resultNum >= mod) {
  1023. uu = u;
  1024. uint c = 0;
  1025. uint [] s = mod.data;
  1026. uint i = 0;
  1027. do {
  1028. uint a = s [i];
  1029. if (((a += c) < c) | ((* (uu++) -= a) > ~a))
  1030. c = 1;
  1031. else
  1032. c = 0;
  1033. i++;
  1034. } while (uu < uuE);
  1035. }
  1036. }
  1037. }
  1038. } while ((mask >>= 1) > 0);
  1039. mask = 0x80000000;
  1040. } while (wPos-- > 0);
  1041. return resultNum;
  1042. }
  1043. private unsafe BigInteger OddModTwoPow (BigInteger exp)
  1044. {
  1045. uint [] wkspace = new uint [mod.length << 1 + 1];
  1046. BigInteger resultNum = Montgomery.ToMont ((BigInteger)2, this.mod);
  1047. resultNum = new BigInteger (resultNum, mod.length << 1 +1);
  1048. uint mPrime = Montgomery.Inverse (mod.data [0]);
  1049. //
  1050. // TODO: eat small bits, the ones we can do with no modular reduction
  1051. //
  1052. uint pos = (uint)exp.bitCount () - 2;
  1053. do {
  1054. Kernel.SquarePositive (resultNum, ref wkspace);
  1055. resultNum = Montgomery.Reduce (resultNum, mod, mPrime);
  1056. if (exp.testBit (pos)) {
  1057. //
  1058. // resultNum = (resultNum * 2) % mod
  1059. //
  1060. fixed (uint* u = resultNum.data) {
  1061. //
  1062. // Double
  1063. //
  1064. uint* uu = u;
  1065. uint* uuE = u + resultNum.length;
  1066. uint x, carry = 0;
  1067. while (uu < uuE) {
  1068. x = *uu;
  1069. *uu = (x << 1) | carry;
  1070. carry = x >> (32 - 1);
  1071. uu++;
  1072. }
  1073. // subtraction inlined because we know it is square
  1074. if (carry != 0 || resultNum >= mod) {
  1075. fixed (uint* s = mod.data) {
  1076. uu = u;
  1077. uint c = 0;
  1078. uint* ss = s;
  1079. do {
  1080. uint a = *ss++;
  1081. if (((a += c) < c) | ((* (uu++) -= a) > ~a))
  1082. c = 1;
  1083. else
  1084. c = 0;
  1085. } while (uu < uuE);
  1086. }
  1087. }
  1088. }
  1089. }
  1090. } while (pos-- > 0);
  1091. resultNum = Montgomery.Reduce (resultNum, mod, mPrime);
  1092. return resultNum;
  1093. }
  1094. #endregion
  1095. }
  1096. public sealed class Montgomery {
  1097. public static uint Inverse (uint n)
  1098. {
  1099. uint y = n, z;
  1100. while ((z = n * y) != 1)
  1101. y *= 2 - z;
  1102. return (uint)-y;
  1103. }
  1104. public static BigInteger ToMont (BigInteger n, BigInteger m)
  1105. {
  1106. n.Normalize (); m.Normalize ();
  1107. n <<= (int)m.length * 32;
  1108. n %= m;
  1109. return n;
  1110. }
  1111. public static unsafe BigInteger Reduce (BigInteger n, BigInteger m, uint mPrime)
  1112. {
  1113. BigInteger A = n;
  1114. fixed (uint* a = A.data, mm = m.data) {
  1115. for (uint i = 0; i < m.length; i++) {
  1116. // The mod here is taken care of by the CPU,
  1117. // since the multiply will overflow.
  1118. uint u_i = a [0] * mPrime /* % 2^32 */;
  1119. //
  1120. // A += u_i * m;
  1121. // A >>= 32
  1122. //
  1123. // mP = Position in mod
  1124. // aSP = the source of bits from a
  1125. // aDP = destination for bits
  1126. uint* mP = mm, aSP = a, aDP = a;
  1127. ulong c = (ulong)u_i * (ulong)*(mP++) + *(aSP++);
  1128. c >>= 32;
  1129. uint j = 1;
  1130. // Multiply and add
  1131. for (; j < m.length; j++) {
  1132. c += (ulong)u_i * (ulong)*(mP++) + *(aSP++);
  1133. *(aDP++) = (uint)c;
  1134. c >>= 32;
  1135. }
  1136. // Account for carry
  1137. // TODO: use a better loop here, we dont need the ulong stuff
  1138. for (; j < A.length; j++) {
  1139. c += *(aSP++);
  1140. *(aDP++) = (uint)c;
  1141. c >>= 32;
  1142. if (c == 0) {j++; break;}
  1143. }
  1144. // Copy the rest
  1145. for (; j < A.length; j++) {
  1146. *(aDP++) = *(aSP++);
  1147. }
  1148. *(aDP++) = (uint)c;
  1149. }
  1150. while (A.length > 1 && a [A.length-1] == 0) A.length--;
  1151. }
  1152. if (A >= m) Kernel.MinusEq (A, m);
  1153. return A;
  1154. }
  1155. public static BigInteger Reduce (BigInteger n, BigInteger m)
  1156. {
  1157. return Reduce (n, m, Inverse (m.data [0]));
  1158. }
  1159. }
  1160. /// <summary>
  1161. /// Low level functions for the BigInteger
  1162. /// </summary>
  1163. private sealed class Kernel {
  1164. #region Addition/Subtraction
  1165. /// <summary>
  1166. /// Adds two numbers with the same sign.
  1167. /// </summary>
  1168. /// <param name="bi1">A BigInteger</param>
  1169. /// <param name="bi2">A BigInteger</param>
  1170. /// <returns>bi1 + bi2</returns>
  1171. public static BigInteger AddSameSign (BigInteger bi1, BigInteger bi2)
  1172. {
  1173. uint [] x, y;
  1174. uint yMax, xMax, i = 0;
  1175. // x should be bigger
  1176. if (bi1.length < bi2.length) {
  1177. x = bi2.data;
  1178. xMax = bi2.length;
  1179. y = bi1.data;
  1180. yMax = bi1.length;
  1181. } else {
  1182. x = bi1.data;
  1183. xMax = bi1.length;
  1184. y = bi2.data;
  1185. yMax = bi2.length;
  1186. }
  1187. BigInteger result = new BigInteger (Sign.Positive, xMax + 1);
  1188. uint [] r = result.data;
  1189. ulong sum = 0;
  1190. // Add common parts of both numbers
  1191. do {
  1192. sum = ((ulong)x [i]) + ((ulong)y [i]) + sum;
  1193. r [i] = (uint)sum;
  1194. sum >>= 32;
  1195. } while (++i < yMax);
  1196. // Copy remainder of longer number while carry propagation is required
  1197. bool carry = (sum != 0);
  1198. if (carry) {
  1199. if (i < xMax) {
  1200. do
  1201. carry = ((r [i] = x [i] + 1) == 0);
  1202. while (++i < xMax && carry);
  1203. }
  1204. if (carry) {
  1205. r [i] = 1;
  1206. result.length = ++i;
  1207. return result;
  1208. }
  1209. }
  1210. // Copy the rest
  1211. if (i < xMax) {
  1212. do
  1213. r [i] = x [i];
  1214. while (++i < xMax);
  1215. }
  1216. result.Normalize ();
  1217. return result;
  1218. }
  1219. public static BigInteger Subtract (BigInteger big, BigInteger small)
  1220. {
  1221. BigInteger result = new BigInteger (Sign.Positive, big.length);
  1222. uint [] r = result.data, b = big.data, s = small.data;
  1223. uint i = 0, c = 0;
  1224. do {
  1225. uint x = s [i];
  1226. if (((x += c) < c) | ((r [i] = b [i] - x) > ~x))
  1227. c = 1;
  1228. else
  1229. c = 0;
  1230. } while (++i < small.length);
  1231. if (i == big.length) goto fixup;
  1232. if (c == 1) {
  1233. do
  1234. r [i] = b [i] - 1;
  1235. while (b [i++] == 0 && i < big.length);
  1236. if (i == big.length) goto fixup;
  1237. }
  1238. do
  1239. r [i] = b [i];
  1240. while (++i < big.length);
  1241. fixup:
  1242. result.Normalize ();
  1243. return result;
  1244. }
  1245. public static void MinusEq (BigInteger big, BigInteger small)
  1246. {
  1247. uint [] b = big.data, s = small.data;
  1248. uint i = 0, c = 0;
  1249. do {
  1250. uint x = s [i];
  1251. if (((x += c) < c) | ((b [i] -= x) > ~x))
  1252. c = 1;
  1253. else
  1254. c = 0;
  1255. } while (++i < small.length);
  1256. if (i == big.length) goto fixup;
  1257. if (c == 1) {
  1258. do
  1259. b [i]--;
  1260. while (b [i++] == 0 && i < big.length);
  1261. }
  1262. fixup:
  1263. // Normalize length
  1264. while (big.length > 0 && big.data [big.length-1] == 0) big.length--;
  1265. // Check for zero
  1266. if (big.length == 0)
  1267. big.length++;
  1268. }
  1269. public static void PlusEq (BigInteger bi1, BigInteger bi2)
  1270. {
  1271. uint [] x, y;
  1272. uint yMax, xMax, i = 0;
  1273. bool flag = false;
  1274. // x should be bigger
  1275. if (bi1.length < bi2.length){
  1276. flag = true;
  1277. x = bi2.data;
  1278. xMax = bi2.length;
  1279. y = bi1.data;
  1280. yMax = bi1.length;
  1281. } else {
  1282. x = bi1.data;
  1283. xMax = bi1.length;
  1284. y = bi2.data;
  1285. yMax = bi2.length;
  1286. }
  1287. uint [] r = bi1.data;
  1288. ulong sum = 0;
  1289. // Add common parts of both numbers
  1290. do {
  1291. sum += ((ulong)x [i]) + ((ulong)y [i]);
  1292. r [i] = (uint)sum;
  1293. sum >>= 32;
  1294. } while (++i < yMax);
  1295. // Copy remainder of longer number while carry propagation is required
  1296. bool carry = (sum != 0);
  1297. if (carry){
  1298. if (i < xMax) {
  1299. do
  1300. carry = ((r [i] = x [i] + 1) == 0);
  1301. while (++i < xMax && carry);
  1302. }
  1303. if (carry) {
  1304. r [i] = 1;
  1305. bi1.length = ++i;
  1306. return;
  1307. }
  1308. }
  1309. // Copy the rest
  1310. if (flag && i < xMax - 1) {
  1311. do
  1312. r [i] = x [i];
  1313. while (++i < xMax);
  1314. }
  1315. bi1.length = xMax + 1;
  1316. bi1.Normalize ();
  1317. }
  1318. #endregion
  1319. #region Compare
  1320. /// <summary>
  1321. /// Compares two BigInteger
  1322. /// </summary>
  1323. /// <param name="bi1">A BigInteger</param>
  1324. /// <param name="bi2">A BigInteger</param>
  1325. /// <returns>The sign of bi1 - bi2</returns>
  1326. public static Sign Compare (BigInteger bi1, BigInteger bi2)
  1327. {
  1328. //
  1329. // Step 1. Compare the lengths
  1330. //
  1331. uint l1 = bi1.length, l2 = bi2.length;
  1332. while (l1 > 0 && bi1.data [l1-1] == 0) l1--;
  1333. while (l2 > 0 && bi2.data [l2-1] == 0) l2--;
  1334. if (l1 == 0 && l2 == 0) return Sign.Zero;
  1335. // bi1 len < bi2 len
  1336. if (l1 < l2) return Sign.Negative;
  1337. // bi1 len > bi2 len
  1338. else if (l1 > l2) return Sign.Positive;
  1339. //
  1340. // Step 2. Compare the bits
  1341. //
  1342. uint pos = l1 - 1;
  1343. while (pos != 0 && bi1.data [pos] == bi2.data [pos]) pos--;
  1344. if (bi1.data [pos] < bi2.data [pos])
  1345. return Sign.Negative;
  1346. else if (bi1.data [pos] > bi2.data [pos])
  1347. return Sign.Positive;
  1348. else
  1349. return Sign.Zero;
  1350. }
  1351. #endregion
  1352. #region Division
  1353. #region Dword
  1354. /// <summary>
  1355. /// Performs n / d and n % d in one operation.
  1356. /// </summary>
  1357. /// <param name="n">A BigInteger, upon exit this will hold n / d</param>
  1358. /// <param name="d">The divisor</param>
  1359. /// <returns>n % d</returns>
  1360. public static uint SingleByteDivideInPlace (BigInteger n, uint d)
  1361. {
  1362. ulong r = 0;
  1363. uint i = n.length;
  1364. while (i-- > 0) {
  1365. r <<= 32;
  1366. r |= n.data [i];
  1367. n.data [i] = (uint)(r / d);
  1368. r %= d;
  1369. }
  1370. n.Normalize ();
  1371. return (uint)r;
  1372. }
  1373. public static uint DwordMod (BigInteger n, uint d)
  1374. {
  1375. ulong r = 0;
  1376. uint i = n.length;
  1377. while (i-- > 0) {
  1378. r <<= 32;
  1379. r |= n.data [i];
  1380. r %= d;
  1381. }
  1382. return (uint)r;
  1383. }
  1384. public static BigInteger DwordDiv (BigInteger n, uint d)
  1385. {
  1386. BigInteger ret = new BigInteger (Sign.Positive, n.length);
  1387. ulong r = 0;
  1388. uint i = n.length;
  1389. while (i-- > 0) {
  1390. r <<= 32;
  1391. r |= n.data [i];
  1392. ret.data [i] = (uint)(r / d);
  1393. r %= d;
  1394. }
  1395. ret.Normalize ();
  1396. return ret;
  1397. }
  1398. public static BigInteger [] DwordDivMod (BigInteger n, uint d)
  1399. {
  1400. BigInteger ret = new BigInteger (Sign.Positive , n.length);
  1401. ulong r = 0;
  1402. uint i = n.length;
  1403. while (i-- > 0) {
  1404. r <<= 32;
  1405. r |= n.data [i];
  1406. ret.data [i] = (uint)(r / d);
  1407. r %= d;
  1408. }
  1409. ret.Normalize ();
  1410. BigInteger rem = (uint)r;
  1411. return new BigInteger [] {ret, rem};
  1412. }
  1413. #endregion
  1414. #region BigNum
  1415. public static BigInteger [] multiByteDivide (BigInteger bi1, BigInteger bi2)
  1416. {
  1417. if (Kernel.Compare (bi1, bi2) == Sign.Negative)
  1418. return new BigInteger [2] { 0, new BigInteger (bi1) };
  1419. bi1.Normalize (); bi2.Normalize ();
  1420. if (bi2.length == 1)
  1421. return DwordDivMod (bi1, bi2.data [0]);
  1422. uint remainderLen = bi1.length + 1;
  1423. int divisorLen = (int)bi2.length + 1;
  1424. uint mask = 0x80000000;
  1425. uint val = bi2.data [bi2.length - 1];
  1426. int shift = 0;
  1427. int resultPos = (int)bi1.length - (int)bi2.length;
  1428. while (mask != 0 && (val & mask) == 0) {
  1429. shift++; mask >>= 1;
  1430. }
  1431. BigInteger quot = new BigInteger (Sign.Positive, bi1.length - bi2.length + 1);
  1432. BigInteger rem = (bi1 << shift);
  1433. uint [] remainder = rem.data;
  1434. bi2 = bi2 << shift;
  1435. int j = (int)(remainderLen - bi2.length);
  1436. int pos = (int)remainderLen - 1;
  1437. uint firstDivisorByte = bi2.data [bi2.length-1];
  1438. ulong secondDivisorByte = bi2.data [bi2.length-2];
  1439. while (j > 0) {
  1440. ulong dividend = ((ulong)remainder [pos] << 32) + (ulong)remainder [pos-1];
  1441. ulong q_hat = dividend / (ulong)firstDivisorByte;
  1442. ulong r_hat = dividend % (ulong)firstDivisorByte;
  1443. do {
  1444. if (q_hat == 0x100000000 ||
  1445. (q_hat * secondDivisorByte) > ((r_hat << 32) + remainder [pos-2])) {
  1446. q_hat--;
  1447. r_hat += (ulong)firstDivisorByte;
  1448. if (r_hat < 0x100000000)
  1449. continue;
  1450. }
  1451. break;
  1452. } while (true);
  1453. //
  1454. // At this point, q_hat is either exact, or one too large
  1455. // (more likely to be exact) so, we attempt to multiply the
  1456. // divisor by q_hat, if we get a borrow, we just subtract
  1457. // one from q_hat and add the divisor back.
  1458. //
  1459. uint t;
  1460. uint dPos = 0;
  1461. int nPos = pos - divisorLen + 1;
  1462. ulong mc = 0;
  1463. uint uint_q_hat = (uint)q_hat;
  1464. do {
  1465. mc += (ulong)bi2.data [dPos] * (ulong)uint_q_hat;
  1466. t = remainder [nPos];
  1467. remainder [nPos] -= (uint)mc;
  1468. mc >>= 32;
  1469. if (remainder [nPos] > t) mc++;
  1470. dPos++; nPos++;
  1471. } while (dPos < divisorLen);
  1472. nPos = pos - divisorLen + 1;
  1473. dPos = 0;
  1474. // Overestimate
  1475. if (mc != 0) {
  1476. uint_q_hat--;
  1477. ulong sum = 0;
  1478. do {
  1479. sum = ((ulong)remainder [nPos]) + ((ulong)bi2.data [dPos]) + sum;
  1480. remainder [nPos] = (uint)sum;
  1481. sum >>= 32;
  1482. dPos++; nPos++;
  1483. } while (dPos < divisorLen);
  1484. }
  1485. quot.data [resultPos--] = (uint)uint_q_hat;
  1486. pos--;
  1487. j--;
  1488. }
  1489. quot.Normalize ();
  1490. rem.Normalize ();
  1491. BigInteger [] ret = new BigInteger [2] { quot, rem };
  1492. if (shift != 0)
  1493. ret [1] >>= shift;
  1494. return ret;
  1495. }
  1496. #endregion
  1497. #endregion
  1498. #region Shift
  1499. public static BigInteger LeftShift (BigInteger bi, int n)
  1500. {
  1501. if (n == 0) return new BigInteger (bi, bi.length + 1);
  1502. int w = n >> 5;
  1503. n &= ((1 << 5) - 1);
  1504. BigInteger ret = new BigInteger (Sign.Positive, bi.length + 1 + (uint)w);
  1505. uint i = 0, l = bi.length;
  1506. if (n != 0) {
  1507. uint x, carry = 0;
  1508. while (i < l) {
  1509. x = bi.data [i];
  1510. ret.data [i + w] = (x << n) | carry;
  1511. carry = x >> (32 - n);
  1512. i++;
  1513. }
  1514. ret.data [i + w] = carry;
  1515. } else {
  1516. while (i < l) {
  1517. ret.data [i + w] = bi.data [i];
  1518. i++;
  1519. }
  1520. }
  1521. ret.Normalize ();
  1522. return ret;
  1523. }
  1524. public static BigInteger RightShift (BigInteger bi, int n)
  1525. {
  1526. if (n == 0) return new BigInteger (bi);
  1527. int w = n >> 5;
  1528. int s = n & ((1 << 5) - 1);
  1529. BigInteger ret = new BigInteger (Sign.Positive, bi.length - (uint)w + 1);
  1530. uint l = (uint)ret.data.Length - 1;
  1531. if (s != 0) {
  1532. uint x, carry = 0;
  1533. while (l-- > 0) {
  1534. x = bi.data [l + w];
  1535. ret.data [l] = (x >> n) | carry;
  1536. carry = x << (32 - n);
  1537. }
  1538. } else {
  1539. while (l-- > 0)
  1540. ret.data [l] = bi.data [l + w];
  1541. }
  1542. ret.Normalize ();
  1543. return ret;
  1544. }
  1545. #endregion
  1546. #region Multiply
  1547. public static BigInteger MultiplyByDword (BigInteger n, uint f)
  1548. {
  1549. BigInteger ret = new BigInteger (Sign.Positive, n.length + 1);
  1550. uint i = 0;
  1551. ulong c = 0;
  1552. do {
  1553. c += (ulong)n.data [i] * (ulong)f;
  1554. ret.data [i] = (uint)c;
  1555. c >>= 32;
  1556. } while (++i < n.length);
  1557. ret.data [i] = (uint)c;
  1558. ret.Normalize ();
  1559. return ret;
  1560. }
  1561. /// <summary>
  1562. /// Multiplies the data in x [xOffset:xOffset+xLen] by
  1563. /// y [yOffset:yOffset+yLen] and puts it into
  1564. /// d [dOffset:dOffset+xLen+yLen].
  1565. /// </summary>
  1566. /// <remarks>
  1567. /// This code is unsafe! It is the caller's responsibility to make
  1568. /// sure that it is safe to access x [xOffset:xOffset+xLen],
  1569. /// y [yOffset:yOffset+yLen], and d [dOffset:dOffset+xLen+yLen].
  1570. /// </remarks>
  1571. public static unsafe void Multiply (uint [] x, uint xOffset, uint xLen, uint [] y, uint yOffset, uint yLen, uint [] d, uint dOffset)
  1572. {
  1573. fixed (uint* x

Large files files are truncated, but you can click here to view the full file