mcs/class/corlib/Mono.Math/BigInteger.cs

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