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