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