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                 internal sealed class Montgomery {
1566
1567                         private Montgomery ()
1568                         {
1569                         }
1570
1571                         public static uint Inverse (uint n)
1572                         {
1573                                 uint y = n, z;
1574
1575                                 while ((z = n * y) != 1)
1576                                         y *= 2 - z;
1577
1578                                 return (uint)-y;
1579                         }
1580
1581                         public static BigInteger ToMont (BigInteger n, BigInteger m)
1582                         {
1583                                 n.Normalize (); m.Normalize ();
1584
1585                                 n <<= (int)m.length * 32;
1586                                 n %= m;
1587                                 return n;
1588                         }
1589
1590                         public static unsafe BigInteger Reduce (BigInteger n, BigInteger m, uint mPrime)
1591                         {
1592                                 BigInteger A = n;
1593                                 fixed (uint* a = A.data, mm = m.data) {
1594                                         for (uint i = 0; i < m.length; i++) {
1595                                                 // The mod here is taken care of by the CPU,
1596                                                 // since the multiply will overflow.
1597                                                 uint u_i = a [0] * mPrime /* % 2^32 */;
1598
1599                                                 //
1600                                                 // A += u_i * m;
1601                                                 // A >>= 32
1602                                                 //
1603
1604                                                 // mP = Position in mod
1605                                                 // aSP = the source of bits from a
1606                                                 // aDP = destination for bits
1607                                                 uint* mP = mm, aSP = a, aDP = a;
1608
1609                                                 ulong c = (ulong)u_i * ((ulong)*(mP++)) + *(aSP++);
1610                                                 c >>= 32;
1611                                                 uint j = 1;
1612
1613                                                 // Multiply and add
1614                                                 for (; j < m.length; j++) {
1615                                                         c += (ulong)u_i * (ulong)*(mP++) + *(aSP++);
1616                                                         *(aDP++) = (uint)c;
1617                                                         c >>= 32;
1618                                                 }
1619
1620                                                 // Account for carry
1621                                                 // TODO: use a better loop here, we dont need the ulong stuff
1622                                                 for (; j < A.length; j++) {
1623                                                         c += *(aSP++);
1624                                                         *(aDP++) = (uint)c;
1625                                                         c >>= 32;
1626                                                         if (c == 0) {j++; break;}
1627                                                 }
1628                                                 // Copy the rest
1629                                                 for (; j < A.length; j++) {
1630                                                         *(aDP++) = *(aSP++);
1631                                                 }
1632
1633                                                 *(aDP++) = (uint)c;
1634                                         }
1635
1636                                         while (A.length > 1 && a [A.length-1] == 0) A.length--;
1637
1638                                 }
1639                                 if (A >= m) Kernel.MinusEq (A, m);
1640
1641                                 return A;
1642                         }
1643 #if _NOT_USED_
1644                         public static BigInteger Reduce (BigInteger n, BigInteger m)
1645                         {
1646                                 return Reduce (n, m, Inverse (m.data [0]));
1647                         }
1648 #endif
1649                 }
1650
1651                 /// <summary>
1652                 /// Low level functions for the BigInteger
1653                 /// </summary>
1654                 private sealed class Kernel {
1655
1656                         #region Addition/Subtraction
1657
1658                         /// <summary>
1659                         /// Adds two numbers with the same sign.
1660                         /// </summary>
1661                         /// <param name="bi1">A BigInteger</param>
1662                         /// <param name="bi2">A BigInteger</param>
1663                         /// <returns>bi1 + bi2</returns>
1664                         public static BigInteger AddSameSign (BigInteger bi1, BigInteger bi2)
1665                         {
1666                                 uint [] x, y;
1667                                 uint yMax, xMax, i = 0;
1668
1669                                 // x should be bigger
1670                                 if (bi1.length < bi2.length) {
1671                                         x = bi2.data;
1672                                         xMax = bi2.length;
1673                                         y = bi1.data;
1674                                         yMax = bi1.length;
1675                                 } else {
1676                                         x = bi1.data;
1677                                         xMax = bi1.length;
1678                                         y = bi2.data;
1679                                         yMax = bi2.length;
1680                                 }
1681
1682                                 BigInteger result = new BigInteger (Sign.Positive, xMax + 1);
1683
1684                                 uint [] r = result.data;
1685
1686                                 ulong sum = 0;
1687
1688                                 // Add common parts of both numbers
1689                                 do {
1690                                         sum = ((ulong)x [i]) + ((ulong)y [i]) + sum;
1691                                         r [i] = (uint)sum;
1692                                         sum >>= 32;
1693                                 } while (++i < yMax);
1694
1695                                 // Copy remainder of longer number while carry propagation is required
1696                                 bool carry = (sum != 0);
1697
1698                                 if (carry) {
1699
1700                                         if (i < xMax) {
1701                                                 do
1702                                                         carry = ((r [i] = x [i] + 1) == 0);
1703                                                 while (++i < xMax && carry);
1704                                         }
1705
1706                                         if (carry) {
1707                                                 r [i] = 1;
1708                                                 result.length = ++i;
1709                                                 return result;
1710                                         }
1711                                 }
1712
1713                                 // Copy the rest
1714                                 if (i < xMax) {
1715                                         do
1716                                                 r [i] = x [i];
1717                                         while (++i < xMax);
1718                                 }
1719
1720                                 result.Normalize ();
1721                                 return result;
1722                         }
1723
1724                         public static BigInteger Subtract (BigInteger big, BigInteger small)
1725                         {
1726                                 BigInteger result = new BigInteger (Sign.Positive, big.length);
1727
1728                                 uint [] r = result.data, b = big.data, s = small.data;
1729                                 uint i = 0, c = 0;
1730
1731                                 do {
1732
1733                                         uint x = s [i];
1734                                         if (((x += c) < c) | ((r [i] = b [i] - x) > ~x))
1735                                                 c = 1;
1736                                         else
1737                                                 c = 0;
1738
1739                                 } while (++i < small.length);
1740
1741                                 if (i == big.length) goto fixup;
1742
1743                                 if (c == 1) {
1744                                         do
1745                                                 r [i] = b [i] - 1;
1746                                         while (b [i++] == 0 && i < big.length);
1747
1748                                         if (i == big.length) goto fixup;
1749                                 }
1750
1751                                 do
1752                                         r [i] = b [i];
1753                                 while (++i < big.length);
1754
1755                                 fixup:
1756
1757                                         result.Normalize ();
1758                                 return result;
1759                         }
1760
1761                         public static void MinusEq (BigInteger big, BigInteger small)
1762                         {
1763                                 uint [] b = big.data, s = small.data;
1764                                 uint i = 0, c = 0;
1765
1766                                 do {
1767                                         uint x = s [i];
1768                                         if (((x += c) < c) | ((b [i] -= x) > ~x))
1769                                                 c = 1;
1770                                         else
1771                                                 c = 0;
1772                                 } while (++i < small.length);
1773
1774                                 if (i == big.length) goto fixup;
1775
1776                                 if (c == 1) {
1777                                         do
1778                                                 b [i]--;
1779                                         while (b [i++] == 0 && i < big.length);
1780                                 }
1781
1782                                 fixup:
1783
1784                                         // Normalize length
1785                                         while (big.length > 0 && big.data [big.length-1] == 0) big.length--;
1786
1787                                 // Check for zero
1788                                 if (big.length == 0)
1789                                         big.length++;
1790
1791                         }
1792
1793                         public static void PlusEq (BigInteger bi1, BigInteger bi2)
1794                         {
1795                                 uint [] x, y;
1796                                 uint yMax, xMax, i = 0;
1797                                 bool flag = false;
1798
1799                                 // x should be bigger
1800                                 if (bi1.length < bi2.length){
1801                                         flag = true;
1802                                         x = bi2.data;
1803                                         xMax = bi2.length;
1804                                         y = bi1.data;
1805                                         yMax = bi1.length;
1806                                 } else {
1807                                         x = bi1.data;
1808                                         xMax = bi1.length;
1809                                         y = bi2.data;
1810                                         yMax = bi2.length;
1811                                 }
1812
1813                                 uint [] r = bi1.data;
1814
1815                                 ulong sum = 0;
1816
1817                                 // Add common parts of both numbers
1818                                 do {
1819                                         sum += ((ulong)x [i]) + ((ulong)y [i]);
1820                                         r [i] = (uint)sum;
1821                                         sum >>= 32;
1822                                 } while (++i < yMax);
1823
1824                                 // Copy remainder of longer number while carry propagation is required
1825                                 bool carry = (sum != 0);
1826
1827                                 if (carry){
1828
1829                                         if (i < xMax) {
1830                                                 do
1831                                                         carry = ((r [i] = x [i] + 1) == 0);
1832                                                 while (++i < xMax && carry);
1833                                         }
1834
1835                                         if (carry) {
1836                                                 r [i] = 1;
1837                                                 bi1.length = ++i;
1838                                                 return;
1839                                         }
1840                                 }
1841
1842                                 // Copy the rest
1843                                 if (flag && i < xMax - 1) {
1844                                         do
1845                                                 r [i] = x [i];
1846                                         while (++i < xMax);
1847                                 }
1848
1849                                 bi1.length = xMax + 1;
1850                                 bi1.Normalize ();
1851                         }
1852
1853                         #endregion
1854
1855                         #region Compare
1856
1857                         /// <summary>
1858                         /// Compares two BigInteger
1859                         /// </summary>
1860                         /// <param name="bi1">A BigInteger</param>
1861                         /// <param name="bi2">A BigInteger</param>
1862                         /// <returns>The sign of bi1 - bi2</returns>
1863                         public static Sign Compare (BigInteger bi1, BigInteger bi2)
1864                         {
1865                                 //
1866                                 // Step 1. Compare the lengths
1867                                 //
1868                                 uint l1 = bi1.length, l2 = bi2.length;
1869
1870                                 while (l1 > 0 && bi1.data [l1-1] == 0) l1--;
1871                                 while (l2 > 0 && bi2.data [l2-1] == 0) l2--;
1872
1873                                 if (l1 == 0 && l2 == 0) return Sign.Zero;
1874
1875                                 // bi1 len < bi2 len
1876                                 if (l1 < l2) return Sign.Negative;
1877                                 // bi1 len > bi2 len
1878                                 else if (l1 > l2) return Sign.Positive;
1879
1880                                 //
1881                                 // Step 2. Compare the bits
1882                                 //
1883
1884                                 uint pos = l1 - 1;
1885
1886                                 while (pos != 0 && bi1.data [pos] == bi2.data [pos]) pos--;
1887
1888                                 if (bi1.data [pos] < bi2.data [pos])
1889                                         return Sign.Negative;
1890                                 else if (bi1.data [pos] > bi2.data [pos])
1891                                         return Sign.Positive;
1892                                 else
1893                                         return Sign.Zero;
1894                         }
1895
1896                         #endregion
1897
1898                         #region Division
1899
1900                         #region Dword
1901
1902                         /// <summary>
1903                         /// Performs n / d and n % d in one operation.
1904                         /// </summary>
1905                         /// <param name="n">A BigInteger, upon exit this will hold n / d</param>
1906                         /// <param name="d">The divisor</param>
1907                         /// <returns>n % d</returns>
1908                         public static uint SingleByteDivideInPlace (BigInteger n, uint d)
1909                         {
1910                                 ulong r = 0;
1911                                 uint i = n.length;
1912
1913                                 while (i-- > 0) {
1914                                         r <<= 32;
1915                                         r |= n.data [i];
1916                                         n.data [i] = (uint)(r / d);
1917                                         r %= d;
1918                                 }
1919                                 n.Normalize ();
1920
1921                                 return (uint)r;
1922                         }
1923
1924                         public static uint DwordMod (BigInteger n, uint d)
1925                         {
1926                                 ulong r = 0;
1927                                 uint i = n.length;
1928
1929                                 while (i-- > 0) {
1930                                         r <<= 32;
1931                                         r |= n.data [i];
1932                                         r %= d;
1933                                 }
1934
1935                                 return (uint)r;
1936                         }
1937
1938                         public static BigInteger DwordDiv (BigInteger n, uint d)
1939                         {
1940                                 BigInteger ret = new BigInteger (Sign.Positive, n.length);
1941
1942                                 ulong r = 0;
1943                                 uint i = n.length;
1944
1945                                 while (i-- > 0) {
1946                                         r <<= 32;
1947                                         r |= n.data [i];
1948                                         ret.data [i] = (uint)(r / d);
1949                                         r %= d;
1950                                 }
1951                                 ret.Normalize ();
1952
1953                                 return ret;
1954                         }
1955
1956                         public static BigInteger [] DwordDivMod (BigInteger n, uint d)
1957                         {
1958                                 BigInteger ret = new BigInteger (Sign.Positive , n.length);
1959
1960                                 ulong r = 0;
1961                                 uint i = n.length;
1962
1963                                 while (i-- > 0) {
1964                                         r <<= 32;
1965                                         r |= n.data [i];
1966                                         ret.data [i] = (uint)(r / d);
1967                                         r %= d;
1968                                 }
1969                                 ret.Normalize ();
1970
1971                                 BigInteger rem = (uint)r;
1972
1973                                 return new BigInteger [] {ret, rem};
1974                         }
1975
1976                                 #endregion
1977
1978                         #region BigNum
1979
1980                         public static BigInteger [] multiByteDivide (BigInteger bi1, BigInteger bi2)
1981                         {
1982                                 if (Kernel.Compare (bi1, bi2) == Sign.Negative)
1983                                         return new BigInteger [2] { 0, new BigInteger (bi1) };
1984
1985                                 bi1.Normalize (); bi2.Normalize ();
1986
1987                                 if (bi2.length == 1)
1988                                         return DwordDivMod (bi1, bi2.data [0]);
1989
1990                                 uint remainderLen = bi1.length + 1;
1991                                 int divisorLen = (int)bi2.length + 1;
1992
1993                                 uint mask = 0x80000000;
1994                                 uint val = bi2.data [bi2.length - 1];
1995                                 int shift = 0;
1996                                 int resultPos = (int)bi1.length - (int)bi2.length;
1997
1998                                 while (mask != 0 && (val & mask) == 0) {
1999                                         shift++; mask >>= 1;
2000                                 }
2001
2002                                 BigInteger quot = new BigInteger (Sign.Positive, bi1.length - bi2.length + 1);
2003                                 BigInteger rem = (bi1 << shift);
2004
2005                                 uint [] remainder = rem.data;
2006
2007                                 bi2 = bi2 << shift;
2008
2009                                 int j = (int)(remainderLen - bi2.length);
2010                                 int pos = (int)remainderLen - 1;
2011
2012                                 uint firstDivisorByte = bi2.data [bi2.length-1];
2013                                 ulong secondDivisorByte = bi2.data [bi2.length-2];
2014
2015                                 while (j > 0) {
2016                                         ulong dividend = ((ulong)remainder [pos] << 32) + (ulong)remainder [pos-1];
2017
2018                                         ulong q_hat = dividend / (ulong)firstDivisorByte;
2019                                         ulong r_hat = dividend % (ulong)firstDivisorByte;
2020
2021                                         do {
2022
2023                                                 if (q_hat == 0x100000000 ||
2024                                                         (q_hat * secondDivisorByte) > ((r_hat << 32) + remainder [pos-2])) {
2025                                                         q_hat--;
2026                                                         r_hat += (ulong)firstDivisorByte;
2027
2028                                                         if (r_hat < 0x100000000)
2029                                                                 continue;
2030                                                 }
2031                                                 break;
2032                                         } while (true);
2033
2034                                         //
2035                                         // At this point, q_hat is either exact, or one too large
2036                                         // (more likely to be exact) so, we attempt to multiply the
2037                                         // divisor by q_hat, if we get a borrow, we just subtract
2038                                         // one from q_hat and add the divisor back.
2039                                         //
2040
2041                                         uint t;
2042                                         uint dPos = 0;
2043                                         int nPos = pos - divisorLen + 1;
2044                                         ulong mc = 0;
2045                                         uint uint_q_hat = (uint)q_hat;
2046                                         do {
2047                                                 mc += (ulong)bi2.data [dPos] * (ulong)uint_q_hat;
2048                                                 t = remainder [nPos];
2049                                                 remainder [nPos] -= (uint)mc;
2050                                                 mc >>= 32;
2051                                                 if (remainder [nPos] > t) mc++;
2052                                                 dPos++; nPos++;
2053                                         } while (dPos < divisorLen);
2054
2055                                         nPos = pos - divisorLen + 1;
2056                                         dPos = 0;
2057
2058                                         // Overestimate
2059                                         if (mc != 0) {
2060                                                 uint_q_hat--;
2061                                                 ulong sum = 0;
2062
2063                                                 do {
2064                                                         sum = ((ulong)remainder [nPos]) + ((ulong)bi2.data [dPos]) + sum;
2065                                                         remainder [nPos] = (uint)sum;
2066                                                         sum >>= 32;
2067                                                         dPos++; nPos++;
2068                                                 } while (dPos < divisorLen);
2069
2070                                         }
2071
2072                                         quot.data [resultPos--] = (uint)uint_q_hat;
2073
2074                                         pos--;
2075                                         j--;
2076                                 }
2077
2078                                 quot.Normalize ();
2079                                 rem.Normalize ();
2080                                 BigInteger [] ret = new BigInteger [2] { quot, rem };
2081
2082                                 if (shift != 0)
2083                                         ret [1] >>= shift;
2084
2085                                 return ret;
2086                         }
2087
2088                         #endregion
2089
2090                         #endregion
2091
2092                         #region Shift
2093                         public static BigInteger LeftShift (BigInteger bi, int n)
2094                         {
2095                                 if (n == 0) return new BigInteger (bi, bi.length + 1);
2096
2097                                 int w = n >> 5;
2098                                 n &= ((1 << 5) - 1);
2099
2100                                 BigInteger ret = new BigInteger (Sign.Positive, bi.length + 1 + (uint)w);
2101
2102                                 uint i = 0, l = bi.length;
2103                                 if (n != 0) {
2104                                         uint x, carry = 0;
2105                                         while (i < l) {
2106                                                 x = bi.data [i];
2107                                                 ret.data [i + w] = (x << n) | carry;
2108                                                 carry = x >> (32 - n);
2109                                                 i++;
2110                                         }
2111                                         ret.data [i + w] = carry;
2112                                 } else {
2113                                         while (i < l) {
2114                                                 ret.data [i + w] = bi.data [i];
2115                                                 i++;
2116                                         }
2117                                 }
2118
2119                                 ret.Normalize ();
2120                                 return ret;
2121                         }
2122
2123                         public static BigInteger RightShift (BigInteger bi, int n)
2124                         {
2125                                 if (n == 0) return new BigInteger (bi);
2126
2127                                 int w = n >> 5;
2128                                 int s = n & ((1 << 5) - 1);
2129
2130                                 BigInteger ret = new BigInteger (Sign.Positive, bi.length - (uint)w + 1);
2131                                 uint l = (uint)ret.data.Length - 1;
2132
2133                                 if (s != 0) {
2134
2135                                         uint x, carry = 0;
2136
2137                                         while (l-- > 0) {
2138                                                 x = bi.data [l + w];
2139                                                 ret.data [l] = (x >> n) | carry;
2140                                                 carry = x << (32 - n);
2141                                         }
2142                                 } else {
2143                                         while (l-- > 0)
2144                                                 ret.data [l] = bi.data [l + w];
2145
2146                                 }
2147                                 ret.Normalize ();
2148                                 return ret;
2149                         }
2150
2151                         #endregion
2152
2153                         #region Multiply
2154
2155                         public static BigInteger MultiplyByDword (BigInteger n, uint f)
2156                         {
2157                                 BigInteger ret = new BigInteger (Sign.Positive, n.length + 1);
2158
2159                                 uint i = 0;
2160                                 ulong c = 0;
2161
2162                                 do {
2163                                         c += (ulong)n.data [i] * (ulong)f;
2164                                         ret.data [i] = (uint)c;
2165                                         c >>= 32;
2166                                 } while (++i < n.length);
2167                                 ret.data [i] = (uint)c;
2168                                 ret.Normalize ();
2169                                 return ret;
2170
2171                         }
2172
2173                         /// <summary>
2174                         /// Multiplies the data in x [xOffset:xOffset+xLen] by
2175                         /// y [yOffset:yOffset+yLen] and puts it into
2176                         /// d [dOffset:dOffset+xLen+yLen].
2177                         /// </summary>
2178                         /// <remarks>
2179                         /// This code is unsafe! It is the caller's responsibility to make
2180                         /// sure that it is safe to access x [xOffset:xOffset+xLen],
2181                         /// y [yOffset:yOffset+yLen], and d [dOffset:dOffset+xLen+yLen].
2182                         /// </remarks>
2183                         public static unsafe void Multiply (uint [] x, uint xOffset, uint xLen, uint [] y, uint yOffset, uint yLen, uint [] d, uint dOffset)
2184                         {
2185                                 fixed (uint* xx = x, yy = y, dd = d) {
2186                                         uint* xP = xx + xOffset,
2187                                                 xE = xP + xLen,
2188                                                 yB = yy + yOffset,
2189                                                 yE = yB + yLen,
2190                                                 dB = dd + dOffset;
2191
2192                                         for (; xP < xE; xP++, dB++) {
2193
2194                                                 if (*xP == 0) continue;
2195
2196                                                 ulong mcarry = 0;
2197
2198                                                 uint* dP = dB;
2199                                                 for (uint* yP = yB; yP < yE; yP++, dP++) {
2200                                                         mcarry += ((ulong)*xP * (ulong)*yP) + (ulong)*dP;
2201
2202                                                         *dP = (uint)mcarry;
2203                                                         mcarry >>= 32;
2204                                                 }
2205
2206                                                 if (mcarry != 0)
2207                                                         *dP = (uint)mcarry;
2208                                         }
2209                                 }
2210                         }
2211
2212                         /// <summary>
2213                         /// Multiplies the data in x [xOffset:xOffset+xLen] by
2214                         /// y [yOffset:yOffset+yLen] and puts the low mod words into
2215                         /// d [dOffset:dOffset+mod].
2216                         /// </summary>
2217                         /// <remarks>
2218                         /// This code is unsafe! It is the caller's responsibility to make
2219                         /// sure that it is safe to access x [xOffset:xOffset+xLen],
2220                         /// y [yOffset:yOffset+yLen], and d [dOffset:dOffset+mod].
2221                         /// </remarks>
2222                         public static unsafe void MultiplyMod2p32pmod (uint [] x, int xOffset, int xLen, uint [] y, int yOffest, int yLen, uint [] d, int dOffset, int mod)
2223                         {
2224                                 fixed (uint* xx = x, yy = y, dd = d) {
2225                                         uint* xP = xx + xOffset,
2226                                                 xE = xP + xLen,
2227                                                 yB = yy + yOffest,
2228                                                 yE = yB + yLen,
2229                                                 dB = dd + dOffset,
2230                                                 dE = dB + mod;
2231
2232                                         for (; xP < xE; xP++, dB++) {
2233
2234                                                 if (*xP == 0) continue;
2235
2236                                                 ulong mcarry = 0;
2237                                                 uint* dP = dB;
2238                                                 for (uint* yP = yB; yP < yE && dP < dE; yP++, dP++) {
2239                                                         mcarry += ((ulong)*xP * (ulong)*yP) + (ulong)*dP;
2240
2241                                                         *dP = (uint)mcarry;
2242                                                         mcarry >>= 32;
2243                                                 }
2244
2245                                                 if (mcarry != 0 && dP < dE)
2246                                                         *dP = (uint)mcarry;
2247                                         }
2248                                 }
2249                         }
2250
2251                         public static unsafe void SquarePositive (BigInteger bi, ref uint [] wkSpace)
2252                         {
2253                                 uint [] t = wkSpace;
2254                                 wkSpace = bi.data;
2255                                 uint [] d = bi.data;
2256                                 uint dl = bi.length;
2257                                 bi.data = t;
2258
2259                                 fixed (uint* dd = d, tt = t) {
2260
2261                                         uint* ttE = tt + t.Length;
2262                                         // Clear the dest
2263                                         for (uint* ttt = tt; ttt < ttE; ttt++)
2264                                                 *ttt = 0;
2265
2266                                         uint* dP = dd, tP = tt;
2267
2268                                         for (uint i = 0; i < dl; i++, dP++) {
2269                                                 if (*dP == 0)
2270                                                         continue;
2271
2272                                                 ulong mcarry = 0;
2273                                                 uint bi1val = *dP;
2274
2275                                                 uint* dP2 = dP + 1, tP2 = tP + 2*i + 1;
2276
2277                                                 for (uint j = i + 1; j < dl; j++, tP2++, dP2++) {
2278                                                         // k = i + j
2279                                                         mcarry += ((ulong)bi1val * (ulong)*dP2) + *tP2;
2280
2281                                                         *tP2 = (uint)mcarry;
2282                                                         mcarry >>= 32;
2283                                                 }
2284
2285                                                 if (mcarry != 0)
2286                                                         *tP2 = (uint)mcarry;
2287                                         }
2288
2289                                         // Double t. Inlined for speed.
2290
2291                                         tP = tt;
2292
2293                                         uint x, carry = 0;
2294                                         while (tP < ttE) {
2295                                                 x = *tP;
2296                                                 *tP = (x << 1) | carry;
2297                                                 carry = x >> (32 - 1);
2298                                                 tP++;
2299                                         }
2300                                         if (carry != 0) *tP = carry;
2301
2302                                         // Add in the diagnals
2303
2304                                         dP = dd;
2305                                         tP = tt;
2306                                         for (uint* dE = dP + dl; (dP < dE); dP++, tP++) {
2307                                                 ulong val = (ulong)*dP * (ulong)*dP + *tP;
2308                                                 *tP = (uint)val;
2309                                                 val >>= 32;
2310                                                 *(++tP) += (uint)val;
2311                                                 if (*tP < (uint)val) {
2312                                                         uint* tP3 = tP;
2313                                                         // Account for the first carry
2314                                                         (*++tP3)++;
2315
2316                                                         // Keep adding until no carry
2317                                                         while ((*tP3++) == 0)
2318                                                                 (*tP3)++;
2319                                                 }
2320
2321                                         }
2322
2323                                         bi.length <<= 1;
2324
2325                                         // Normalize length
2326                                         while (tt [bi.length-1] == 0 && bi.length > 1) bi.length--;
2327
2328                                 }
2329                         }
2330
2331 /*
2332  * Never called in BigInteger (and part of a private class)
2333  *                      public static bool Double (uint [] u, int l)
2334                         {
2335                                 uint x, carry = 0;
2336                                 uint i = 0;
2337                                 while (i < l) {
2338                                         x = u [i];
2339                                         u [i] = (x << 1) | carry;
2340                                         carry = x >> (32 - 1);
2341                                         i++;
2342                                 }
2343                                 if (carry != 0) u [l] = carry;
2344                                 return carry != 0;
2345                         }*/
2346
2347                         #endregion
2348
2349                         #region Number Theory
2350
2351                         public static BigInteger gcd (BigInteger a, BigInteger b)
2352                         {
2353                                 BigInteger x = a;
2354                                 BigInteger y = b;
2355
2356                                 BigInteger g = y;
2357
2358                                 while (x.length > 1) {
2359                                         g = x;
2360                                         x = y % x;
2361                                         y = g;
2362
2363                                 }
2364                                 if (x == 0) return g;
2365
2366                                 // TODO: should we have something here if we can convert to long?
2367
2368                                 //
2369                                 // Now we can just do it with single precision. I am using the binary gcd method,
2370                                 // as it should be faster.
2371                                 //
2372
2373                                 uint yy = x.data [0];
2374                                 uint xx = y % yy;
2375
2376                                 int t = 0;
2377
2378                                 while (((xx | yy) & 1) == 0) {
2379                                         xx >>= 1; yy >>= 1; t++;
2380                                 }
2381                                 while (xx != 0) {
2382                                         while ((xx & 1) == 0) xx >>= 1;
2383                                         while ((yy & 1) == 0) yy >>= 1;
2384                                         if (xx >= yy)
2385                                                 xx = (xx - yy) >> 1;
2386                                         else
2387                                                 yy = (yy - xx) >> 1;
2388                                 }
2389
2390                                 return yy << t;
2391                         }
2392
2393                         public static uint modInverse (BigInteger bi, uint modulus)
2394                         {
2395                                 uint a = modulus, b = bi % modulus;
2396                                 uint p0 = 0, p1 = 1;
2397
2398                                 while (b != 0) {
2399                                         if (b == 1)
2400                                                 return p1;
2401                                         p0 += (a / b) * p1;
2402                                         a %= b;
2403
2404                                         if (a == 0)
2405                                                 break;
2406                                         if (a == 1)
2407                                                 return modulus-p0;
2408
2409                                         p1 += (b / a) * p0;
2410                                         b %= a;
2411
2412                                 }
2413                                 return 0;
2414                         }
2415
2416                         public static BigInteger modInverse (BigInteger bi, BigInteger modulus)
2417                         {
2418                                 if (modulus.length == 1) return modInverse (bi, modulus.data [0]);
2419
2420                                 BigInteger [] p = { 0, 1 };
2421                                 BigInteger [] q = new BigInteger [2];    // quotients
2422                                 BigInteger [] r = { 0, 0 };             // remainders
2423
2424                                 int step = 0;
2425
2426                                 BigInteger a = modulus;
2427                                 BigInteger b = bi;
2428
2429                                 ModulusRing mr = new ModulusRing (modulus);
2430
2431                                 while (b != 0) {
2432
2433                                         if (step > 1) {
2434
2435                                                 BigInteger pval = mr.Difference (p [0], p [1] * q [0]);
2436                                                 p [0] = p [1]; p [1] = pval;
2437                                         }
2438
2439                                         BigInteger [] divret = multiByteDivide (a, b);
2440
2441                                         q [0] = q [1]; q [1] = divret [0];
2442                                         r [0] = r [1]; r [1] = divret [1];
2443                                         a = b;
2444                                         b = divret [1];
2445
2446                                         step++;
2447                                 }
2448
2449                                 if (r [0] != 1)
2450                                         throw (new ArithmeticException ("No inverse!"));
2451
2452                                 return mr.Difference (p [0], p [1] * q [0]);
2453
2454                         }
2455                         #endregion
2456                 }
2457         }
2458 }