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