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, 2007 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                         if (inData.Length == 0)
 212                                 inData = new byte [1];
 213                         length = (uint)inData.Length >> 2;
 214                         int leftOver = inData.Length & 0x3;
 215
 216                         // length not multiples of 4
 217                         if (leftOver != 0) length++;
 218
 219                         data = new uint [length];
 220
 221                         for (int i = inData.Length - 1, j = 0; i >= 3; i -= 4, j++) {
 222                                 data [j] = (uint)(
 223                                         (inData [i-3] << (3*8)) |
 224                                         (inData [i-2] << (2*8)) |
 225                                         (inData [i-1] << (1*8)) |
 226                                         (inData [i])
 227                                         );
 228                         }
 229
 230                         switch (leftOver) {
 231                         case 1: data [length-1] = (uint)inData [0]; break;
 232                         case 2: data [length-1] = (uint)((inData [0] << 8) | inData [1]); break;
 233                         case 3: data [length-1] = (uint)((inData [0] << 16) | (inData [1] << 8) | inData [2]); break;
 234                         }
 235
 236                         this.Normalize ();
 237                 }
 238
 239 #if !INSIDE_CORLIB
 240                 [CLSCompliant (false)]
 241 #endif
 242                 public BigInteger (uint [] inData)
 243                 {
 244                         if (inData.Length == 0)
 245                                 inData = new uint [1];
 246                         length = (uint)inData.Length;
 247
 248                         data = new uint [length];
 249
 250                         for (int i = (int)length - 1, j = 0; i >= 0; i--, j++)
 251                                 data [j] = inData [i];
 252
 253                         this.Normalize ();
 254                 }
 255
 256 #if !INSIDE_CORLIB
 257                 [CLSCompliant (false)]
 258 #endif
 259                 public BigInteger (uint ui)
 260                 {
 261                         data = new uint [] {ui};
 262                 }
 263
 264 #if !INSIDE_CORLIB
 265                 [CLSCompliant (false)]
 266 #endif
 267                 public BigInteger (ulong ul)
 268                 {
 269                         data = new uint [2] { (uint)ul, (uint)(ul >> 32)};
 270                         length = 2;
 271
 272                         this.Normalize ();
 273                 }
 274
 275 #if !INSIDE_CORLIB
 276                 [CLSCompliant (false)]
 277 #endif
 278                 public static implicit operator BigInteger (uint value)
 279                 {
 280                         return (new BigInteger (value));
 281                 }
 282
 283                 public static implicit operator BigInteger (int value)
 284                 {
 285                         if (value < 0) throw new ArgumentOutOfRangeException ("value");
 286                         return (new BigInteger ((uint)value));
 287                 }
 288
 289 #if !INSIDE_CORLIB
 290                 [CLSCompliant (false)]
 291 #endif
 292                 public static implicit operator BigInteger (ulong value)
 293                 {
 294                         return (new BigInteger (value));
 295                 }
 296
 297                 /* This is the BigInteger.Parse method I use. This method works
 298                 because BigInteger.ToString returns the input I gave to Parse. */
 299                 public static BigInteger Parse (string number)
 300                 {
 301                         if (number == null)
 302                                 throw new ArgumentNullException ("number");
 303
 304                         int i = 0, len = number.Length;
 305                         char c;
 306                         bool digits_seen = false;
 307                         BigInteger val = new BigInteger (0);
 308                         if (number [i] == '+') {
 309                                 i++;
 310                         }
 311                         else if (number [i] == '-') {
 312                                 throw new FormatException (WouldReturnNegVal);
 313                         }
 314
 315                         for (; i < len; i++) {
 316                                 c = number [i];
 317                                 if (c == '\0') {
 318                                         i = len;
 319                                         continue;
 320                                 }
 321                                 if (c >= '0' && c <= '9') {
 322                                         val = val * 10 + (c - '0');
 323                                         digits_seen = true;
 324                                 }
 325                                 else {
 326                                         if (Char.IsWhiteSpace (c)) {
 327                                                 for (i++; i < len; i++) {
 328                                                         if (!Char.IsWhiteSpace (number [i]))
 329                                                                 throw new FormatException ();
 330                                                 }
 331                                                 break;
 332                                         }
 333                                         else
 334                                                 throw new FormatException ();
 335                                 }
 336                         }
 337                         if (!digits_seen)
 338                                 throw new FormatException ();
 339                         return val;
 340                 }
 341
 342                 #endregion
 343
 344                 #region Operators
 345
 346                 public static BigInteger operator + (BigInteger bi1, BigInteger bi2)
 347                 {
 348                         if (bi1 == 0)
 349                                 return new BigInteger (bi2);
 350                         else if (bi2 == 0)
 351                                 return new BigInteger (bi1);
 352                         else
 353                                 return Kernel.AddSameSign (bi1, bi2);
 354                 }
 355
 356                 public static BigInteger operator - (BigInteger bi1, BigInteger bi2)
 357                 {
 358                         if (bi2 == 0)
 359                                 return new BigInteger (bi1);
 360
 361                         if (bi1 == 0)
 362                                 throw new ArithmeticException (WouldReturnNegVal);
 363
 364                         switch (Kernel.Compare (bi1, bi2)) {
 365
 366                                 case Sign.Zero:
 367                                         return 0;
 368
 369                                 case Sign.Positive:
 370                                         return Kernel.Subtract (bi1, bi2);
 371
 372                                 case Sign.Negative:
 373                                         throw new ArithmeticException (WouldReturnNegVal);
 374                                 default:
 375                                         throw new Exception ();
 376                         }
 377                 }
 378
 379                 public static int operator % (BigInteger bi, int i)
 380                 {
 381                         if (i > 0)
 382                                 return (int)Kernel.DwordMod (bi, (uint)i);
 383                         else
 384                                 return -(int)Kernel.DwordMod (bi, (uint)-i);
 385                 }
 386
 387 #if !INSIDE_CORLIB
 388                 [CLSCompliant (false)]
 389 #endif
 390                 public static uint operator % (BigInteger bi, uint ui)
 391                 {
 392                         return Kernel.DwordMod (bi, (uint)ui);
 393                 }
 394
 395                 public static BigInteger operator % (BigInteger bi1, BigInteger bi2)
 396                 {
 397                         return Kernel.multiByteDivide (bi1, bi2)[1];
 398                 }
 399
 400                 public static BigInteger operator / (BigInteger bi, int i)
 401                 {
 402                         if (i > 0)
 403                                 return Kernel.DwordDiv (bi, (uint)i);
 404
 405                         throw new ArithmeticException (WouldReturnNegVal);
 406                 }
 407
 408                 public static BigInteger operator / (BigInteger bi1, BigInteger bi2)
 409                 {
 410                         return Kernel.multiByteDivide (bi1, bi2)[0];
 411                 }
 412
 413                 public static BigInteger operator * (BigInteger bi1, BigInteger bi2)
 414                 {
 415                         if (bi1 == 0 || bi2 == 0) return 0;
 416
 417                         //
 418                         // Validate pointers
 419                         //
 420                         if (bi1.data.Length < bi1.length) throw new IndexOutOfRangeException ("bi1 out of range");
 421                         if (bi2.data.Length < bi2.length) throw new IndexOutOfRangeException ("bi2 out of range");
 422
 423                         BigInteger ret = new BigInteger (Sign.Positive, bi1.length + bi2.length);
 424
 425                         Kernel.Multiply (bi1.data, 0, bi1.length, bi2.data, 0, bi2.length, ret.data, 0);
 426
 427                         ret.Normalize ();
 428                         return ret;
 429                 }
 430
 431                 public static BigInteger operator * (BigInteger bi, int i)
 432                 {
 433                         if (i < 0) throw new ArithmeticException (WouldReturnNegVal);
 434                         if (i == 0) return 0;
 435                         if (i == 1) return new BigInteger (bi);
 436
 437                         return Kernel.MultiplyByDword (bi, (uint)i);
 438                 }
 439
 440                 public static BigInteger operator << (BigInteger bi1, int shiftVal)
 441                 {
 442                         return Kernel.LeftShift (bi1, shiftVal);
 443                 }
 444
 445                 public static BigInteger operator >> (BigInteger bi1, int shiftVal)
 446                 {
 447                         return Kernel.RightShift (bi1, shiftVal);
 448                 }
 449
 450                 #endregion
 451
 452                 #region Friendly names for operators
 453
 454                 // with names suggested by FxCop 1.30
 455
 456                 public static BigInteger Add (BigInteger bi1, BigInteger bi2)
 457                 {
 458                         return (bi1 + bi2);
 459                 }
 460
 461                 public static BigInteger Subtract (BigInteger bi1, BigInteger bi2)
 462                 {
 463                         return (bi1 - bi2);
 464                 }
 465
 466                 public static int Modulus (BigInteger bi, int i)
 467                 {
 468                         return (bi % i);
 469                 }
 470
 471 #if !INSIDE_CORLIB
 472                 [CLSCompliant (false)]
 473 #endif
 474                 public static uint Modulus (BigInteger bi, uint ui)
 475                 {
 476                         return (bi % ui);
 477                 }
 478
 479                 public static BigInteger Modulus (BigInteger bi1, BigInteger bi2)
 480                 {
 481                         return (bi1 % bi2);
 482                 }
 483
 484                 public static BigInteger Divid (BigInteger bi, int i)
 485                 {
 486                         return (bi / i);
 487                 }
 488
 489                 public static BigInteger Divid (BigInteger bi1, BigInteger bi2)
 490                 {
 491                         return (bi1 / bi2);
 492                 }
 493
 494                 public static BigInteger Multiply (BigInteger bi1, BigInteger bi2)
 495                 {
 496                         return (bi1 * bi2);
 497                 }
 498
 499                 public static BigInteger Multiply (BigInteger bi, int i)
 500                 {
 501                         return (bi * i);
 502                 }
 503
 504                 #endregion
 505
 506                 #region Random
 507                 private static RandomNumberGenerator rng;
 508                 private static RandomNumberGenerator Rng {
 509                         get {
 510                                 if (rng == null)
 511                                         rng = RandomNumberGenerator.Create ();
 512                                 return rng;
 513                         }
 514                 }
 515
 516                 /// <summary>
 517                 /// Generates a new, random BigInteger of the specified length.
 518                 /// </summary>
 519                 /// <param name="bits">The number of bits for the new number.</param>
 520                 /// <param name="rng">A random number generator to use to obtain the bits.</param>
 521                 /// <returns>A random number of the specified length.</returns>
 522                 public static BigInteger GenerateRandom (int bits, RandomNumberGenerator rng)
 523                 {
 524                         int dwords = bits >> 5;
 525                         int remBits = bits & 0x1F;
 526
 527                         if (remBits != 0)
 528                                 dwords++;
 529
 530                         BigInteger ret = new BigInteger (Sign.Positive, (uint)dwords + 1);
 531                         byte [] random = new byte [dwords << 2];
 532
 533                         rng.GetBytes (random);
 534                         Buffer.BlockCopy (random, 0, ret.data, 0, (int)dwords << 2);
 535
 536                         if (remBits != 0) {
 537                                 uint mask = (uint)(0x01 << (remBits-1));
 538                                 ret.data [dwords-1] |= mask;
 539
 540                                 mask = (uint)(0xFFFFFFFF >> (32 - remBits));
 541                                 ret.data [dwords-1] &= mask;
 542                         }
 543                         else
 544                                 ret.data [dwords-1] |= 0x80000000;
 545
 546                         ret.Normalize ();
 547                         return ret;
 548                 }
 549
 550                 /// <summary>
 551                 /// Generates a new, random BigInteger of the specified length using the default RNG crypto service provider.
 552                 /// </summary>
 553                 /// <param name="bits">The number of bits for the new number.</param>
 554                 /// <returns>A random number of the specified length.</returns>
 555                 public static BigInteger GenerateRandom (int bits)
 556                 {
 557                         return GenerateRandom (bits, Rng);
 558                 }
 559
 560                 /// <summary>
 561                 /// Randomizes the bits in "this" from the specified RNG.
 562                 /// </summary>
 563                 /// <param name="rng">A RNG.</param>
 564                 public void Randomize (RandomNumberGenerator rng)
 565                 {
 566                         if (this == 0)
 567                                 return;
 568
 569                         int bits = this.BitCount ();
 570                         int dwords = bits >> 5;
 571                         int remBits = bits & 0x1F;
 572
 573                         if (remBits != 0)
 574                                 dwords++;
 575
 576                         byte [] random = new byte [dwords << 2];
 577
 578                         rng.GetBytes (random);
 579                         Buffer.BlockCopy (random, 0, data, 0, (int)dwords << 2);
 580
 581                         if (remBits != 0) {
 582                                 uint mask = (uint)(0x01 << (remBits-1));
 583                                 data [dwords-1] |= mask;
 584
 585                                 mask = (uint)(0xFFFFFFFF >> (32 - remBits));
 586                                 data [dwords-1] &= mask;
 587                         }
 588
 589                         else
 590                                 data [dwords-1] |= 0x80000000;
 591
 592                         Normalize ();
 593                 }
 594
 595                 /// <summary>
 596                 /// Randomizes the bits in "this" from the default RNG.
 597                 /// </summary>
 598                 public void Randomize ()
 599                 {
 600                         Randomize (Rng);
 601                 }
 602
 603                 #endregion
 604
 605                 #region Bitwise
 606
 607                 public int BitCount ()
 608                 {
 609                         this.Normalize ();
 610
 611                         uint value = data [length - 1];
 612                         uint mask = 0x80000000;
 613                         uint bits = 32;
 614
 615                         while (bits > 0 && (value & mask) == 0) {
 616                                 bits--;
 617                                 mask >>= 1;
 618                         }
 619                         bits += ((length - 1) << 5);
 620
 621                         return (int)bits;
 622                 }
 623
 624                 /// <summary>
 625                 /// Tests if the specified bit is 1.
 626                 /// </summary>
 627                 /// <param name="bitNum">The bit to test. The least significant bit is 0.</param>
 628                 /// <returns>True if bitNum is set to 1, else false.</returns>
 629 #if !INSIDE_CORLIB
 630                 [CLSCompliant (false)]
 631 #endif
 632                 public bool TestBit (uint bitNum)
 633                 {
 634                         uint bytePos = bitNum >> 5;             // divide by 32
 635                         byte bitPos = (byte)(bitNum & 0x1F);    // get the lowest 5 bits
 636
 637                         uint mask = (uint)1 << bitPos;
 638                         return ((this.data [bytePos] & mask) != 0);
 639                 }
 640
 641                 public bool TestBit (int bitNum)
 642                 {
 643                         if (bitNum < 0) throw new IndexOutOfRangeException ("bitNum out of range");
 644
 645                         uint bytePos = (uint)bitNum >> 5;             // divide by 32
 646                         byte bitPos = (byte)(bitNum & 0x1F);    // get the lowest 5 bits
 647
 648                         uint mask = (uint)1 << bitPos;
 649                         return ((this.data [bytePos] | mask) == this.data [bytePos]);
 650                 }
 651
 652 #if !INSIDE_CORLIB
 653                 [CLSCompliant (false)]
 654 #endif
 655                 public void SetBit (uint bitNum)
 656                 {
 657                         SetBit (bitNum, true);
 658                 }
 659
 660 #if !INSIDE_CORLIB
 661                 [CLSCompliant (false)]
 662 #endif
 663                 public void ClearBit (uint bitNum)
 664                 {
 665                         SetBit (bitNum, false);
 666                 }
 667
 668 #if !INSIDE_CORLIB
 669                 [CLSCompliant (false)]
 670 #endif
 671                 public void SetBit (uint bitNum, bool value)
 672                 {
 673                         uint bytePos = bitNum >> 5;             // divide by 32
 674
 675                         if (bytePos < this.length) {
 676                                 uint mask = (uint)1 << (int)(bitNum & 0x1F);
 677                                 if (value)
 678                                         this.data [bytePos] |= mask;
 679                                 else
 680                                         this.data [bytePos] &= ~mask;
 681                         }
 682                 }
 683
 684                 public int LowestSetBit ()
 685                 {
 686                         if (this == 0) return -1;
 687                         int i = 0;
 688                         while (!TestBit (i)) i++;
 689                         return i;
 690                 }
 691
 692                 public byte[] GetBytes ()
 693                 {
 694                         if (this == 0) return new byte [1];
 695
 696                         int numBits = BitCount ();
 697                         int numBytes = numBits >> 3;
 698                         if ((numBits & 0x7) != 0)
 699                                 numBytes++;
 700
 701                         byte [] result = new byte [numBytes];
 702
 703                         int numBytesInWord = numBytes & 0x3;
 704                         if (numBytesInWord == 0) numBytesInWord = 4;
 705
 706                         int pos = 0;
 707                         for (int i = (int)length - 1; i >= 0; i--) {
 708                                 uint val = data [i];
 709                                 for (int j = numBytesInWord - 1; j >= 0; j--) {
 710                                         result [pos+j] = (byte)(val & 0xFF);
 711                                         val >>= 8;
 712                                 }
 713                                 pos += numBytesInWord;
 714                                 numBytesInWord = 4;
 715                         }
 716                         return result;
 717                 }
 718
 719                 #endregion
 720
 721                 #region Compare
 722
 723 #if !INSIDE_CORLIB
 724                 [CLSCompliant (false)]
 725 #endif
 726                 public static bool operator == (BigInteger bi1, uint ui)
 727                 {
 728                         if (bi1.length != 1) bi1.Normalize ();
 729                         return bi1.length == 1 && bi1.data [0] == ui;
 730                 }
 731
 732 #if !INSIDE_CORLIB
 733                 [CLSCompliant (false)]
 734 #endif
 735                 public static bool operator != (BigInteger bi1, uint ui)
 736                 {
 737                         if (bi1.length != 1) bi1.Normalize ();
 738                         return !(bi1.length == 1 && bi1.data [0] == ui);
 739                 }
 740
 741                 public static bool operator == (BigInteger bi1, BigInteger bi2)
 742                 {
 743                         // we need to compare with null
 744                         if ((bi1 as object) == (bi2 as object))
 745                                 return true;
 746                         if (null == bi1 || null == bi2)
 747                                 return false;
 748                         return Kernel.Compare (bi1, bi2) == 0;
 749                 }
 750
 751                 public static bool operator != (BigInteger bi1, BigInteger bi2)
 752                 {
 753                         // we need to compare with null
 754                         if ((bi1 as object) == (bi2 as object))
 755                                 return false;
 756                         if (null == bi1 || null == bi2)
 757                                 return true;
 758                         return Kernel.Compare (bi1, bi2) != 0;
 759                 }
 760
 761                 public static bool operator > (BigInteger bi1, BigInteger bi2)
 762                 {
 763                         return Kernel.Compare (bi1, bi2) > 0;
 764                 }
 765
 766                 public static bool operator < (BigInteger bi1, BigInteger bi2)
 767                 {
 768                         return Kernel.Compare (bi1, bi2) < 0;
 769                 }
 770
 771                 public static bool operator >= (BigInteger bi1, BigInteger bi2)
 772                 {
 773                         return Kernel.Compare (bi1, bi2) >= 0;
 774                 }
 775
 776                 public static bool operator <= (BigInteger bi1, BigInteger bi2)
 777                 {
 778                         return Kernel.Compare (bi1, bi2) <= 0;
 779                 }
 780
 781                 public Sign Compare (BigInteger bi)
 782                 {
 783                         return Kernel.Compare (this, bi);
 784                 }
 785
 786                 #endregion
 787
 788                 #region Formatting
 789
 790 #if !INSIDE_CORLIB
 791                 [CLSCompliant (false)]
 792 #endif
 793                 public string ToString (uint radix)
 794                 {
 795                         return ToString (radix, "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ");
 796                 }
 797
 798 #if !INSIDE_CORLIB
 799                 [CLSCompliant (false)]
 800 #endif
 801                 public string ToString (uint radix, string characterSet)
 802                 {
 803                         if (characterSet.Length < radix)
 804                                 throw new ArgumentException ("charSet length less than radix", "characterSet");
 805                         if (radix == 1)
 806                                 throw new ArgumentException ("There is no such thing as radix one notation", "radix");
 807
 808                         if (this == 0) return "0";
 809                         if (this == 1) return "1";
 810
 811                         string result = "";
 812
 813                         BigInteger a = new BigInteger (this);
 814
 815                         while (a != 0) {
 816                                 uint rem = Kernel.SingleByteDivideInPlace (a, radix);
 817                                 result = characterSet [(int) rem] + result;
 818                         }
 819
 820                         return result;
 821                 }
 822
 823                 #endregion
 824
 825                 #region Misc
 826
 827                 /// <summary>
 828                 ///     Normalizes this by setting the length to the actual number of
 829                 ///     uints used in data and by setting the sign to Sign.Zero if the
 830                 ///     value of this is 0.
 831                 /// </summary>
 832                 private void Normalize ()
 833                 {
 834                         // Normalize length
 835                         while (length > 0 && data [length-1] == 0) length--;
 836
 837                         // Check for zero
 838                         if (length == 0)
 839                                 length++;
 840                 }
 841
 842                 public void Clear ()
 843                 {
 844                         for (int i=0; i < length; i++)
 845                                 data [i] = 0x00;
 846                 }
 847
 848                 #endregion
 849
 850                 #region Object Impl
 851
 852                 public override int GetHashCode ()
 853                 {
 854                         uint val = 0;
 855
 856                         for (uint i = 0; i < this.length; i++)
 857                                 val ^= this.data [i];
 858
 859                         return (int)val;
 860                 }
 861
 862                 public override string ToString ()
 863                 {
 864                         return ToString (10);
 865                 }
 866
 867                 public override bool Equals (object o)
 868                 {
 869                         if (o == null)
 870                                 return false;
 871                         if (o is int)
 872                                 return (int)o >= 0 && this == (uint)o;
 873
 874                         BigInteger bi = o as BigInteger;
 875                         if (bi == null)
 876                                 return false;
 877
 878                         return Kernel.Compare (this, bi) == 0;
 879                 }
 880
 881                 #endregion
 882
 883                 #region Number Theory
 884
 885                 public BigInteger GCD (BigInteger bi)
 886                 {
 887                         return Kernel.gcd (this, bi);
 888                 }
 889
 890                 public BigInteger ModInverse (BigInteger modulus)
 891                 {
 892                         return Kernel.modInverse (this, modulus);
 893                 }
 894
 895                 public BigInteger ModPow (BigInteger exp, BigInteger n)
 896                 {
 897                         ModulusRing mr = new ModulusRing (n);
 898                         return mr.Pow (this, exp);
 899                 }
 900
 901                 #endregion
 902
 903                 #region Prime Testing
 904
 905                 public bool IsProbablePrime ()
 906                 {
 907                         // can we use our small-prime table ?
 908                         if (this <= smallPrimes[smallPrimes.Length - 1]) {
 909                                 for (int p = 0; p < smallPrimes.Length; p++) {
 910                                         if (this == smallPrimes[p])
 911                                                 return true;
 912                                 }
 913                                 // the list is complete, so it's not a prime
 914                                 return false;
 915                         }
 916
 917                         // otherwise check if we can divide by one of the small primes
 918                         for (int p = 0; p < smallPrimes.Length; p++) {
 919                                 if (this % smallPrimes[p] == 0)
 920                                         return false;
 921                         }
 922                         // the last step is to confirm the "large" prime with the SPP or Miller-Rabin test
 923                         return PrimalityTests.Test (this, Prime.ConfidenceFactor.Medium);
 924                 }
 925
 926                 #endregion
 927
 928                 #region Prime Number Generation
 929
 930                 /// <summary>
 931                 /// Generates the smallest prime >= bi
 932                 /// </summary>
 933                 /// <param name="bi">A BigInteger</param>
 934                 /// <returns>The smallest prime >= bi. More mathematically, if bi is prime: bi, else Prime [PrimePi [bi] + 1].</returns>
 935                 public static BigInteger NextHighestPrime (BigInteger bi)
 936                 {
 937                         NextPrimeFinder npf = new NextPrimeFinder ();
 938                         return npf.GenerateNewPrime (0, bi);
 939                 }
 940
 941                 public static BigInteger GeneratePseudoPrime (int bits)
 942                 {
 943                         SequentialSearchPrimeGeneratorBase sspg = new SequentialSearchPrimeGeneratorBase ();
 944                         return sspg.GenerateNewPrime (bits);
 945                 }
 946
 947                 /// <summary>
 948                 /// Increments this by two
 949                 /// </summary>
 950                 public void Incr2 ()
 951                 {
 952                         int i = 0;
 953
 954                         data [0] += 2;
 955
 956                         // If there was no carry, nothing to do
 957                         if (data [0] < 2) {
 958
 959                                 // Account for the first carry
 960                                 data [++i]++;
 961
 962                                 // Keep adding until no carry
 963                                 while (data [i++] == 0x0)
 964                                         data [i]++;
 965
 966                                 // See if we increased the data length
 967                                 if (length == (uint)i)
 968                                         length++;
 969                         }
 970                 }
 971
 972                 #endregion
 973
 974 #if INSIDE_CORLIB
 975                 internal
 976 #else
 977                 public
 978 #endif
 979                 sealed class ModulusRing {
 980
 981                         BigInteger mod, constant;
 982
 983                         public ModulusRing (BigInteger modulus)
 984                         {
 985                                 this.mod = modulus;
 986
 987                                 // calculate constant = b^ (2k) / m
 988                                 uint i = mod.length << 1;
 989
 990                                 constant = new BigInteger (Sign.Positive, i + 1);
 991                                 constant.data [i] = 0x00000001;
 992
 993                                 constant = constant / mod;
 994                         }
 995
 996                         public void BarrettReduction (BigInteger x)
 997                         {
 998                                 BigInteger n = mod;
 999                                 uint k = n.length,
1000                                         kPlusOne = k+1,
1001                                         kMinusOne = k-1;
1002
1003                                 // x < mod, so nothing to do.
1004                                 if (x.length < k) return;
1005
1006                                 BigInteger q3;
1007
1008                                 //
1009                                 // Validate pointers
1010                                 //
1011                                 if (x.data.Length < x.length) throw new IndexOutOfRangeException ("x out of range");
1012
1013                                 // q1 = x / b^ (k-1)
1014                                 // q2 = q1 * constant
1015                                 // q3 = q2 / b^ (k+1), Needs to be accessed with an offset of kPlusOne
1016
1017                                 // TODO: We should the method in HAC p 604 to do this (14.45)
1018                                 q3 = new BigInteger (Sign.Positive, x.length - kMinusOne + constant.length);
1019                                 Kernel.Multiply (x.data, kMinusOne, x.length - kMinusOne, constant.data, 0, constant.length, q3.data, 0);
1020
1021                                 // r1 = x mod b^ (k+1)
1022                                 // i.e. keep the lowest (k+1) words
1023
1024                                 uint lengthToCopy = (x.length > kPlusOne) ? kPlusOne : x.length;
1025
1026                                 x.length = lengthToCopy;
1027                                 x.Normalize ();
1028
1029                                 // r2 = (q3 * n) mod b^ (k+1)
1030                                 // partial multiplication of q3 and n
1031
1032                                 BigInteger r2 = new BigInteger (Sign.Positive, kPlusOne);
1033                                 Kernel.MultiplyMod2p32pmod (q3.data, (int)kPlusOne, (int)q3.length - (int)kPlusOne, n.data, 0, (int)n.length, r2.data, 0, (int)kPlusOne);
1034
1035                                 r2.Normalize ();
1036
1037                                 if (r2 <= x) {
1038                                         Kernel.MinusEq (x, r2);
1039                                 } else {
1040                                         BigInteger val = new BigInteger (Sign.Positive, kPlusOne + 1);
1041                                         val.data [kPlusOne] = 0x00000001;
1042
1043                                         Kernel.MinusEq (val, r2);
1044                                         Kernel.PlusEq (x, val);
1045                                 }
1046
1047                                 while (x >= n)
1048                                         Kernel.MinusEq (x, n);
1049                         }
1050
1051                         public BigInteger Multiply (BigInteger a, BigInteger b)
1052                         {
1053                                 if (a == 0 || b == 0) return 0;
1054
1055                                 if (a > mod)
1056                                         a %= mod;
1057
1058                                 if (b > mod)
1059                                         b %= mod;
1060
1061                                 BigInteger ret = a * b;
1062                                 BarrettReduction (ret);
1063
1064                                 return ret;
1065                         }
1066
1067                         public BigInteger Difference (BigInteger a, BigInteger b)
1068                         {
1069                                 Sign cmp = Kernel.Compare (a, b);
1070                                 BigInteger diff;
1071
1072                                 switch (cmp) {
1073                                         case Sign.Zero:
1074                                                 return 0;
1075                                         case Sign.Positive:
1076                                                 diff = a - b; break;
1077                                         case Sign.Negative:
1078                                                 diff = b - a; break;
1079                                         default:
1080                                                 throw new Exception ();
1081                                 }
1082
1083                                 if (diff >= mod) {
1084                                         if (diff.length >= mod.length << 1)
1085                                                 diff %= mod;
1086                                         else
1087                                                 BarrettReduction (diff);
1088                                 }
1089                                 if (cmp == Sign.Negative)
1090                                         diff = mod - diff;
1091                                 return diff;
1092                         }
1093 #if true
1094                         public BigInteger Pow (BigInteger a, BigInteger k)
1095                         {
1096                                 BigInteger b = new BigInteger (1);
1097                                 if (k == 0)
1098                                         return b;
1099
1100                                 BigInteger A = a;
1101                                 if (k.TestBit (0))
1102                                         b = a;
1103
1104                                 for (int i = 1; i < k.BitCount (); i++) {
1105                                         A = Multiply (A, A);
1106                                         if (k.TestBit (i))
1107                                                 b = Multiply (A, b);
1108                                 }
1109                                 return b;
1110                         }
1111 #else
1112                         public BigInteger Pow (BigInteger b, BigInteger exp)
1113                         {
1114                                 if ((mod.data [0] & 1) == 1) return OddPow (b, exp);
1115                                 else return EvenPow (b, exp);
1116                         }
1117
1118                         public BigInteger EvenPow (BigInteger b, BigInteger exp)
1119                         {
1120                                 BigInteger resultNum = new BigInteger ((BigInteger)1, mod.length << 1);
1121                                 BigInteger tempNum = new BigInteger (b % mod, mod.length << 1);  // ensures (tempNum * tempNum) < b^ (2k)
1122
1123                                 uint totalBits = (uint)exp.BitCount ();
1124
1125                                 uint [] wkspace = new uint [mod.length << 1];
1126
1127                                 // perform squaring and multiply exponentiation
1128                                 for (uint pos = 0; pos < totalBits; pos++) {
1129                                         if (exp.TestBit (pos)) {
1130
1131                                                 Array.Clear (wkspace, 0, wkspace.Length);
1132                                                 Kernel.Multiply (resultNum.data, 0, resultNum.length, tempNum.data, 0, tempNum.length, wkspace, 0);
1133                                                 resultNum.length += tempNum.length;
1134                                                 uint [] t = wkspace;
1135                                                 wkspace = resultNum.data;
1136                                                 resultNum.data = t;
1137
1138                                                 BarrettReduction (resultNum);
1139                                         }
1140
1141                                         Kernel.SquarePositive (tempNum, ref wkspace);
1142                                         BarrettReduction (tempNum);
1143
1144                                         if (tempNum == 1) {
1145                                                 return resultNum;
1146                                         }
1147                                 }
1148
1149                                 return resultNum;
1150                         }
1151
1152                         private BigInteger OddPow (BigInteger b, BigInteger exp)
1153                         {
1154                                 BigInteger resultNum = new BigInteger (Montgomery.ToMont (1, mod), mod.length << 1);
1155                                 BigInteger tempNum = new BigInteger (Montgomery.ToMont (b, mod), mod.length << 1);  // ensures (tempNum * tempNum) < b^ (2k)
1156                                 uint mPrime = Montgomery.Inverse (mod.data [0]);
1157                                 uint totalBits = (uint)exp.BitCount ();
1158
1159                                 uint [] wkspace = new uint [mod.length << 1];
1160
1161                                 // perform squaring and multiply exponentiation
1162                                 for (uint pos = 0; pos < totalBits; pos++) {
1163                                         if (exp.TestBit (pos)) {
1164
1165                                                 Array.Clear (wkspace, 0, wkspace.Length);
1166                                                 Kernel.Multiply (resultNum.data, 0, resultNum.length, tempNum.data, 0, tempNum.length, wkspace, 0);
1167                                                 resultNum.length += tempNum.length;
1168                                                 uint [] t = wkspace;
1169                                                 wkspace = resultNum.data;
1170                                                 resultNum.data = t;
1171
1172                                                 Montgomery.Reduce (resultNum, mod, mPrime);
1173                                         }
1174
1175                                         // the value of tempNum is required in the last loop
1176                                         if (pos < totalBits - 1) {
1177                                                 Kernel.SquarePositive (tempNum, ref wkspace);
1178                                                 Montgomery.Reduce (tempNum, mod, mPrime);
1179                                         }
1180                                 }
1181
1182                                 Montgomery.Reduce (resultNum, mod, mPrime);
1183                                 return resultNum;
1184                         }
1185 #endif
1186                         #region Pow Small Base
1187
1188                         // TODO: Make tests for this, not really needed b/c prime stuff
1189                         // checks it, but still would be nice
1190 #if !INSIDE_CORLIB
1191                         [CLSCompliant (false)]
1192 #endif
1193 #if true
1194                         public BigInteger Pow (uint b, BigInteger exp)
1195                         {
1196                                 return Pow (new BigInteger (b), exp);
1197                         }
1198 #else
1199                         public BigInteger Pow (uint b, BigInteger exp)
1200                         {
1201 //                              if (b != 2) {
1202                                         if ((mod.data [0] & 1) == 1)
1203                                                 return OddPow (b, exp);
1204                                         else
1205                                                 return EvenPow (b, exp);
1206 /* buggy in some cases (like the well tested primes)
1207                                 } else {
1208                                         if ((mod.data [0] & 1) == 1)
1209                                                 return OddModTwoPow (exp);
1210                                         else
1211                                                 return EvenModTwoPow (exp);
1212                                 }*/
1213                         }
1214
1215                         private unsafe BigInteger OddPow (uint b, BigInteger exp)
1216                         {
1217                                 exp.Normalize ();
1218                                 uint [] wkspace = new uint [mod.length << 1 + 1];
1219
1220                                 BigInteger resultNum = Montgomery.ToMont ((BigInteger)b, this.mod);
1221                                 resultNum = new BigInteger (resultNum, mod.length << 1 +1);
1222
1223                                 uint mPrime = Montgomery.Inverse (mod.data [0]);
1224
1225                                 int bc = exp.BitCount () - 2;
1226                                 uint pos = (bc > 1 ? (uint) bc : 1);
1227
1228                                 //
1229                                 // We know that the first itr will make the val b
1230                                 //
1231
1232                                 do {
1233                                         //
1234                                         // r = r ^ 2 % m
1235                                         //
1236                                         Kernel.SquarePositive (resultNum, ref wkspace);
1237                                         resultNum = Montgomery.Reduce (resultNum, mod, mPrime);
1238
1239                                         if (exp.TestBit (pos)) {
1240
1241                                                 //
1242                                                 // r = r * b % m
1243                                                 //
1244
1245                                                 // TODO: Is Unsafe really speeding things up?
1246                                                 fixed (uint* u = resultNum.data) {
1247
1248                                                         uint i = 0;
1249                                                         ulong mc = 0;
1250
1251                                                         do {
1252                                                                 mc += (ulong)u [i] * (ulong)b;
1253                                                                 u [i] = (uint)mc;
1254                                                                 mc >>= 32;
1255                                                         } while (++i < resultNum.length);
1256
1257                                                         if (resultNum.length < mod.length) {
1258                                                                 if (mc != 0) {
1259                                                                         u [i] = (uint)mc;
1260                                                                         resultNum.length++;
1261                                                                         while (resultNum >= mod)
1262                                                                                 Kernel.MinusEq (resultNum, mod);
1263                                                                 }
1264                                                         } else if (mc != 0) {
1265
1266                                                                 //
1267                                                                 // First, we estimate the quotient by dividing
1268                                                                 // the first part of each of the numbers. Then
1269                                                                 // we correct this, if necessary, with a subtraction.
1270                                                                 //
1271
1272                                                                 uint cc = (uint)mc;
1273
1274                                                                 // We would rather have this estimate overshoot,
1275                                                                 // so we add one to the divisor
1276                                                                 uint divEstimate;
1277                                                                 if (mod.data [mod.length - 1] < UInt32.MaxValue) {
1278                                                                         divEstimate = (uint) ((((ulong)cc << 32) | (ulong) u [i -1]) /
1279                                                                                 (mod.data [mod.length-1] + 1));
1280                                                                 }
1281                                                                 else {
1282                                                                         // guess but don't divide by 0
1283                                                                         divEstimate = (uint) ((((ulong)cc << 32) | (ulong) u [i -1]) /
1284                                                                                 (mod.data [mod.length-1]));
1285                                                                 }
1286
1287                                                                 uint t;
1288
1289                                                                 i = 0;
1290                                                                 mc = 0;
1291                                                                 do {
1292                                                                         mc += (ulong)mod.data [i] * (ulong)divEstimate;
1293                                                                         t = u [i];
1294                                                                         u [i] -= (uint)mc;
1295                                                                         mc >>= 32;
1296                                                                         if (u [i] > t) mc++;
1297                                                                         i++;
1298                                                                 } while (i < resultNum.length);
1299                                                                 cc -= (uint)mc;
1300
1301                                                                 if (cc != 0) {
1302
1303                                                                         uint sc = 0, j = 0;
1304                                                                         uint [] s = mod.data;
1305                                                                         do {
1306                                                                                 uint a = s [j];
1307                                                                                 if (((a += sc) < sc) | ((u [j] -= a) > ~a)) sc = 1;
1308                                                                                 else sc = 0;
1309                                                                                 j++;
1310                                                                         } while (j < resultNum.length);
1311                                                                         cc -= sc;
1312                                                                 }
1313                                                                 while (resultNum >= mod)
1314                                                                         Kernel.MinusEq (resultNum, mod);
1315                                                         } else {
1316                                                                 while (resultNum >= mod)
1317                                                                         Kernel.MinusEq (resultNum, mod);
1318                                                         }
1319                                                 }
1320                                         }
1321                                 } while (pos-- > 0);
1322
1323                                 resultNum = Montgomery.Reduce (resultNum, mod, mPrime);
1324                                 return resultNum;
1325
1326                         }
1327
1328                         private unsafe BigInteger EvenPow (uint b, BigInteger exp)
1329                         {
1330                                 exp.Normalize ();
1331                                 uint [] wkspace = new uint [mod.length << 1 + 1];
1332                                 BigInteger resultNum = new BigInteger ((BigInteger)b, mod.length << 1 + 1);
1333
1334                                 uint pos = (uint)exp.BitCount () - 2;
1335
1336                                 //
1337                                 // We know that the first itr will make the val b
1338                                 //
1339
1340                                 do {
1341                                         //
1342                                         // r = r ^ 2 % m
1343                                         //
1344                                         Kernel.SquarePositive (resultNum, ref wkspace);
1345                                         if (!(resultNum.length < mod.length))
1346                                                 BarrettReduction (resultNum);
1347
1348                                         if (exp.TestBit (pos)) {
1349
1350                                                 //
1351                                                 // r = r * b % m
1352                                                 //
1353
1354                                                 // TODO: Is Unsafe really speeding things up?
1355                                                 fixed (uint* u = resultNum.data) {
1356
1357                                                         uint i = 0;
1358                                                         ulong mc = 0;
1359
1360                                                         do {
1361                                                                 mc += (ulong)u [i] * (ulong)b;
1362                                                                 u [i] = (uint)mc;
1363                                                                 mc >>= 32;
1364                                                         } while (++i < resultNum.length);
1365
1366                                                         if (resultNum.length < mod.length) {
1367                                                                 if (mc != 0) {
1368                                                                         u [i] = (uint)mc;
1369                                                                         resultNum.length++;
1370                                                                         while (resultNum >= mod)
1371                                                                                 Kernel.MinusEq (resultNum, mod);
1372                                                                 }
1373                                                         } else if (mc != 0) {
1374
1375                                                                 //
1376                                                                 // First, we estimate the quotient by dividing
1377                                                                 // the first part of each of the numbers. Then
1378                                                                 // we correct this, if necessary, with a subtraction.
1379                                                                 //
1380
1381                                                                 uint cc = (uint)mc;
1382
1383                                                                 // We would rather have this estimate overshoot,
1384                                                                 // so we add one to the divisor
1385                                                                 uint divEstimate = (uint) ((((ulong)cc << 32) | (ulong) u [i -1]) /
1386                                                                         (mod.data [mod.length-1] + 1));
1387
1388                                                                 uint t;
1389
1390                                                                 i = 0;
1391                                                                 mc = 0;
1392                                                                 do {
1393                                                                         mc += (ulong)mod.data [i] * (ulong)divEstimate;
1394                                                                         t = u [i];
1395                                                                         u [i] -= (uint)mc;
1396                                                                         mc >>= 32;
1397                                                                         if (u [i] > t) mc++;
1398                                                                         i++;
1399                                                                 } while (i < resultNum.length);
1400                                                                 cc -= (uint)mc;
1401
1402                                                                 if (cc != 0) {
1403
1404                                                                         uint sc = 0, j = 0;
1405                                                                         uint [] s = mod.data;
1406                                                                         do {
1407                                                                                 uint a = s [j];
1408                                                                                 if (((a += sc) < sc) | ((u [j] -= a) > ~a)) sc = 1;
1409                                                                                 else sc = 0;
1410                                                                                 j++;
1411                                                                         } while (j < resultNum.length);
1412                                                                         cc -= sc;
1413                                                                 }
1414                                                                 while (resultNum >= mod)
1415                                                                         Kernel.MinusEq (resultNum, mod);
1416                                                         } else {
1417                                                                 while (resultNum >= mod)
1418                                                                         Kernel.MinusEq (resultNum, mod);
1419                                                         }
1420                                                 }
1421                                         }
1422                                 } while (pos-- > 0);
1423
1424                                 return resultNum;
1425                         }
1426 #endif
1427 /* known to be buggy in some cases */
1428 #if false
1429                         private unsafe BigInteger EvenModTwoPow (BigInteger exp)
1430                         {
1431                                 exp.Normalize ();
1432                                 uint [] wkspace = new uint [mod.length << 1 + 1];
1433
1434                                 BigInteger resultNum = new BigInteger (2, mod.length << 1 +1);
1435
1436                                 uint value = exp.data [exp.length - 1];
1437                                 uint mask = 0x80000000;
1438
1439                                 // Find the first bit of the exponent
1440                                 while ((value & mask) == 0)
1441                                         mask >>= 1;
1442
1443                                 //
1444                                 // We know that the first itr will make the val 2,
1445                                 // so eat one bit of the exponent
1446                                 //
1447                                 mask >>= 1;
1448
1449                                 uint wPos = exp.length - 1;
1450
1451                                 do {
1452                                         value = exp.data [wPos];
1453                                         do {
1454                                                 Kernel.SquarePositive (resultNum, ref wkspace);
1455                                                 if (resultNum.length >= mod.length)
1456                                                         BarrettReduction (resultNum);
1457
1458                                                 if ((value & mask) != 0) {
1459                                                         //
1460                                                         // resultNum = (resultNum * 2) % mod
1461                                                         //
1462
1463                                                         fixed (uint* u = resultNum.data) {
1464                                                                 //
1465                                                                 // Double
1466                                                                 //
1467                                                                 uint* uu = u;
1468                                                                 uint* uuE = u + resultNum.length;
1469                                                                 uint x, carry = 0;
1470                                                                 while (uu < uuE) {
1471                                                                         x = *uu;
1472                                                                         *uu = (x << 1) | carry;
1473                                                                         carry = x >> (32 - 1);
1474                                                                         uu++;
1475                                                                 }
1476
1477                                                                 // subtraction inlined because we know it is square
1478                                                                 if (carry != 0 || resultNum >= mod) {
1479                                                                         uu = u;
1480                                                                         uint c = 0;
1481                                                                         uint [] s = mod.data;
1482                                                                         uint i = 0;
1483                                                                         do {
1484                                                                                 uint a = s [i];
1485                                                                                 if (((a += c) < c) | ((* (uu++) -= a) > ~a))
1486                                                                                         c = 1;
1487                                                                                 else
1488                                                                                         c = 0;
1489                                                                                 i++;
1490                                                                         } while (uu < uuE);
1491                                                                 }
1492                                                         }
1493                                                 }
1494                                         } while ((mask >>= 1) > 0);
1495                                         mask = 0x80000000;
1496                                 } while (wPos-- > 0);
1497
1498                                 return resultNum;
1499                         }
1500
1501                         private unsafe BigInteger OddModTwoPow (BigInteger exp)
1502                         {
1503
1504                                 uint [] wkspace = new uint [mod.length << 1 + 1];
1505
1506                                 BigInteger resultNum = Montgomery.ToMont ((BigInteger)2, this.mod);
1507                                 resultNum = new BigInteger (resultNum, mod.length << 1 +1);
1508
1509                                 uint mPrime = Montgomery.Inverse (mod.data [0]);
1510
1511                                 //
1512                                 // TODO: eat small bits, the ones we can do with no modular reduction
1513                                 //
1514                                 uint pos = (uint)exp.BitCount () - 2;
1515
1516                                 do {
1517                                         Kernel.SquarePositive (resultNum, ref wkspace);
1518                                         resultNum = Montgomery.Reduce (resultNum, mod, mPrime);
1519
1520                                         if (exp.TestBit (pos)) {
1521                                                 //
1522                                                 // resultNum = (resultNum * 2) % mod
1523                                                 //
1524
1525                                                 fixed (uint* u = resultNum.data) {
1526                                                         //
1527                                                         // Double
1528                                                         //
1529                                                         uint* uu = u;
1530                                                         uint* uuE = u + resultNum.length;
1531                                                         uint x, carry = 0;
1532                                                         while (uu < uuE) {
1533                                                                 x = *uu;
1534                                                                 *uu = (x << 1) | carry;
1535                                                                 carry = x >> (32 - 1);
1536                                                                 uu++;
1537                                                         }
1538
1539                                                         // subtraction inlined because we know it is square
1540                                                         if (carry != 0 || resultNum >= mod) {
1541                                                                 fixed (uint* s = mod.data) {
1542                                                                         uu = u;
1543                                                                         uint c = 0;
1544                                                                         uint* ss = s;
1545                                                                         do {
1546                                                                                 uint a = *ss++;
1547                                                                                 if (((a += c) < c) | ((* (uu++) -= a) > ~a))
1548                                                                                         c = 1;
1549                                                                                 else
1550                                                                                         c = 0;
1551                                                                         } while (uu < uuE);
1552                                                                 }
1553                                                         }
1554                                                 }
1555                                         }
1556                                 } while (pos-- > 0);
1557
1558                                 resultNum = Montgomery.Reduce (resultNum, mod, mPrime);
1559                                 return resultNum;
1560                         }
1561 #endif
1562                         #endregion
1563                 }
1564
1565                 /// <summary>
1566                 /// Low level functions for the BigInteger
1567                 /// </summary>
1568                 private sealed class Kernel {
1569
1570                         #region Addition/Subtraction
1571
1572                         /// <summary>
1573                         /// Adds two numbers with the same sign.
1574                         /// </summary>
1575                         /// <param name="bi1">A BigInteger</param>
1576                         /// <param name="bi2">A BigInteger</param>
1577                         /// <returns>bi1 + bi2</returns>
1578                         public static BigInteger AddSameSign (BigInteger bi1, BigInteger bi2)
1579                         {
1580                                 uint [] x, y;
1581                                 uint yMax, xMax, i = 0;
1582
1583                                 // x should be bigger
1584                                 if (bi1.length < bi2.length) {
1585                                         x = bi2.data;
1586                                         xMax = bi2.length;
1587                                         y = bi1.data;
1588                                         yMax = bi1.length;
1589                                 } else {
1590                                         x = bi1.data;
1591                                         xMax = bi1.length;
1592                                         y = bi2.data;
1593                                         yMax = bi2.length;
1594                                 }
1595
1596                                 BigInteger result = new BigInteger (Sign.Positive, xMax + 1);
1597
1598                                 uint [] r = result.data;
1599
1600                                 ulong sum = 0;
1601
1602                                 // Add common parts of both numbers
1603                                 do {
1604                                         sum = ((ulong)x [i]) + ((ulong)y [i]) + sum;
1605                                         r [i] = (uint)sum;
1606                                         sum >>= 32;
1607                                 } while (++i < yMax);
1608
1609                                 // Copy remainder of longer number while carry propagation is required
1610                                 bool carry = (sum != 0);
1611
1612                                 if (carry) {
1613
1614                                         if (i < xMax) {
1615                                                 do
1616                                                         carry = ((r [i] = x [i] + 1) == 0);
1617                                                 while (++i < xMax && carry);
1618                                         }
1619
1620                                         if (carry) {
1621                                                 r [i] = 1;
1622                                                 result.length = ++i;
1623                                                 return result;
1624                                         }
1625                                 }
1626
1627                                 // Copy the rest
1628                                 if (i < xMax) {
1629                                         do
1630                                                 r [i] = x [i];
1631                                         while (++i < xMax);
1632                                 }
1633
1634                                 result.Normalize ();
1635                                 return result;
1636                         }
1637
1638                         public static BigInteger Subtract (BigInteger big, BigInteger small)
1639                         {
1640                                 BigInteger result = new BigInteger (Sign.Positive, big.length);
1641
1642                                 uint [] r = result.data, b = big.data, s = small.data;
1643                                 uint i = 0, c = 0;
1644
1645                                 do {
1646
1647                                         uint x = s [i];
1648                                         if (((x += c) < c) | ((r [i] = b [i] - x) > ~x))
1649                                                 c = 1;
1650                                         else
1651                                                 c = 0;
1652
1653                                 } while (++i < small.length);
1654
1655                                 if (i == big.length) goto fixup;
1656
1657                                 if (c == 1) {
1658                                         do
1659                                                 r [i] = b [i] - 1;
1660                                         while (b [i++] == 0 && i < big.length);
1661
1662                                         if (i == big.length) goto fixup;
1663                                 }
1664
1665                                 do
1666                                         r [i] = b [i];
1667                                 while (++i < big.length);
1668
1669                                 fixup:
1670
1671                                         result.Normalize ();
1672                                 return result;
1673                         }
1674
1675                         public static void MinusEq (BigInteger big, BigInteger small)
1676                         {
1677                                 uint [] b = big.data, s = small.data;
1678                                 uint i = 0, c = 0;
1679
1680                                 do {
1681                                         uint x = s [i];
1682                                         if (((x += c) < c) | ((b [i] -= x) > ~x))
1683                                                 c = 1;
1684                                         else
1685                                                 c = 0;
1686                                 } while (++i < small.length);
1687
1688                                 if (i == big.length) goto fixup;
1689
1690                                 if (c == 1) {
1691                                         do
1692                                                 b [i]--;
1693                                         while (b [i++] == 0 && i < big.length);
1694                                 }
1695
1696                                 fixup:
1697
1698                                         // Normalize length
1699                                         while (big.length > 0 && big.data [big.length-1] == 0) big.length--;
1700
1701                                 // Check for zero
1702                                 if (big.length == 0)
1703                                         big.length++;
1704
1705                         }
1706
1707                         public static void PlusEq (BigInteger bi1, BigInteger bi2)
1708                         {
1709                                 uint [] x, y;
1710                                 uint yMax, xMax, i = 0;
1711                                 bool flag = false;
1712
1713                                 // x should be bigger
1714                                 if (bi1.length < bi2.length){
1715                                         flag = true;
1716                                         x = bi2.data;
1717                                         xMax = bi2.length;
1718                                         y = bi1.data;
1719                                         yMax = bi1.length;
1720                                 } else {
1721                                         x = bi1.data;
1722                                         xMax = bi1.length;
1723                                         y = bi2.data;
1724                                         yMax = bi2.length;
1725                                 }
1726
1727                                 uint [] r = bi1.data;
1728
1729                                 ulong sum = 0;
1730
1731                                 // Add common parts of both numbers
1732                                 do {
1733                                         sum += ((ulong)x [i]) + ((ulong)y [i]);
1734                                         r [i] = (uint)sum;
1735                                         sum >>= 32;
1736                                 } while (++i < yMax);
1737
1738                                 // Copy remainder of longer number while carry propagation is required
1739                                 bool carry = (sum != 0);
1740
1741                                 if (carry){
1742
1743                                         if (i < xMax) {
1744                                                 do
1745                                                         carry = ((r [i] = x [i] + 1) == 0);
1746                                                 while (++i < xMax && carry);
1747                                         }
1748
1749                                         if (carry) {
1750                                                 r [i] = 1;
1751                                                 bi1.length = ++i;
1752                                                 return;
1753                                         }
1754                                 }
1755
1756                                 // Copy the rest
1757                                 if (flag && i < xMax - 1) {
1758                                         do
1759                                                 r [i] = x [i];
1760                                         while (++i < xMax);
1761                                 }
1762
1763                                 bi1.length = xMax + 1;
1764                                 bi1.Normalize ();
1765                         }
1766
1767                         #endregion
1768
1769                         #region Compare
1770
1771                         /// <summary>
1772                         /// Compares two BigInteger
1773                         /// </summary>
1774                         /// <param name="bi1">A BigInteger</param>
1775                         /// <param name="bi2">A BigInteger</param>
1776                         /// <returns>The sign of bi1 - bi2</returns>
1777                         public static Sign Compare (BigInteger bi1, BigInteger bi2)
1778                         {
1779                                 //
1780                                 // Step 1. Compare the lengths
1781                                 //
1782                                 uint l1 = bi1.length, l2 = bi2.length;
1783
1784                                 while (l1 > 0 && bi1.data [l1-1] == 0) l1--;
1785                                 while (l2 > 0 && bi2.data [l2-1] == 0) l2--;
1786
1787                                 if (l1 == 0 && l2 == 0) return Sign.Zero;
1788
1789                                 // bi1 len < bi2 len
1790                                 if (l1 < l2) return Sign.Negative;
1791                                 // bi1 len > bi2 len
1792                                 else if (l1 > l2) return Sign.Positive;
1793
1794                                 //
1795                                 // Step 2. Compare the bits
1796                                 //
1797
1798                                 uint pos = l1 - 1;
1799
1800                                 while (pos != 0 && bi1.data [pos] == bi2.data [pos]) pos--;
1801
1802                                 if (bi1.data [pos] < bi2.data [pos])
1803                                         return Sign.Negative;
1804                                 else if (bi1.data [pos] > bi2.data [pos])
1805                                         return Sign.Positive;
1806                                 else
1807                                         return Sign.Zero;
1808                         }
1809
1810                         #endregion
1811
1812                         #region Division
1813
1814                         #region Dword
1815
1816                         /// <summary>
1817                         /// Performs n / d and n % d in one operation.
1818                         /// </summary>
1819                         /// <param name="n">A BigInteger, upon exit this will hold n / d</param>
1820                         /// <param name="d">The divisor</param>
1821                         /// <returns>n % d</returns>
1822                         public static uint SingleByteDivideInPlace (BigInteger n, uint d)
1823                         {
1824                                 ulong r = 0;
1825                                 uint i = n.length;
1826
1827                                 while (i-- > 0) {
1828                                         r <<= 32;
1829                                         r |= n.data [i];
1830                                         n.data [i] = (uint)(r / d);
1831                                         r %= d;
1832                                 }
1833                                 n.Normalize ();
1834
1835                                 return (uint)r;
1836                         }
1837
1838                         public static uint DwordMod (BigInteger n, uint d)
1839                         {
1840                                 ulong r = 0;
1841                                 uint i = n.length;
1842
1843                                 while (i-- > 0) {
1844                                         r <<= 32;
1845                                         r |= n.data [i];
1846                                         r %= d;
1847                                 }
1848
1849                                 return (uint)r;
1850                         }
1851
1852                         public static BigInteger DwordDiv (BigInteger n, uint d)
1853                         {
1854                                 BigInteger ret = new BigInteger (Sign.Positive, n.length);
1855
1856                                 ulong r = 0;
1857                                 uint i = n.length;
1858
1859                                 while (i-- > 0) {
1860                                         r <<= 32;
1861                                         r |= n.data [i];
1862                                         ret.data [i] = (uint)(r / d);
1863                                         r %= d;
1864                                 }
1865                                 ret.Normalize ();
1866
1867                                 return ret;
1868                         }
1869
1870                         public static BigInteger [] DwordDivMod (BigInteger n, uint d)
1871                         {
1872                                 BigInteger ret = new BigInteger (Sign.Positive , n.length);
1873
1874                                 ulong r = 0;
1875                                 uint i = n.length;
1876
1877                                 while (i-- > 0) {
1878                                         r <<= 32;
1879                                         r |= n.data [i];
1880                                         ret.data [i] = (uint)(r / d);
1881                                         r %= d;
1882                                 }
1883                                 ret.Normalize ();
1884
1885                                 BigInteger rem = (uint)r;
1886
1887                                 return new BigInteger [] {ret, rem};
1888                         }
1889
1890                                 #endregion
1891
1892                         #region BigNum
1893
1894                         public static BigInteger [] multiByteDivide (BigInteger bi1, BigInteger bi2)
1895                         {
1896                                 if (Kernel.Compare (bi1, bi2) == Sign.Negative)
1897                                         return new BigInteger [2] { 0, new BigInteger (bi1) };
1898
1899                                 bi1.Normalize (); bi2.Normalize ();
1900
1901                                 if (bi2.length == 1)
1902                                         return DwordDivMod (bi1, bi2.data [0]);
1903
1904                                 uint remainderLen = bi1.length + 1;
1905                                 int divisorLen = (int)bi2.length + 1;
1906
1907                                 uint mask = 0x80000000;
1908                                 uint val = bi2.data [bi2.length - 1];
1909                                 int shift = 0;
1910                                 int resultPos = (int)bi1.length - (int)bi2.length;
1911
1912                                 while (mask != 0 && (val & mask) == 0) {
1913                                         shift++; mask >>= 1;
1914                                 }
1915
1916                                 BigInteger quot = new BigInteger (Sign.Positive, bi1.length - bi2.length + 1);
1917                                 BigInteger rem = (bi1 << shift);
1918
1919                                 uint [] remainder = rem.data;
1920
1921                                 bi2 = bi2 << shift;
1922
1923                                 int j = (int)(remainderLen - bi2.length);
1924                                 int pos = (int)remainderLen - 1;
1925
1926                                 uint firstDivisorByte = bi2.data [bi2.length-1];
1927                                 ulong secondDivisorByte = bi2.data [bi2.length-2];
1928
1929                                 while (j > 0) {
1930                                         ulong dividend = ((ulong)remainder [pos] << 32) + (ulong)remainder [pos-1];
1931
1932                                         ulong q_hat = dividend / (ulong)firstDivisorByte;
1933                                         ulong r_hat = dividend % (ulong)firstDivisorByte;
1934
1935                                         do {
1936
1937                                                 if (q_hat == 0x100000000 ||
1938                                                         (q_hat * secondDivisorByte) > ((r_hat << 32) + remainder [pos-2])) {
1939                                                         q_hat--;
1940                                                         r_hat += (ulong)firstDivisorByte;
1941
1942                                                         if (r_hat < 0x100000000)
1943                                                                 continue;
1944                                                 }
1945                                                 break;
1946                                         } while (true);
1947
1948                                         //
1949                                         // At this point, q_hat is either exact, or one too large
1950                                         // (more likely to be exact) so, we attempt to multiply the
1951                                         // divisor by q_hat, if we get a borrow, we just subtract
1952                                         // one from q_hat and add the divisor back.
1953                                         //
1954
1955                                         uint t;
1956                                         uint dPos = 0;
1957                                         int nPos = pos - divisorLen + 1;
1958                                         ulong mc = 0;
1959                                         uint uint_q_hat = (uint)q_hat;
1960                                         do {
1961                                                 mc += (ulong)bi2.data [dPos] * (ulong)uint_q_hat;
1962                                                 t = remainder [nPos];
1963                                                 remainder [nPos] -= (uint)mc;
1964                                                 mc >>= 32;
1965                                                 if (remainder [nPos] > t) mc++;
1966                                                 dPos++; nPos++;
1967                                         } while (dPos < divisorLen);
1968
1969                                         nPos = pos - divisorLen + 1;
1970                                         dPos = 0;
1971
1972                                         // Overestimate
1973                                         if (mc != 0) {
1974                                                 uint_q_hat--;
1975                                                 ulong sum = 0;
1976
1977                                                 do {
1978                                                         sum = ((ulong)remainder [nPos]) + ((ulong)bi2.data [dPos]) + sum;
1979                                                         remainder [nPos] = (uint)sum;
1980                                                         sum >>= 32;
1981                                                         dPos++; nPos++;
1982                                                 } while (dPos < divisorLen);
1983
1984                                         }
1985
1986                                         quot.data [resultPos--] = (uint)uint_q_hat;
1987
1988                                         pos--;
1989                                         j--;
1990                                 }
1991
1992                                 quot.Normalize ();
1993                                 rem.Normalize ();
1994                                 BigInteger [] ret = new BigInteger [2] { quot, rem };
1995
1996                                 if (shift != 0)
1997                                         ret [1] >>= shift;
1998
1999                                 return ret;
2000                         }
2001
2002                         #endregion
2003
2004                         #endregion
2005
2006                         #region Shift
2007                         public static BigInteger LeftShift (BigInteger bi, int n)
2008                         {
2009                                 if (n == 0) return new BigInteger (bi, bi.length + 1);
2010
2011                                 int w = n >> 5;
2012                                 n &= ((1 << 5) - 1);
2013
2014                                 BigInteger ret = new BigInteger (Sign.Positive, bi.length + 1 + (uint)w);
2015
2016                                 uint i = 0, l = bi.length;
2017                                 if (n != 0) {
2018                                         uint x, carry = 0;
2019                                         while (i < l) {
2020                                                 x = bi.data [i];
2021                                                 ret.data [i + w] = (x << n) | carry;
2022                                                 carry = x >> (32 - n);
2023                                                 i++;
2024                                         }
2025                                         ret.data [i + w] = carry;
2026                                 } else {
2027                                         while (i < l) {
2028                                                 ret.data [i + w] = bi.data [i];
2029                                                 i++;
2030                                         }
2031                                 }
2032
2033                                 ret.Normalize ();
2034                                 return ret;
2035                         }
2036
2037                         public static BigInteger RightShift (BigInteger bi, int n)
2038                         {
2039                                 if (n == 0) return new BigInteger (bi);
2040
2041                                 int w = n >> 5;
2042                                 int s = n & ((1 << 5) - 1);
2043
2044                                 BigInteger ret = new BigInteger (Sign.Positive, bi.length - (uint)w + 1);
2045                                 uint l = (uint)ret.data.Length - 1;
2046
2047                                 if (s != 0) {
2048
2049                                         uint x, carry = 0;
2050
2051                                         while (l-- > 0) {
2052                                                 x = bi.data [l + w];
2053                                                 ret.data [l] = (x >> n) | carry;
2054                                                 carry = x << (32 - n);
2055                                         }
2056                                 } else {
2057                                         while (l-- > 0)
2058                                                 ret.data [l] = bi.data [l + w];
2059
2060                                 }
2061                                 ret.Normalize ();
2062                                 return ret;
2063                         }
2064
2065                         #endregion
2066
2067                         #region Multiply
2068
2069                         public static BigInteger MultiplyByDword (BigInteger n, uint f)
2070                         {
2071                                 BigInteger ret = new BigInteger (Sign.Positive, n.length + 1);
2072
2073                                 uint i = 0;
2074                                 ulong c = 0;
2075
2076                                 do {
2077                                         c += (ulong)n.data [i] * (ulong)f;
2078                                         ret.data [i] = (uint)c;
2079                                         c >>= 32;
2080                                 } while (++i < n.length);
2081                                 ret.data [i] = (uint)c;
2082                                 ret.Normalize ();
2083                                 return ret;
2084
2085                         }
2086
2087                         /// <summary>
2088                         /// Multiplies the data in x [xOffset:xOffset+xLen] by
2089                         /// y [yOffset:yOffset+yLen] and puts it into
2090                         /// d [dOffset:dOffset+xLen+yLen].
2091                         /// </summary>
2092                         /// <remarks>
2093                         /// This code is unsafe! It is the caller's responsibility to make
2094                         /// sure that it is safe to access x [xOffset:xOffset+xLen],
2095                         /// y [yOffset:yOffset+yLen], and d [dOffset:dOffset+xLen+yLen].
2096                         /// </remarks>
2097                         public static unsafe void Multiply (uint [] x, uint xOffset, uint xLen, uint [] y, uint yOffset, uint yLen, uint [] d, uint dOffset)
2098                         {
2099                                 fixed (uint* xx = x, yy = y, dd = d) {
2100                                         uint* xP = xx + xOffset,
2101                                                 xE = xP + xLen,
2102                                                 yB = yy + yOffset,
2103                                                 yE = yB + yLen,
2104                                                 dB = dd + dOffset;
2105
2106                                         for (; xP < xE; xP++, dB++) {
2107
2108                                                 if (*xP == 0) continue;
2109
2110                                                 ulong mcarry = 0;
2111
2112                                                 uint* dP = dB;
2113                                                 for (uint* yP = yB; yP < yE; yP++, dP++) {
2114                                                         mcarry += ((ulong)*xP * (ulong)*yP) + (ulong)*dP;
2115
2116                                                         *dP = (uint)mcarry;
2117                                                         mcarry >>= 32;
2118                                                 }
2119
2120                                                 if (mcarry != 0)
2121                                                         *dP = (uint)mcarry;
2122                                         }
2123                                 }
2124                         }
2125
2126                         /// <summary>
2127                         /// Multiplies the data in x [xOffset:xOffset+xLen] by
2128                         /// y [yOffset:yOffset+yLen] and puts the low mod words into
2129                         /// d [dOffset:dOffset+mod].
2130                         /// </summary>
2131                         /// <remarks>
2132                         /// This code is unsafe! It is the caller's responsibility to make
2133                         /// sure that it is safe to access x [xOffset:xOffset+xLen],
2134                         /// y [yOffset:yOffset+yLen], and d [dOffset:dOffset+mod].
2135                         /// </remarks>
2136                         public static unsafe void MultiplyMod2p32pmod (uint [] x, int xOffset, int xLen, uint [] y, int yOffest, int yLen, uint [] d, int dOffset, int mod)
2137                         {
2138                                 fixed (uint* xx = x, yy = y, dd = d) {
2139                                         uint* xP = xx + xOffset,
2140                                                 xE = xP + xLen,
2141                                                 yB = yy + yOffest,
2142                                                 yE = yB + yLen,
2143                                                 dB = dd + dOffset,
2144                                                 dE = dB + mod;
2145
2146                                         for (; xP < xE; xP++, dB++) {
2147
2148                                                 if (*xP == 0) continue;
2149
2150                                                 ulong mcarry = 0;
2151                                                 uint* dP = dB;
2152                                                 for (uint* yP = yB; yP < yE && dP < dE; yP++, dP++) {
2153                                                         mcarry += ((ulong)*xP * (ulong)*yP) + (ulong)*dP;
2154
2155                                                         *dP = (uint)mcarry;
2156                                                         mcarry >>= 32;
2157                                                 }
2158
2159                                                 if (mcarry != 0 && dP < dE)
2160                                                         *dP = (uint)mcarry;
2161                                         }
2162                                 }
2163                         }
2164
2165                         public static unsafe void SquarePositive (BigInteger bi, ref uint [] wkSpace)
2166                         {
2167                                 uint [] t = wkSpace;
2168                                 wkSpace = bi.data;
2169                                 uint [] d = bi.data;
2170                                 uint dl = bi.length;
2171                                 bi.data = t;
2172
2173                                 fixed (uint* dd = d, tt = t) {
2174
2175                                         uint* ttE = tt + t.Length;
2176                                         // Clear the dest
2177                                         for (uint* ttt = tt; ttt < ttE; ttt++)
2178                                                 *ttt = 0;
2179
2180                                         uint* dP = dd, tP = tt;
2181
2182                                         for (uint i = 0; i < dl; i++, dP++) {
2183                                                 if (*dP == 0)
2184                                                         continue;
2185
2186                                                 ulong mcarry = 0;
2187                                                 uint bi1val = *dP;
2188
2189                                                 uint* dP2 = dP + 1, tP2 = tP + 2*i + 1;
2190
2191                                                 for (uint j = i + 1; j < dl; j++, tP2++, dP2++) {
2192                                                         // k = i + j
2193                                                         mcarry += ((ulong)bi1val * (ulong)*dP2) + *tP2;
2194
2195                                                         *tP2 = (uint)mcarry;
2196                                                         mcarry >>= 32;
2197                                                 }
2198
2199                                                 if (mcarry != 0)
2200                                                         *tP2 = (uint)mcarry;
2201                                         }
2202
2203                                         // Double t. Inlined for speed.
2204
2205                                         tP = tt;
2206
2207                                         uint x, carry = 0;
2208                                         while (tP < ttE) {
2209                                                 x = *tP;
2210                                                 *tP = (x << 1) | carry;
2211                                                 carry = x >> (32 - 1);
2212                                                 tP++;
2213                                         }
2214                                         if (carry != 0) *tP = carry;
2215
2216                                         // Add in the diagnals
2217
2218                                         dP = dd;
2219                                         tP = tt;
2220                                         for (uint* dE = dP + dl; (dP < dE); dP++, tP++) {
2221                                                 ulong val = (ulong)*dP * (ulong)*dP + *tP;
2222                                                 *tP = (uint)val;
2223                                                 val >>= 32;
2224                                                 *(++tP) += (uint)val;
2225                                                 if (*tP < (uint)val) {
2226                                                         uint* tP3 = tP;
2227                                                         // Account for the first carry
2228                                                         (*++tP3)++;
2229
2230                                                         // Keep adding until no carry
2231                                                         while ((*tP3++) == 0)
2232                                                                 (*tP3)++;
2233                                                 }
2234
2235                                         }
2236
2237                                         bi.length <<= 1;
2238
2239                                         // Normalize length
2240                                         while (tt [bi.length-1] == 0 && bi.length > 1) bi.length--;
2241
2242                                 }
2243                         }
2244
2245 /*
2246  * Never called in BigInteger (and part of a private class)
2247  *                      public static bool Double (uint [] u, int l)
2248                         {
2249                                 uint x, carry = 0;
2250                                 uint i = 0;
2251                                 while (i < l) {
2252                                         x = u [i];
2253                                         u [i] = (x << 1) | carry;
2254                                         carry = x >> (32 - 1);
2255                                         i++;
2256                                 }
2257                                 if (carry != 0) u [l] = carry;
2258                                 return carry != 0;
2259                         }*/
2260
2261                         #endregion
2262
2263                         #region Number Theory
2264
2265                         public static BigInteger gcd (BigInteger a, BigInteger b)
2266                         {
2267                                 BigInteger x = a;
2268                                 BigInteger y = b;
2269
2270                                 BigInteger g = y;
2271
2272                                 while (x.length > 1) {
2273                                         g = x;
2274                                         x = y % x;
2275                                         y = g;
2276
2277                                 }
2278                                 if (x == 0) return g;
2279
2280                                 // TODO: should we have something here if we can convert to long?
2281
2282                                 //
2283                                 // Now we can just do it with single precision. I am using the binary gcd method,
2284                                 // as it should be faster.
2285                                 //
2286
2287                                 uint yy = x.data [0];
2288                                 uint xx = y % yy;
2289
2290                                 int t = 0;
2291
2292                                 while (((xx | yy) & 1) == 0) {
2293                                         xx >>= 1; yy >>= 1; t++;
2294                                 }
2295                                 while (xx != 0) {
2296                                         while ((xx & 1) == 0) xx >>= 1;
2297                                         while ((yy & 1) == 0) yy >>= 1;
2298                                         if (xx >= yy)
2299                                                 xx = (xx - yy) >> 1;
2300                                         else
2301                                                 yy = (yy - xx) >> 1;
2302                                 }
2303
2304                                 return yy << t;
2305                         }
2306
2307                         public static uint modInverse (BigInteger bi, uint modulus)
2308                         {
2309                                 uint a = modulus, b = bi % modulus;
2310                                 uint p0 = 0, p1 = 1;
2311
2312                                 while (b != 0) {
2313                                         if (b == 1)
2314                                                 return p1;
2315                                         p0 += (a / b) * p1;
2316                                         a %= b;
2317
2318                                         if (a == 0)
2319                                                 break;
2320                                         if (a == 1)
2321                                                 return modulus-p0;
2322
2323                                         p1 += (b / a) * p0;
2324                                         b %= a;
2325
2326                                 }
2327                                 return 0;
2328                         }
2329
2330                         public static BigInteger modInverse (BigInteger bi, BigInteger modulus)
2331                         {
2332                                 if (modulus.length == 1) return modInverse (bi, modulus.data [0]);
2333
2334                                 BigInteger [] p = { 0, 1 };
2335                                 BigInteger [] q = new BigInteger [2];    // quotients
2336                                 BigInteger [] r = { 0, 0 };             // remainders
2337
2338                                 int step = 0;
2339
2340                                 BigInteger a = modulus;
2341                                 BigInteger b = bi;
2342
2343                                 ModulusRing mr = new ModulusRing (modulus);
2344
2345                                 while (b != 0) {
2346
2347                                         if (step > 1) {
2348
2349                                                 BigInteger pval = mr.Difference (p [0], p [1] * q [0]);
2350                                                 p [0] = p [1]; p [1] = pval;
2351                                         }
2352
2353                                         BigInteger [] divret = multiByteDivide (a, b);
2354
2355                                         q [0] = q [1]; q [1] = divret [0];
2356                                         r [0] = r [1]; r [1] = divret [1];
2357                                         a = b;
2358                                         b = divret [1];
2359
2360                                         step++;
2361                                 }
2362
2363                                 if (r [0] != 1)
2364                                         throw (new ArithmeticException ("No inverse!"));
2365
2366                                 return mr.Difference (p [0], p [1] * q [0]);
2367
2368                         }
2369                         #endregion
2370                 }
2371         }
2372 }