2 // BigInteger.cs - Big Integer implementation
7 // Sebastien Pouliot <sebastien@ximian.com>
8 // Pieter Philippaerts <Pieter@mentalis.org>
10 // Copyright (c) 2003 Ben Maurer
11 // All rights reserved
13 // Copyright (c) 2002 Chew Keong TAN
14 // All rights reserved.
16 // Copyright (C) 2004, 2007 Novell, Inc (http://www.novell.com)
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:
26 // The above copyright notice and this permission notice shall be
27 // included in all copies or substantial portions of the Software.
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.
39 using System.Security.Cryptography;
40 using Mono.Math.Prime.Generator;
41 using Mono.Math.Prime;
55 /// The Length of this BigInteger
60 /// The data for this BigInteger
69 /// Default length of a BigInteger in bytes
71 const uint DEFAULT_LEN = 20;
74 /// Table of primes below 2000.
78 /// This table was generated using Mathematica 4.1 using the following function:
82 /// PrimeTable [x_] := Prime [Range [1, PrimePi [x]]]
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,
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,
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,
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,
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,
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
156 public enum Sign : int {
162 #region Exception Messages
163 const string WouldReturnNegVal = "Operation would return a negative value";
172 data = new uint [DEFAULT_LEN];
173 this.length = DEFAULT_LEN;
177 [CLSCompliant (false)]
179 public BigInteger (Sign sign, uint len)
181 this.data = new uint [len];
185 public BigInteger (BigInteger bi)
187 this.data = (uint [])bi.data.Clone ();
188 this.length = bi.length;
192 [CLSCompliant (false)]
194 public BigInteger (BigInteger bi, uint len)
197 this.data = new uint [len];
199 for (uint i = 0; i < bi.length; i++)
200 this.data [i] = bi.data [i];
202 this.length = bi.length;
209 public BigInteger (byte [] inData)
211 length = (uint)inData.Length >> 2;
212 int leftOver = inData.Length & 0x3;
214 // length not multiples of 4
215 if (leftOver != 0) length++;
217 data = new uint [length];
219 for (int i = inData.Length - 1, j = 0; i >= 3; i -= 4, j++) {
221 (inData [i-3] << (3*8)) |
222 (inData [i-2] << (2*8)) |
223 (inData [i-1] << (1*8)) |
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;
238 [CLSCompliant (false)]
240 public BigInteger (uint [] inData)
242 length = (uint)inData.Length;
244 data = new uint [length];
246 for (int i = (int)length - 1, j = 0; i >= 0; i--, j++)
247 data [j] = inData [i];
253 [CLSCompliant (false)]
255 public BigInteger (uint ui)
257 data = new uint [] {ui};
261 [CLSCompliant (false)]
263 public BigInteger (ulong ul)
265 data = new uint [2] { (uint)ul, (uint)(ul >> 32)};
272 [CLSCompliant (false)]
274 public static implicit operator BigInteger (uint value)
276 return (new BigInteger (value));
279 public static implicit operator BigInteger (int value)
281 if (value < 0) throw new ArgumentOutOfRangeException ("value");
282 return (new BigInteger ((uint)value));
286 [CLSCompliant (false)]
288 public static implicit operator BigInteger (ulong value)
290 return (new BigInteger (value));
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)
298 throw new ArgumentNullException ("number");
300 int i = 0, len = number.Length;
302 bool digits_seen = false;
303 BigInteger val = new BigInteger (0);
304 if (number [i] == '+') {
307 else if (number [i] == '-') {
308 throw new FormatException (WouldReturnNegVal);
311 for (; i < len; i++) {
317 if (c >= '0' && c <= '9') {
318 val = val * 10 + (c - '0');
322 if (Char.IsWhiteSpace (c)) {
323 for (i++; i < len; i++) {
324 if (!Char.IsWhiteSpace (number [i]))
325 throw new FormatException ();
330 throw new FormatException ();
334 throw new FormatException ();
342 public static BigInteger operator + (BigInteger bi1, BigInteger bi2)
345 return new BigInteger (bi2);
347 return new BigInteger (bi1);
349 return Kernel.AddSameSign (bi1, bi2);
352 public static BigInteger operator - (BigInteger bi1, BigInteger bi2)
355 return new BigInteger (bi1);
358 throw new ArithmeticException (WouldReturnNegVal);
360 switch (Kernel.Compare (bi1, bi2)) {
366 return Kernel.Subtract (bi1, bi2);
369 throw new ArithmeticException (WouldReturnNegVal);
371 throw new Exception ();
375 public static int operator % (BigInteger bi, int i)
378 return (int)Kernel.DwordMod (bi, (uint)i);
380 return -(int)Kernel.DwordMod (bi, (uint)-i);
384 [CLSCompliant (false)]
386 public static uint operator % (BigInteger bi, uint ui)
388 return Kernel.DwordMod (bi, (uint)ui);
391 public static BigInteger operator % (BigInteger bi1, BigInteger bi2)
393 return Kernel.multiByteDivide (bi1, bi2)[1];
396 public static BigInteger operator / (BigInteger bi, int i)
399 return Kernel.DwordDiv (bi, (uint)i);
401 throw new ArithmeticException (WouldReturnNegVal);
404 public static BigInteger operator / (BigInteger bi1, BigInteger bi2)
406 return Kernel.multiByteDivide (bi1, bi2)[0];
409 public static BigInteger operator * (BigInteger bi1, BigInteger bi2)
411 if (bi1 == 0 || bi2 == 0) return 0;
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");
419 BigInteger ret = new BigInteger (Sign.Positive, bi1.length + bi2.length);
421 Kernel.Multiply (bi1.data, 0, bi1.length, bi2.data, 0, bi2.length, ret.data, 0);
427 public static BigInteger operator * (BigInteger bi, int i)
429 if (i < 0) throw new ArithmeticException (WouldReturnNegVal);
430 if (i == 0) return 0;
431 if (i == 1) return new BigInteger (bi);
433 return Kernel.MultiplyByDword (bi, (uint)i);
436 public static BigInteger operator << (BigInteger bi1, int shiftVal)
438 return Kernel.LeftShift (bi1, shiftVal);
441 public static BigInteger operator >> (BigInteger bi1, int shiftVal)
443 return Kernel.RightShift (bi1, shiftVal);
448 #region Friendly names for operators
450 // with names suggested by FxCop 1.30
452 public static BigInteger Add (BigInteger bi1, BigInteger bi2)
457 public static BigInteger Subtract (BigInteger bi1, BigInteger bi2)
462 public static int Modulus (BigInteger bi, int i)
468 [CLSCompliant (false)]
470 public static uint Modulus (BigInteger bi, uint ui)
475 public static BigInteger Modulus (BigInteger bi1, BigInteger bi2)
480 public static BigInteger Divid (BigInteger bi, int i)
485 public static BigInteger Divid (BigInteger bi1, BigInteger bi2)
490 public static BigInteger Multiply (BigInteger bi1, BigInteger bi2)
495 public static BigInteger Multiply (BigInteger bi, int i)
503 private static RandomNumberGenerator rng;
504 private static RandomNumberGenerator Rng {
507 rng = RandomNumberGenerator.Create ();
513 /// Generates a new, random BigInteger of the specified length.
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)
520 int dwords = bits >> 5;
521 int remBits = bits & 0x1F;
526 BigInteger ret = new BigInteger (Sign.Positive, (uint)dwords + 1);
527 byte [] random = new byte [dwords << 2];
529 rng.GetBytes (random);
530 Buffer.BlockCopy (random, 0, ret.data, 0, (int)dwords << 2);
533 uint mask = (uint)(0x01 << (remBits-1));
534 ret.data [dwords-1] |= mask;
536 mask = (uint)(0xFFFFFFFF >> (32 - remBits));
537 ret.data [dwords-1] &= mask;
540 ret.data [dwords-1] |= 0x80000000;
547 /// Generates a new, random BigInteger of the specified length using the default RNG crypto service provider.
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)
553 return GenerateRandom (bits, Rng);
557 /// Randomizes the bits in "this" from the specified RNG.
559 /// <param name="rng">A RNG.</param>
560 public void Randomize (RandomNumberGenerator rng)
565 int bits = this.BitCount ();
566 int dwords = bits >> 5;
567 int remBits = bits & 0x1F;
572 byte [] random = new byte [dwords << 2];
574 rng.GetBytes (random);
575 Buffer.BlockCopy (random, 0, data, 0, (int)dwords << 2);
578 uint mask = (uint)(0x01 << (remBits-1));
579 data [dwords-1] |= mask;
581 mask = (uint)(0xFFFFFFFF >> (32 - remBits));
582 data [dwords-1] &= mask;
586 data [dwords-1] |= 0x80000000;
592 /// Randomizes the bits in "this" from the default RNG.
594 public void Randomize ()
603 public int BitCount ()
607 uint value = data [length - 1];
608 uint mask = 0x80000000;
611 while (bits > 0 && (value & mask) == 0) {
615 bits += ((length - 1) << 5);
621 /// Tests if the specified bit is 1.
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>
626 [CLSCompliant (false)]
628 public bool TestBit (uint bitNum)
630 uint bytePos = bitNum >> 5; // divide by 32
631 byte bitPos = (byte)(bitNum & 0x1F); // get the lowest 5 bits
633 uint mask = (uint)1 << bitPos;
634 return ((this.data [bytePos] & mask) != 0);
637 public bool TestBit (int bitNum)
639 if (bitNum < 0) throw new IndexOutOfRangeException ("bitNum out of range");
641 uint bytePos = (uint)bitNum >> 5; // divide by 32
642 byte bitPos = (byte)(bitNum & 0x1F); // get the lowest 5 bits
644 uint mask = (uint)1 << bitPos;
645 return ((this.data [bytePos] | mask) == this.data [bytePos]);
649 [CLSCompliant (false)]
651 public void SetBit (uint bitNum)
653 SetBit (bitNum, true);
657 [CLSCompliant (false)]
659 public void ClearBit (uint bitNum)
661 SetBit (bitNum, false);
665 [CLSCompliant (false)]
667 public void SetBit (uint bitNum, bool value)
669 uint bytePos = bitNum >> 5; // divide by 32
671 if (bytePos < this.length) {
672 uint mask = (uint)1 << (int)(bitNum & 0x1F);
674 this.data [bytePos] |= mask;
676 this.data [bytePos] &= ~mask;
680 public int LowestSetBit ()
682 if (this == 0) return -1;
684 while (!TestBit (i)) i++;
688 public byte[] GetBytes ()
690 if (this == 0) return new byte [1];
692 int numBits = BitCount ();
693 int numBytes = numBits >> 3;
694 if ((numBits & 0x7) != 0)
697 byte [] result = new byte [numBytes];
699 int numBytesInWord = numBytes & 0x3;
700 if (numBytesInWord == 0) numBytesInWord = 4;
703 for (int i = (int)length - 1; i >= 0; i--) {
705 for (int j = numBytesInWord - 1; j >= 0; j--) {
706 result [pos+j] = (byte)(val & 0xFF);
709 pos += numBytesInWord;
720 [CLSCompliant (false)]
722 public static bool operator == (BigInteger bi1, uint ui)
724 if (bi1.length != 1) bi1.Normalize ();
725 return bi1.length == 1 && bi1.data [0] == ui;
729 [CLSCompliant (false)]
731 public static bool operator != (BigInteger bi1, uint ui)
733 if (bi1.length != 1) bi1.Normalize ();
734 return !(bi1.length == 1 && bi1.data [0] == ui);
737 public static bool operator == (BigInteger bi1, BigInteger bi2)
739 // we need to compare with null
740 if ((bi1 as object) == (bi2 as object))
742 if (null == bi1 || null == bi2)
744 return Kernel.Compare (bi1, bi2) == 0;
747 public static bool operator != (BigInteger bi1, BigInteger bi2)
749 // we need to compare with null
750 if ((bi1 as object) == (bi2 as object))
752 if (null == bi1 || null == bi2)
754 return Kernel.Compare (bi1, bi2) != 0;
757 public static bool operator > (BigInteger bi1, BigInteger bi2)
759 return Kernel.Compare (bi1, bi2) > 0;
762 public static bool operator < (BigInteger bi1, BigInteger bi2)
764 return Kernel.Compare (bi1, bi2) < 0;
767 public static bool operator >= (BigInteger bi1, BigInteger bi2)
769 return Kernel.Compare (bi1, bi2) >= 0;
772 public static bool operator <= (BigInteger bi1, BigInteger bi2)
774 return Kernel.Compare (bi1, bi2) <= 0;
777 public Sign Compare (BigInteger bi)
779 return Kernel.Compare (this, bi);
787 [CLSCompliant (false)]
789 public string ToString (uint radix)
791 return ToString (radix, "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ");
795 [CLSCompliant (false)]
797 public string ToString (uint radix, string characterSet)
799 if (characterSet.Length < radix)
800 throw new ArgumentException ("charSet length less than radix", "characterSet");
802 throw new ArgumentException ("There is no such thing as radix one notation", "radix");
804 if (this == 0) return "0";
805 if (this == 1) return "1";
809 BigInteger a = new BigInteger (this);
812 uint rem = Kernel.SingleByteDivideInPlace (a, radix);
813 result = characterSet [(int) rem] + result;
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.
828 private void Normalize ()
831 while (length > 0 && data [length-1] == 0) length--;
840 for (int i=0; i < length; i++)
848 public override int GetHashCode ()
852 for (uint i = 0; i < this.length; i++)
853 val ^= this.data [i];
858 public override string ToString ()
860 return ToString (10);
863 public override bool Equals (object o)
868 return (int)o >= 0 && this == (uint)o;
870 BigInteger bi = o as BigInteger;
874 return Kernel.Compare (this, bi) == 0;
879 #region Number Theory
881 public BigInteger GCD (BigInteger bi)
883 return Kernel.gcd (this, bi);
886 public BigInteger ModInverse (BigInteger modulus)
888 return Kernel.modInverse (this, modulus);
891 public BigInteger ModPow (BigInteger exp, BigInteger n)
893 ModulusRing mr = new ModulusRing (n);
894 return mr.Pow (this, exp);
899 #region Prime Testing
901 public bool IsProbablePrime ()
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])
909 // the list is complete, so it's not a prime
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)
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);
924 #region Prime Number Generation
927 /// Generates the smallest prime >= bi
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)
933 NextPrimeFinder npf = new NextPrimeFinder ();
934 return npf.GenerateNewPrime (0, bi);
937 public static BigInteger GeneratePseudoPrime (int bits)
939 SequentialSearchPrimeGeneratorBase sspg = new SequentialSearchPrimeGeneratorBase ();
940 return sspg.GenerateNewPrime (bits);
944 /// Increments this by two
952 // If there was no carry, nothing to do
955 // Account for the first carry
958 // Keep adding until no carry
959 while (data [i++] == 0x0)
962 // See if we increased the data length
963 if (length == (uint)i)
975 sealed class ModulusRing {
977 BigInteger mod, constant;
979 public ModulusRing (BigInteger modulus)
983 // calculate constant = b^ (2k) / m
984 uint i = mod.length << 1;
986 constant = new BigInteger (Sign.Positive, i + 1);
987 constant.data [i] = 0x00000001;
989 constant = constant / mod;
992 public void BarrettReduction (BigInteger x)
999 // x < mod, so nothing to do.
1000 if (x.length < k) return;
1005 // Validate pointers
1007 if (x.data.Length < x.length) throw new IndexOutOfRangeException ("x out of range");
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
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);
1017 // r1 = x mod b^ (k+1)
1018 // i.e. keep the lowest (k+1) words
1020 uint lengthToCopy = (x.length > kPlusOne) ? kPlusOne : x.length;
1022 x.length = lengthToCopy;
1025 // r2 = (q3 * n) mod b^ (k+1)
1026 // partial multiplication of q3 and n
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);
1034 Kernel.MinusEq (x, r2);
1036 BigInteger val = new BigInteger (Sign.Positive, kPlusOne + 1);
1037 val.data [kPlusOne] = 0x00000001;
1039 Kernel.MinusEq (val, r2);
1040 Kernel.PlusEq (x, val);
1044 Kernel.MinusEq (x, n);
1047 public BigInteger Multiply (BigInteger a, BigInteger b)
1049 if (a == 0 || b == 0) return 0;
1057 BigInteger ret = new BigInteger (a * b);
1058 BarrettReduction (ret);
1063 public BigInteger Difference (BigInteger a, BigInteger b)
1065 Sign cmp = Kernel.Compare (a, b);
1072 diff = a - b; break;
1074 diff = b - a; break;
1076 throw new Exception ();
1080 if (diff.length >= mod.length << 1)
1083 BarrettReduction (diff);
1085 if (cmp == Sign.Negative)
1090 public BigInteger Pow (BigInteger a, BigInteger k)
1092 BigInteger b = new BigInteger (1);
1100 for (int i = 1; i < k.BitCount (); i++) {
1101 A = Multiply (A, A);
1103 b = Multiply (A, b);
1108 public BigInteger Pow (BigInteger b, BigInteger exp)
1110 if ((mod.data [0] & 1) == 1) return OddPow (b, exp);
1111 else return EvenPow (b, exp);
1114 public BigInteger EvenPow (BigInteger b, BigInteger exp)
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)
1119 uint totalBits = (uint)exp.BitCount ();
1121 uint [] wkspace = new uint [mod.length << 1];
1123 // perform squaring and multiply exponentiation
1124 for (uint pos = 0; pos < totalBits; pos++) {
1125 if (exp.TestBit (pos)) {
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;
1134 BarrettReduction (resultNum);
1137 Kernel.SquarePositive (tempNum, ref wkspace);
1138 BarrettReduction (tempNum);
1148 private BigInteger OddPow (BigInteger b, BigInteger exp)
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 ();
1155 uint [] wkspace = new uint [mod.length << 1];
1157 // perform squaring and multiply exponentiation
1158 for (uint pos = 0; pos < totalBits; pos++) {
1159 if (exp.TestBit (pos)) {
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;
1168 Montgomery.Reduce (resultNum, mod, mPrime);
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);
1178 Montgomery.Reduce (resultNum, mod, mPrime);
1182 #region Pow Small Base
1184 // TODO: Make tests for this, not really needed b/c prime stuff
1185 // checks it, but still would be nice
1187 [CLSCompliant (false)]
1190 public BigInteger Pow (uint b, BigInteger exp)
1192 return Pow (new BigInteger (b), exp);
1195 public BigInteger Pow (uint b, BigInteger exp)
1198 if ((mod.data [0] & 1) == 1)
1199 return OddPow (b, exp);
1201 return EvenPow (b, exp);
1202 /* buggy in some cases (like the well tested primes)
1204 if ((mod.data [0] & 1) == 1)
1205 return OddModTwoPow (exp);
1207 return EvenModTwoPow (exp);
1211 private unsafe BigInteger OddPow (uint b, BigInteger exp)
1214 uint [] wkspace = new uint [mod.length << 1 + 1];
1216 BigInteger resultNum = Montgomery.ToMont ((BigInteger)b, this.mod);
1217 resultNum = new BigInteger (resultNum, mod.length << 1 +1);
1219 uint mPrime = Montgomery.Inverse (mod.data [0]);
1221 int bc = exp.BitCount () - 2;
1222 uint pos = (bc > 1 ? (uint) bc : 1);
1225 // We know that the first itr will make the val b
1232 Kernel.SquarePositive (resultNum, ref wkspace);
1233 resultNum = Montgomery.Reduce (resultNum, mod, mPrime);
1235 if (exp.TestBit (pos)) {
1241 // TODO: Is Unsafe really speeding things up?
1242 fixed (uint* u = resultNum.data) {
1248 mc += (ulong)u [i] * (ulong)b;
1251 } while (++i < resultNum.length);
1253 if (resultNum.length < mod.length) {
1257 while (resultNum >= mod)
1258 Kernel.MinusEq (resultNum, mod);
1260 } else if (mc != 0) {
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.
1270 // We would rather have this estimate overshoot,
1271 // so we add one to the divisor
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));
1278 // guess but don't divide by 0
1279 divEstimate = (uint) ((((ulong)cc << 32) | (ulong) u [i -1]) /
1280 (mod.data [mod.length-1]));
1288 mc += (ulong)mod.data [i] * (ulong)divEstimate;
1292 if (u [i] > t) mc++;
1294 } while (i < resultNum.length);
1300 uint [] s = mod.data;
1303 if (((a += sc) < sc) | ((u [j] -= a) > ~a)) sc = 1;
1306 } while (j < resultNum.length);
1309 while (resultNum >= mod)
1310 Kernel.MinusEq (resultNum, mod);
1312 while (resultNum >= mod)
1313 Kernel.MinusEq (resultNum, mod);
1317 } while (pos-- > 0);
1319 resultNum = Montgomery.Reduce (resultNum, mod, mPrime);
1324 private unsafe BigInteger EvenPow (uint b, BigInteger exp)
1327 uint [] wkspace = new uint [mod.length << 1 + 1];
1328 BigInteger resultNum = new BigInteger ((BigInteger)b, mod.length << 1 + 1);
1330 uint pos = (uint)exp.BitCount () - 2;
1333 // We know that the first itr will make the val b
1340 Kernel.SquarePositive (resultNum, ref wkspace);
1341 if (!(resultNum.length < mod.length))
1342 BarrettReduction (resultNum);
1344 if (exp.TestBit (pos)) {
1350 // TODO: Is Unsafe really speeding things up?
1351 fixed (uint* u = resultNum.data) {
1357 mc += (ulong)u [i] * (ulong)b;
1360 } while (++i < resultNum.length);
1362 if (resultNum.length < mod.length) {
1366 while (resultNum >= mod)
1367 Kernel.MinusEq (resultNum, mod);
1369 } else if (mc != 0) {
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.
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));
1389 mc += (ulong)mod.data [i] * (ulong)divEstimate;
1393 if (u [i] > t) mc++;
1395 } while (i < resultNum.length);
1401 uint [] s = mod.data;
1404 if (((a += sc) < sc) | ((u [j] -= a) > ~a)) sc = 1;
1407 } while (j < resultNum.length);
1410 while (resultNum >= mod)
1411 Kernel.MinusEq (resultNum, mod);
1413 while (resultNum >= mod)
1414 Kernel.MinusEq (resultNum, mod);
1418 } while (pos-- > 0);
1423 /* known to be buggy in some cases */
1425 private unsafe BigInteger EvenModTwoPow (BigInteger exp)
1428 uint [] wkspace = new uint [mod.length << 1 + 1];
1430 BigInteger resultNum = new BigInteger (2, mod.length << 1 +1);
1432 uint value = exp.data [exp.length - 1];
1433 uint mask = 0x80000000;
1435 // Find the first bit of the exponent
1436 while ((value & mask) == 0)
1440 // We know that the first itr will make the val 2,
1441 // so eat one bit of the exponent
1445 uint wPos = exp.length - 1;
1448 value = exp.data [wPos];
1450 Kernel.SquarePositive (resultNum, ref wkspace);
1451 if (resultNum.length >= mod.length)
1452 BarrettReduction (resultNum);
1454 if ((value & mask) != 0) {
1456 // resultNum = (resultNum * 2) % mod
1459 fixed (uint* u = resultNum.data) {
1464 uint* uuE = u + resultNum.length;
1468 *uu = (x << 1) | carry;
1469 carry = x >> (32 - 1);
1473 // subtraction inlined because we know it is square
1474 if (carry != 0 || resultNum >= mod) {
1477 uint [] s = mod.data;
1481 if (((a += c) < c) | ((* (uu++) -= a) > ~a))
1490 } while ((mask >>= 1) > 0);
1492 } while (wPos-- > 0);
1497 private unsafe BigInteger OddModTwoPow (BigInteger exp)
1500 uint [] wkspace = new uint [mod.length << 1 + 1];
1502 BigInteger resultNum = Montgomery.ToMont ((BigInteger)2, this.mod);
1503 resultNum = new BigInteger (resultNum, mod.length << 1 +1);
1505 uint mPrime = Montgomery.Inverse (mod.data [0]);
1508 // TODO: eat small bits, the ones we can do with no modular reduction
1510 uint pos = (uint)exp.BitCount () - 2;
1513 Kernel.SquarePositive (resultNum, ref wkspace);
1514 resultNum = Montgomery.Reduce (resultNum, mod, mPrime);
1516 if (exp.TestBit (pos)) {
1518 // resultNum = (resultNum * 2) % mod
1521 fixed (uint* u = resultNum.data) {
1526 uint* uuE = u + resultNum.length;
1530 *uu = (x << 1) | carry;
1531 carry = x >> (32 - 1);
1535 // subtraction inlined because we know it is square
1536 if (carry != 0 || resultNum >= mod) {
1537 fixed (uint* s = mod.data) {
1543 if (((a += c) < c) | ((* (uu++) -= a) > ~a))
1552 } while (pos-- > 0);
1554 resultNum = Montgomery.Reduce (resultNum, mod, mPrime);
1562 /// Low level functions for the BigInteger
1564 private sealed class Kernel {
1566 #region Addition/Subtraction
1569 /// Adds two numbers with the same sign.
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)
1577 uint yMax, xMax, i = 0;
1579 // x should be bigger
1580 if (bi1.length < bi2.length) {
1592 BigInteger result = new BigInteger (Sign.Positive, xMax + 1);
1594 uint [] r = result.data;
1598 // Add common parts of both numbers
1600 sum = ((ulong)x [i]) + ((ulong)y [i]) + sum;
1603 } while (++i < yMax);
1605 // Copy remainder of longer number while carry propagation is required
1606 bool carry = (sum != 0);
1612 carry = ((r [i] = x [i] + 1) == 0);
1613 while (++i < xMax && carry);
1618 result.length = ++i;
1630 result.Normalize ();
1634 public static BigInteger Subtract (BigInteger big, BigInteger small)
1636 BigInteger result = new BigInteger (Sign.Positive, big.length);
1638 uint [] r = result.data, b = big.data, s = small.data;
1644 if (((x += c) < c) | ((r [i] = b [i] - x) > ~x))
1649 } while (++i < small.length);
1651 if (i == big.length) goto fixup;
1656 while (b [i++] == 0 && i < big.length);
1658 if (i == big.length) goto fixup;
1663 while (++i < big.length);
1667 result.Normalize ();
1671 public static void MinusEq (BigInteger big, BigInteger small)
1673 uint [] b = big.data, s = small.data;
1678 if (((x += c) < c) | ((b [i] -= x) > ~x))
1682 } while (++i < small.length);
1684 if (i == big.length) goto fixup;
1689 while (b [i++] == 0 && i < big.length);
1695 while (big.length > 0 && big.data [big.length-1] == 0) big.length--;
1698 if (big.length == 0)
1703 public static void PlusEq (BigInteger bi1, BigInteger bi2)
1706 uint yMax, xMax, i = 0;
1709 // x should be bigger
1710 if (bi1.length < bi2.length){
1723 uint [] r = bi1.data;
1727 // Add common parts of both numbers
1729 sum += ((ulong)x [i]) + ((ulong)y [i]);
1732 } while (++i < yMax);
1734 // Copy remainder of longer number while carry propagation is required
1735 bool carry = (sum != 0);
1741 carry = ((r [i] = x [i] + 1) == 0);
1742 while (++i < xMax && carry);
1753 if (flag && i < xMax - 1) {
1759 bi1.length = xMax + 1;
1768 /// Compares two BigInteger
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)
1776 // Step 1. Compare the lengths
1778 uint l1 = bi1.length, l2 = bi2.length;
1780 while (l1 > 0 && bi1.data [l1-1] == 0) l1--;
1781 while (l2 > 0 && bi2.data [l2-1] == 0) l2--;
1783 if (l1 == 0 && l2 == 0) return Sign.Zero;
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;
1791 // Step 2. Compare the bits
1796 while (pos != 0 && bi1.data [pos] == bi2.data [pos]) pos--;
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;
1813 /// Performs n / d and n % d in one operation.
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)
1826 n.data [i] = (uint)(r / d);
1834 public static uint DwordMod (BigInteger n, uint d)
1848 public static BigInteger DwordDiv (BigInteger n, uint d)
1850 BigInteger ret = new BigInteger (Sign.Positive, n.length);
1858 ret.data [i] = (uint)(r / d);
1866 public static BigInteger [] DwordDivMod (BigInteger n, uint d)
1868 BigInteger ret = new BigInteger (Sign.Positive , n.length);
1876 ret.data [i] = (uint)(r / d);
1881 BigInteger rem = (uint)r;
1883 return new BigInteger [] {ret, rem};
1890 public static BigInteger [] multiByteDivide (BigInteger bi1, BigInteger bi2)
1892 if (Kernel.Compare (bi1, bi2) == Sign.Negative)
1893 return new BigInteger [2] { 0, new BigInteger (bi1) };
1895 bi1.Normalize (); bi2.Normalize ();
1897 if (bi2.length == 1)
1898 return DwordDivMod (bi1, bi2.data [0]);
1900 uint remainderLen = bi1.length + 1;
1901 int divisorLen = (int)bi2.length + 1;
1903 uint mask = 0x80000000;
1904 uint val = bi2.data [bi2.length - 1];
1906 int resultPos = (int)bi1.length - (int)bi2.length;
1908 while (mask != 0 && (val & mask) == 0) {
1909 shift++; mask >>= 1;
1912 BigInteger quot = new BigInteger (Sign.Positive, bi1.length - bi2.length + 1);
1913 BigInteger rem = (bi1 << shift);
1915 uint [] remainder = rem.data;
1919 int j = (int)(remainderLen - bi2.length);
1920 int pos = (int)remainderLen - 1;
1922 uint firstDivisorByte = bi2.data [bi2.length-1];
1923 ulong secondDivisorByte = bi2.data [bi2.length-2];
1926 ulong dividend = ((ulong)remainder [pos] << 32) + (ulong)remainder [pos-1];
1928 ulong q_hat = dividend / (ulong)firstDivisorByte;
1929 ulong r_hat = dividend % (ulong)firstDivisorByte;
1933 if (q_hat == 0x100000000 ||
1934 (q_hat * secondDivisorByte) > ((r_hat << 32) + remainder [pos-2])) {
1936 r_hat += (ulong)firstDivisorByte;
1938 if (r_hat < 0x100000000)
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.
1953 int nPos = pos - divisorLen + 1;
1955 uint uint_q_hat = (uint)q_hat;
1957 mc += (ulong)bi2.data [dPos] * (ulong)uint_q_hat;
1958 t = remainder [nPos];
1959 remainder [nPos] -= (uint)mc;
1961 if (remainder [nPos] > t) mc++;
1963 } while (dPos < divisorLen);
1965 nPos = pos - divisorLen + 1;
1974 sum = ((ulong)remainder [nPos]) + ((ulong)bi2.data [dPos]) + sum;
1975 remainder [nPos] = (uint)sum;
1978 } while (dPos < divisorLen);
1982 quot.data [resultPos--] = (uint)uint_q_hat;
1990 BigInteger [] ret = new BigInteger [2] { quot, rem };
2003 public static BigInteger LeftShift (BigInteger bi, int n)
2005 if (n == 0) return new BigInteger (bi, bi.length + 1);
2008 n &= ((1 << 5) - 1);
2010 BigInteger ret = new BigInteger (Sign.Positive, bi.length + 1 + (uint)w);
2012 uint i = 0, l = bi.length;
2017 ret.data [i + w] = (x << n) | carry;
2018 carry = x >> (32 - n);
2021 ret.data [i + w] = carry;
2024 ret.data [i + w] = bi.data [i];
2033 public static BigInteger RightShift (BigInteger bi, int n)
2035 if (n == 0) return new BigInteger (bi);
2038 int s = n & ((1 << 5) - 1);
2040 BigInteger ret = new BigInteger (Sign.Positive, bi.length - (uint)w + 1);
2041 uint l = (uint)ret.data.Length - 1;
2048 x = bi.data [l + w];
2049 ret.data [l] = (x >> n) | carry;
2050 carry = x << (32 - n);
2054 ret.data [l] = bi.data [l + w];
2065 public static BigInteger MultiplyByDword (BigInteger n, uint f)
2067 BigInteger ret = new BigInteger (Sign.Positive, n.length + 1);
2073 c += (ulong)n.data [i] * (ulong)f;
2074 ret.data [i] = (uint)c;
2076 } while (++i < n.length);
2077 ret.data [i] = (uint)c;
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].
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].
2093 public static unsafe void Multiply (uint [] x, uint xOffset, uint xLen, uint [] y, uint yOffset, uint yLen, uint [] d, uint dOffset)
2095 fixed (uint* xx = x, yy = y, dd = d) {
2096 uint* xP = xx + xOffset,
2102 for (; xP < xE; xP++, dB++) {
2104 if (*xP == 0) continue;
2109 for (uint* yP = yB; yP < yE; yP++, dP++) {
2110 mcarry += ((ulong)*xP * (ulong)*yP) + (ulong)*dP;
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].
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].
2132 public static unsafe void MultiplyMod2p32pmod (uint [] x, int xOffset, int xLen, uint [] y, int yOffest, int yLen, uint [] d, int dOffset, int mod)
2134 fixed (uint* xx = x, yy = y, dd = d) {
2135 uint* xP = xx + xOffset,
2142 for (; xP < xE; xP++, dB++) {
2144 if (*xP == 0) continue;
2148 for (uint* yP = yB; yP < yE && dP < dE; yP++, dP++) {
2149 mcarry += ((ulong)*xP * (ulong)*yP) + (ulong)*dP;
2155 if (mcarry != 0 && dP < dE)
2161 public static unsafe void SquarePositive (BigInteger bi, ref uint [] wkSpace)
2163 uint [] t = wkSpace;
2165 uint [] d = bi.data;
2166 uint dl = bi.length;
2169 fixed (uint* dd = d, tt = t) {
2171 uint* ttE = tt + t.Length;
2173 for (uint* ttt = tt; ttt < ttE; ttt++)
2176 uint* dP = dd, tP = tt;
2178 for (uint i = 0; i < dl; i++, dP++) {
2185 uint* dP2 = dP + 1, tP2 = tP + 2*i + 1;
2187 for (uint j = i + 1; j < dl; j++, tP2++, dP2++) {
2189 mcarry += ((ulong)bi1val * (ulong)*dP2) + *tP2;
2191 *tP2 = (uint)mcarry;
2196 *tP2 = (uint)mcarry;
2199 // Double t. Inlined for speed.
2206 *tP = (x << 1) | carry;
2207 carry = x >> (32 - 1);
2210 if (carry != 0) *tP = carry;
2212 // Add in the diagnals
2216 for (uint* dE = dP + dl; (dP < dE); dP++, tP++) {
2217 ulong val = (ulong)*dP * (ulong)*dP + *tP;
2220 *(++tP) += (uint)val;
2221 if (*tP < (uint)val) {
2223 // Account for the first carry
2226 // Keep adding until no carry
2227 while ((*tP3++) == 0)
2236 while (tt [bi.length-1] == 0 && bi.length > 1) bi.length--;
2242 * Never called in BigInteger (and part of a private class)
2243 * public static bool Double (uint [] u, int l)
2249 u [i] = (x << 1) | carry;
2250 carry = x >> (32 - 1);
2253 if (carry != 0) u [l] = carry;
2259 #region Number Theory
2261 public static BigInteger gcd (BigInteger a, BigInteger b)
2268 while (x.length > 1) {
2274 if (x == 0) return g;
2276 // TODO: should we have something here if we can convert to long?
2279 // Now we can just do it with single precision. I am using the binary gcd method,
2280 // as it should be faster.
2283 uint yy = x.data [0];
2288 while (((xx | yy) & 1) == 0) {
2289 xx >>= 1; yy >>= 1; t++;
2292 while ((xx & 1) == 0) xx >>= 1;
2293 while ((yy & 1) == 0) yy >>= 1;
2295 xx = (xx - yy) >> 1;
2297 yy = (yy - xx) >> 1;
2303 public static uint modInverse (BigInteger bi, uint modulus)
2305 uint a = modulus, b = bi % modulus;
2306 uint p0 = 0, p1 = 1;
2326 public static BigInteger modInverse (BigInteger bi, BigInteger modulus)
2328 if (modulus.length == 1) return modInverse (bi, modulus.data [0]);
2330 BigInteger [] p = { 0, 1 };
2331 BigInteger [] q = new BigInteger [2]; // quotients
2332 BigInteger [] r = { 0, 0 }; // remainders
2336 BigInteger a = modulus;
2339 ModulusRing mr = new ModulusRing (modulus);
2345 BigInteger pval = mr.Difference (p [0], p [1] * q [0]);
2346 p [0] = p [1]; p [1] = pval;
2349 BigInteger [] divret = multiByteDivide (a, b);
2351 q [0] = q [1]; q [1] = divret [0];
2352 r [0] = r [1]; r [1] = divret [1];
2360 throw (new ArithmeticException ("No inverse!"));
2362 return mr.Difference (p [0], p [1] * q [0]);