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