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 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)
865 if (o == null) return false;
866 if (o is int) return (int)o >= 0 && this == (uint)o;
868 return Kernel.Compare (this, (BigInteger)o) == 0;
873 #region Number Theory
875 public BigInteger GCD (BigInteger bi)
877 return Kernel.gcd (this, bi);
880 public BigInteger ModInverse (BigInteger modulus)
882 return Kernel.modInverse (this, modulus);
885 public BigInteger ModPow (BigInteger exp, BigInteger n)
887 ModulusRing mr = new ModulusRing (n);
888 return mr.Pow (this, exp);
893 #region Prime Testing
895 public bool IsProbablePrime ()
897 if (this < smallPrimes [smallPrimes.Length - 1]) {
898 for (int p = 0; p < smallPrimes.Length; p++) {
899 if (this == smallPrimes [p])
904 for (int p = 0; p < smallPrimes.Length; p++) {
905 if (this % smallPrimes [p] == 0)
909 return PrimalityTests.RabinMillerTest (this, Prime.ConfidenceFactor.Medium);
914 #region Prime Number Generation
917 /// Generates the smallest prime >= bi
919 /// <param name="bi">A BigInteger</param>
920 /// <returns>The smallest prime >= bi. More mathematically, if bi is prime: bi, else Prime [PrimePi [bi] + 1].</returns>
921 public static BigInteger NextHighestPrime (BigInteger bi)
923 NextPrimeFinder npf = new NextPrimeFinder ();
924 return npf.GenerateNewPrime (0, bi);
927 public static BigInteger GeneratePseudoPrime (int bits)
929 SequentialSearchPrimeGeneratorBase sspg = new SequentialSearchPrimeGeneratorBase ();
930 return sspg.GenerateNewPrime (bits);
934 /// Increments this by two
942 // If there was no carry, nothing to do
945 // Account for the first carry
948 // Keep adding until no carry
949 while (data [i++] == 0x0)
952 // See if we increased the data length
953 if (length == (uint)i)
965 sealed class ModulusRing {
967 BigInteger mod, constant;
969 public ModulusRing (BigInteger modulus)
973 // calculate constant = b^ (2k) / m
974 uint i = mod.length << 1;
976 constant = new BigInteger (Sign.Positive, i + 1);
977 constant.data [i] = 0x00000001;
979 constant = constant / mod;
982 public void BarrettReduction (BigInteger x)
989 // x < mod, so nothing to do.
990 if (x.length < k) return;
997 if (x.data.Length < x.length) throw new IndexOutOfRangeException ("x out of range");
1000 // q2 = q1 * constant
1001 // q3 = q2 / b^ (k+1), Needs to be accessed with an offset of kPlusOne
1003 // TODO: We should the method in HAC p 604 to do this (14.45)
1004 q3 = new BigInteger (Sign.Positive, x.length - kMinusOne + constant.length);
1005 Kernel.Multiply (x.data, kMinusOne, x.length - kMinusOne, constant.data, 0, constant.length, q3.data, 0);
1007 // r1 = x mod b^ (k+1)
1008 // i.e. keep the lowest (k+1) words
1010 uint lengthToCopy = (x.length > kPlusOne) ? kPlusOne : x.length;
1012 x.length = lengthToCopy;
1015 // r2 = (q3 * n) mod b^ (k+1)
1016 // partial multiplication of q3 and n
1018 BigInteger r2 = new BigInteger (Sign.Positive, kPlusOne);
1019 Kernel.MultiplyMod2p32pmod (q3.data, (int)kPlusOne, (int)q3.length - (int)kPlusOne, n.data, 0, (int)n.length, r2.data, 0, (int)kPlusOne);
1024 Kernel.MinusEq (x, r2);
1026 BigInteger val = new BigInteger (Sign.Positive, kPlusOne + 1);
1027 val.data [kPlusOne] = 0x00000001;
1029 Kernel.MinusEq (val, r2);
1030 Kernel.PlusEq (x, val);
1034 Kernel.MinusEq (x, n);
1037 public BigInteger Multiply (BigInteger a, BigInteger b)
1039 if (a == 0 || b == 0) return 0;
1041 if (a.length >= mod.length << 1)
1044 if (b.length >= mod.length << 1)
1047 if (a.length >= mod.length)
1048 BarrettReduction (a);
1050 if (b.length >= mod.length)
1051 BarrettReduction (b);
1053 BigInteger ret = new BigInteger (a * b);
1054 BarrettReduction (ret);
1059 public BigInteger Difference (BigInteger a, BigInteger b)
1061 Sign cmp = Kernel.Compare (a, b);
1068 diff = a - b; break;
1070 diff = b - a; break;
1072 throw new Exception ();
1076 if (diff.length >= mod.length << 1)
1079 BarrettReduction (diff);
1081 if (cmp == Sign.Negative)
1086 public BigInteger Pow (BigInteger b, BigInteger exp)
1088 if ((mod.data [0] & 1) == 1) return OddPow (b, exp);
1089 else return EvenPow (b, exp);
1092 public BigInteger EvenPow (BigInteger b, BigInteger exp)
1094 BigInteger resultNum = new BigInteger ((BigInteger)1, mod.length << 1);
1095 BigInteger tempNum = new BigInteger (b % mod, mod.length << 1); // ensures (tempNum * tempNum) < b^ (2k)
1097 uint totalBits = (uint)exp.BitCount ();
1099 uint [] wkspace = new uint [mod.length << 1];
1101 // perform squaring and multiply exponentiation
1102 for (uint pos = 0; pos < totalBits; pos++) {
1103 if (exp.TestBit (pos)) {
1105 Array.Clear (wkspace, 0, wkspace.Length);
1106 Kernel.Multiply (resultNum.data, 0, resultNum.length, tempNum.data, 0, tempNum.length, wkspace, 0);
1107 resultNum.length += tempNum.length;
1108 uint [] t = wkspace;
1109 wkspace = resultNum.data;
1112 BarrettReduction (resultNum);
1115 Kernel.SquarePositive (tempNum, ref wkspace);
1116 BarrettReduction (tempNum);
1126 private BigInteger OddPow (BigInteger b, BigInteger exp)
1128 BigInteger resultNum = new BigInteger (Montgomery.ToMont (1, mod), mod.length << 1);
1129 BigInteger tempNum = new BigInteger (Montgomery.ToMont (b, mod), mod.length << 1); // ensures (tempNum * tempNum) < b^ (2k)
1130 uint mPrime = Montgomery.Inverse (mod.data [0]);
1131 uint totalBits = (uint)exp.BitCount ();
1133 uint [] wkspace = new uint [mod.length << 1];
1135 // perform squaring and multiply exponentiation
1136 for (uint pos = 0; pos < totalBits; pos++) {
1137 if (exp.TestBit (pos)) {
1139 Array.Clear (wkspace, 0, wkspace.Length);
1140 Kernel.Multiply (resultNum.data, 0, resultNum.length, tempNum.data, 0, tempNum.length, wkspace, 0);
1141 resultNum.length += tempNum.length;
1142 uint [] t = wkspace;
1143 wkspace = resultNum.data;
1146 Montgomery.Reduce (resultNum, mod, mPrime);
1149 Kernel.SquarePositive (tempNum, ref wkspace);
1150 Montgomery.Reduce (tempNum, mod, mPrime);
1153 Montgomery.Reduce (resultNum, mod, mPrime);
1157 #region Pow Small Base
1159 // TODO: Make tests for this, not really needed b/c prime stuff
1160 // checks it, but still would be nice
1162 [CLSCompliant (false)]
1164 public BigInteger Pow (uint b, BigInteger exp)
1167 if ((mod.data [0] & 1) == 1)
1168 return OddPow (b, exp);
1170 return EvenPow (b, exp);
1171 /* buggy in some cases (like the well tested primes)
1173 if ((mod.data [0] & 1) == 1)
1174 return OddModTwoPow (exp);
1176 return EvenModTwoPow (exp);
1180 private unsafe BigInteger OddPow (uint b, BigInteger exp)
1183 uint [] wkspace = new uint [mod.length << 1 + 1];
1185 BigInteger resultNum = Montgomery.ToMont ((BigInteger)b, this.mod);
1186 resultNum = new BigInteger (resultNum, mod.length << 1 +1);
1188 uint mPrime = Montgomery.Inverse (mod.data [0]);
1190 uint pos = (uint)exp.BitCount () - 2;
1193 // We know that the first itr will make the val b
1200 Kernel.SquarePositive (resultNum, ref wkspace);
1201 resultNum = Montgomery.Reduce (resultNum, mod, mPrime);
1203 if (exp.TestBit (pos)) {
1209 // TODO: Is Unsafe really speeding things up?
1210 fixed (uint* u = resultNum.data) {
1216 mc += (ulong)u [i] * (ulong)b;
1219 } while (++i < resultNum.length);
1221 if (resultNum.length < mod.length) {
1225 while (resultNum >= mod)
1226 Kernel.MinusEq (resultNum, mod);
1228 } else if (mc != 0) {
1231 // First, we estimate the quotient by dividing
1232 // the first part of each of the numbers. Then
1233 // we correct this, if necessary, with a subtraction.
1238 // We would rather have this estimate overshoot,
1239 // so we add one to the divisor
1241 if (mod.data [mod.length - 1] < UInt32.MaxValue) {
1242 divEstimate = (uint) ((((ulong)cc << 32) | (ulong) u [i -1]) /
1243 (mod.data [mod.length-1] + 1));
1246 // guess but don't divide by 0
1247 divEstimate = (uint) ((((ulong)cc << 32) | (ulong) u [i -1]) /
1248 (mod.data [mod.length-1]));
1256 mc += (ulong)mod.data [i] * (ulong)divEstimate;
1260 if (u [i] > t) mc++;
1262 } while (i < resultNum.length);
1268 uint [] s = mod.data;
1271 if (((a += sc) < sc) | ((u [j] -= a) > ~a)) sc = 1;
1274 } while (j < resultNum.length);
1277 while (resultNum >= mod)
1278 Kernel.MinusEq (resultNum, mod);
1280 while (resultNum >= mod)
1281 Kernel.MinusEq (resultNum, mod);
1285 } while (pos-- > 0);
1287 resultNum = Montgomery.Reduce (resultNum, mod, mPrime);
1292 private unsafe BigInteger EvenPow (uint b, BigInteger exp)
1295 uint [] wkspace = new uint [mod.length << 1 + 1];
1296 BigInteger resultNum = new BigInteger ((BigInteger)b, mod.length << 1 + 1);
1298 uint pos = (uint)exp.BitCount () - 2;
1301 // We know that the first itr will make the val b
1308 Kernel.SquarePositive (resultNum, ref wkspace);
1309 if (!(resultNum.length < mod.length))
1310 BarrettReduction (resultNum);
1312 if (exp.TestBit (pos)) {
1318 // TODO: Is Unsafe really speeding things up?
1319 fixed (uint* u = resultNum.data) {
1325 mc += (ulong)u [i] * (ulong)b;
1328 } while (++i < resultNum.length);
1330 if (resultNum.length < mod.length) {
1334 while (resultNum >= mod)
1335 Kernel.MinusEq (resultNum, mod);
1337 } else if (mc != 0) {
1340 // First, we estimate the quotient by dividing
1341 // the first part of each of the numbers. Then
1342 // we correct this, if necessary, with a subtraction.
1347 // We would rather have this estimate overshoot,
1348 // so we add one to the divisor
1349 uint divEstimate = (uint) ((((ulong)cc << 32) | (ulong) u [i -1]) /
1350 (mod.data [mod.length-1] + 1));
1357 mc += (ulong)mod.data [i] * (ulong)divEstimate;
1361 if (u [i] > t) mc++;
1363 } while (i < resultNum.length);
1369 uint [] s = mod.data;
1372 if (((a += sc) < sc) | ((u [j] -= a) > ~a)) sc = 1;
1375 } while (j < resultNum.length);
1378 while (resultNum >= mod)
1379 Kernel.MinusEq (resultNum, mod);
1381 while (resultNum >= mod)
1382 Kernel.MinusEq (resultNum, mod);
1386 } while (pos-- > 0);
1391 /* known to be buggy in some cases
1392 private unsafe BigInteger EvenModTwoPow (BigInteger exp)
1395 uint [] wkspace = new uint [mod.length << 1 + 1];
1397 BigInteger resultNum = new BigInteger (2, mod.length << 1 +1);
1399 uint value = exp.data [exp.length - 1];
1400 uint mask = 0x80000000;
1402 // Find the first bit of the exponent
1403 while ((value & mask) == 0)
1407 // We know that the first itr will make the val 2,
1408 // so eat one bit of the exponent
1412 uint wPos = exp.length - 1;
1415 value = exp.data [wPos];
1417 Kernel.SquarePositive (resultNum, ref wkspace);
1418 if (resultNum.length >= mod.length)
1419 BarrettReduction (resultNum);
1421 if ((value & mask) != 0) {
1423 // resultNum = (resultNum * 2) % mod
1426 fixed (uint* u = resultNum.data) {
1431 uint* uuE = u + resultNum.length;
1435 *uu = (x << 1) | carry;
1436 carry = x >> (32 - 1);
1440 // subtraction inlined because we know it is square
1441 if (carry != 0 || resultNum >= mod) {
1444 uint [] s = mod.data;
1448 if (((a += c) < c) | ((* (uu++) -= a) > ~a))
1457 } while ((mask >>= 1) > 0);
1459 } while (wPos-- > 0);
1464 private unsafe BigInteger OddModTwoPow (BigInteger exp)
1467 uint [] wkspace = new uint [mod.length << 1 + 1];
1469 BigInteger resultNum = Montgomery.ToMont ((BigInteger)2, this.mod);
1470 resultNum = new BigInteger (resultNum, mod.length << 1 +1);
1472 uint mPrime = Montgomery.Inverse (mod.data [0]);
1475 // TODO: eat small bits, the ones we can do with no modular reduction
1477 uint pos = (uint)exp.BitCount () - 2;
1480 Kernel.SquarePositive (resultNum, ref wkspace);
1481 resultNum = Montgomery.Reduce (resultNum, mod, mPrime);
1483 if (exp.TestBit (pos)) {
1485 // resultNum = (resultNum * 2) % mod
1488 fixed (uint* u = resultNum.data) {
1493 uint* uuE = u + resultNum.length;
1497 *uu = (x << 1) | carry;
1498 carry = x >> (32 - 1);
1502 // subtraction inlined because we know it is square
1503 if (carry != 0 || resultNum >= mod) {
1504 fixed (uint* s = mod.data) {
1510 if (((a += c) < c) | ((* (uu++) -= a) > ~a))
1519 } while (pos-- > 0);
1521 resultNum = Montgomery.Reduce (resultNum, mod, mPrime);
1528 internal sealed class Montgomery {
1530 private Montgomery ()
1534 public static uint Inverse (uint n)
1538 while ((z = n * y) != 1)
1544 public static BigInteger ToMont (BigInteger n, BigInteger m)
1546 n.Normalize (); m.Normalize ();
1548 n <<= (int)m.length * 32;
1553 public static unsafe BigInteger Reduce (BigInteger n, BigInteger m, uint mPrime)
1556 fixed (uint* a = A.data, mm = m.data) {
1557 for (uint i = 0; i < m.length; i++) {
1558 // The mod here is taken care of by the CPU,
1559 // since the multiply will overflow.
1560 uint u_i = a [0] * mPrime /* % 2^32 */;
1567 // mP = Position in mod
1568 // aSP = the source of bits from a
1569 // aDP = destination for bits
1570 uint* mP = mm, aSP = a, aDP = a;
1572 ulong c = (ulong)u_i * ((ulong)*(mP++)) + *(aSP++);
1577 for (; j < m.length; j++) {
1578 c += (ulong)u_i * (ulong)*(mP++) + *(aSP++);
1583 // Account for carry
1584 // TODO: use a better loop here, we dont need the ulong stuff
1585 for (; j < A.length; j++) {
1589 if (c == 0) {j++; break;}
1592 for (; j < A.length; j++) {
1593 *(aDP++) = *(aSP++);
1599 while (A.length > 1 && a [A.length-1] == 0) A.length--;
1602 if (A >= m) Kernel.MinusEq (A, m);
1607 public static BigInteger Reduce (BigInteger n, BigInteger m)
1609 return Reduce (n, m, Inverse (m.data [0]));
1615 /// Low level functions for the BigInteger
1617 private sealed class Kernel {
1619 #region Addition/Subtraction
1622 /// Adds two numbers with the same sign.
1624 /// <param name="bi1">A BigInteger</param>
1625 /// <param name="bi2">A BigInteger</param>
1626 /// <returns>bi1 + bi2</returns>
1627 public static BigInteger AddSameSign (BigInteger bi1, BigInteger bi2)
1630 uint yMax, xMax, i = 0;
1632 // x should be bigger
1633 if (bi1.length < bi2.length) {
1645 BigInteger result = new BigInteger (Sign.Positive, xMax + 1);
1647 uint [] r = result.data;
1651 // Add common parts of both numbers
1653 sum = ((ulong)x [i]) + ((ulong)y [i]) + sum;
1656 } while (++i < yMax);
1658 // Copy remainder of longer number while carry propagation is required
1659 bool carry = (sum != 0);
1665 carry = ((r [i] = x [i] + 1) == 0);
1666 while (++i < xMax && carry);
1671 result.length = ++i;
1683 result.Normalize ();
1687 public static BigInteger Subtract (BigInteger big, BigInteger small)
1689 BigInteger result = new BigInteger (Sign.Positive, big.length);
1691 uint [] r = result.data, b = big.data, s = small.data;
1697 if (((x += c) < c) | ((r [i] = b [i] - x) > ~x))
1702 } while (++i < small.length);
1704 if (i == big.length) goto fixup;
1709 while (b [i++] == 0 && i < big.length);
1711 if (i == big.length) goto fixup;
1716 while (++i < big.length);
1720 result.Normalize ();
1724 public static void MinusEq (BigInteger big, BigInteger small)
1726 uint [] b = big.data, s = small.data;
1731 if (((x += c) < c) | ((b [i] -= x) > ~x))
1735 } while (++i < small.length);
1737 if (i == big.length) goto fixup;
1742 while (b [i++] == 0 && i < big.length);
1748 while (big.length > 0 && big.data [big.length-1] == 0) big.length--;
1751 if (big.length == 0)
1756 public static void PlusEq (BigInteger bi1, BigInteger bi2)
1759 uint yMax, xMax, i = 0;
1762 // x should be bigger
1763 if (bi1.length < bi2.length){
1776 uint [] r = bi1.data;
1780 // Add common parts of both numbers
1782 sum += ((ulong)x [i]) + ((ulong)y [i]);
1785 } while (++i < yMax);
1787 // Copy remainder of longer number while carry propagation is required
1788 bool carry = (sum != 0);
1794 carry = ((r [i] = x [i] + 1) == 0);
1795 while (++i < xMax && carry);
1806 if (flag && i < xMax - 1) {
1812 bi1.length = xMax + 1;
1821 /// Compares two BigInteger
1823 /// <param name="bi1">A BigInteger</param>
1824 /// <param name="bi2">A BigInteger</param>
1825 /// <returns>The sign of bi1 - bi2</returns>
1826 public static Sign Compare (BigInteger bi1, BigInteger bi2)
1829 // Step 1. Compare the lengths
1831 uint l1 = bi1.length, l2 = bi2.length;
1833 while (l1 > 0 && bi1.data [l1-1] == 0) l1--;
1834 while (l2 > 0 && bi2.data [l2-1] == 0) l2--;
1836 if (l1 == 0 && l2 == 0) return Sign.Zero;
1838 // bi1 len < bi2 len
1839 if (l1 < l2) return Sign.Negative;
1840 // bi1 len > bi2 len
1841 else if (l1 > l2) return Sign.Positive;
1844 // Step 2. Compare the bits
1849 while (pos != 0 && bi1.data [pos] == bi2.data [pos]) pos--;
1851 if (bi1.data [pos] < bi2.data [pos])
1852 return Sign.Negative;
1853 else if (bi1.data [pos] > bi2.data [pos])
1854 return Sign.Positive;
1866 /// Performs n / d and n % d in one operation.
1868 /// <param name="n">A BigInteger, upon exit this will hold n / d</param>
1869 /// <param name="d">The divisor</param>
1870 /// <returns>n % d</returns>
1871 public static uint SingleByteDivideInPlace (BigInteger n, uint d)
1879 n.data [i] = (uint)(r / d);
1887 public static uint DwordMod (BigInteger n, uint d)
1901 public static BigInteger DwordDiv (BigInteger n, uint d)
1903 BigInteger ret = new BigInteger (Sign.Positive, n.length);
1911 ret.data [i] = (uint)(r / d);
1919 public static BigInteger [] DwordDivMod (BigInteger n, uint d)
1921 BigInteger ret = new BigInteger (Sign.Positive , n.length);
1929 ret.data [i] = (uint)(r / d);
1934 BigInteger rem = (uint)r;
1936 return new BigInteger [] {ret, rem};
1943 public static BigInteger [] multiByteDivide (BigInteger bi1, BigInteger bi2)
1945 if (Kernel.Compare (bi1, bi2) == Sign.Negative)
1946 return new BigInteger [2] { 0, new BigInteger (bi1) };
1948 bi1.Normalize (); bi2.Normalize ();
1950 if (bi2.length == 1)
1951 return DwordDivMod (bi1, bi2.data [0]);
1953 uint remainderLen = bi1.length + 1;
1954 int divisorLen = (int)bi2.length + 1;
1956 uint mask = 0x80000000;
1957 uint val = bi2.data [bi2.length - 1];
1959 int resultPos = (int)bi1.length - (int)bi2.length;
1961 while (mask != 0 && (val & mask) == 0) {
1962 shift++; mask >>= 1;
1965 BigInteger quot = new BigInteger (Sign.Positive, bi1.length - bi2.length + 1);
1966 BigInteger rem = (bi1 << shift);
1968 uint [] remainder = rem.data;
1972 int j = (int)(remainderLen - bi2.length);
1973 int pos = (int)remainderLen - 1;
1975 uint firstDivisorByte = bi2.data [bi2.length-1];
1976 ulong secondDivisorByte = bi2.data [bi2.length-2];
1979 ulong dividend = ((ulong)remainder [pos] << 32) + (ulong)remainder [pos-1];
1981 ulong q_hat = dividend / (ulong)firstDivisorByte;
1982 ulong r_hat = dividend % (ulong)firstDivisorByte;
1986 if (q_hat == 0x100000000 ||
1987 (q_hat * secondDivisorByte) > ((r_hat << 32) + remainder [pos-2])) {
1989 r_hat += (ulong)firstDivisorByte;
1991 if (r_hat < 0x100000000)
1998 // At this point, q_hat is either exact, or one too large
1999 // (more likely to be exact) so, we attempt to multiply the
2000 // divisor by q_hat, if we get a borrow, we just subtract
2001 // one from q_hat and add the divisor back.
2006 int nPos = pos - divisorLen + 1;
2008 uint uint_q_hat = (uint)q_hat;
2010 mc += (ulong)bi2.data [dPos] * (ulong)uint_q_hat;
2011 t = remainder [nPos];
2012 remainder [nPos] -= (uint)mc;
2014 if (remainder [nPos] > t) mc++;
2016 } while (dPos < divisorLen);
2018 nPos = pos - divisorLen + 1;
2027 sum = ((ulong)remainder [nPos]) + ((ulong)bi2.data [dPos]) + sum;
2028 remainder [nPos] = (uint)sum;
2031 } while (dPos < divisorLen);
2035 quot.data [resultPos--] = (uint)uint_q_hat;
2043 BigInteger [] ret = new BigInteger [2] { quot, rem };
2056 public static BigInteger LeftShift (BigInteger bi, int n)
2058 if (n == 0) return new BigInteger (bi, bi.length + 1);
2061 n &= ((1 << 5) - 1);
2063 BigInteger ret = new BigInteger (Sign.Positive, bi.length + 1 + (uint)w);
2065 uint i = 0, l = bi.length;
2070 ret.data [i + w] = (x << n) | carry;
2071 carry = x >> (32 - n);
2074 ret.data [i + w] = carry;
2077 ret.data [i + w] = bi.data [i];
2086 public static BigInteger RightShift (BigInteger bi, int n)
2088 if (n == 0) return new BigInteger (bi);
2091 int s = n & ((1 << 5) - 1);
2093 BigInteger ret = new BigInteger (Sign.Positive, bi.length - (uint)w + 1);
2094 uint l = (uint)ret.data.Length - 1;
2101 x = bi.data [l + w];
2102 ret.data [l] = (x >> n) | carry;
2103 carry = x << (32 - n);
2107 ret.data [l] = bi.data [l + w];
2118 public static BigInteger MultiplyByDword (BigInteger n, uint f)
2120 BigInteger ret = new BigInteger (Sign.Positive, n.length + 1);
2126 c += (ulong)n.data [i] * (ulong)f;
2127 ret.data [i] = (uint)c;
2129 } while (++i < n.length);
2130 ret.data [i] = (uint)c;
2137 /// Multiplies the data in x [xOffset:xOffset+xLen] by
2138 /// y [yOffset:yOffset+yLen] and puts it into
2139 /// d [dOffset:dOffset+xLen+yLen].
2142 /// This code is unsafe! It is the caller's responsibility to make
2143 /// sure that it is safe to access x [xOffset:xOffset+xLen],
2144 /// y [yOffset:yOffset+yLen], and d [dOffset:dOffset+xLen+yLen].
2146 public static unsafe void Multiply (uint [] x, uint xOffset, uint xLen, uint [] y, uint yOffset, uint yLen, uint [] d, uint dOffset)
2148 fixed (uint* xx = x, yy = y, dd = d) {
2149 uint* xP = xx + xOffset,
2155 for (; xP < xE; xP++, dB++) {
2157 if (*xP == 0) continue;
2162 for (uint* yP = yB; yP < yE; yP++, dP++) {
2163 mcarry += ((ulong)*xP * (ulong)*yP) + (ulong)*dP;
2176 /// Multiplies the data in x [xOffset:xOffset+xLen] by
2177 /// y [yOffset:yOffset+yLen] and puts the low mod words into
2178 /// d [dOffset:dOffset+mod].
2181 /// This code is unsafe! It is the caller's responsibility to make
2182 /// sure that it is safe to access x [xOffset:xOffset+xLen],
2183 /// y [yOffset:yOffset+yLen], and d [dOffset:dOffset+mod].
2185 public static unsafe void MultiplyMod2p32pmod (uint [] x, int xOffset, int xLen, uint [] y, int yOffest, int yLen, uint [] d, int dOffset, int mod)
2187 fixed (uint* xx = x, yy = y, dd = d) {
2188 uint* xP = xx + xOffset,
2195 for (; xP < xE; xP++, dB++) {
2197 if (*xP == 0) continue;
2201 for (uint* yP = yB; yP < yE && dP < dE; yP++, dP++) {
2202 mcarry += ((ulong)*xP * (ulong)*yP) + (ulong)*dP;
2208 if (mcarry != 0 && dP < dE)
2214 public static unsafe void SquarePositive (BigInteger bi, ref uint [] wkSpace)
2216 uint [] t = wkSpace;
2218 uint [] d = bi.data;
2219 uint dl = bi.length;
2222 fixed (uint* dd = d, tt = t) {
2224 uint* ttE = tt + t.Length;
2226 for (uint* ttt = tt; ttt < ttE; ttt++)
2229 uint* dP = dd, tP = tt;
2231 for (uint i = 0; i < dl; i++, dP++) {
2238 uint* dP2 = dP + 1, tP2 = tP + 2*i + 1;
2240 for (uint j = i + 1; j < dl; j++, tP2++, dP2++) {
2242 mcarry += ((ulong)bi1val * (ulong)*dP2) + *tP2;
2244 *tP2 = (uint)mcarry;
2249 *tP2 = (uint)mcarry;
2252 // Double t. Inlined for speed.
2259 *tP = (x << 1) | carry;
2260 carry = x >> (32 - 1);
2263 if (carry != 0) *tP = carry;
2265 // Add in the diagnals
2269 for (uint* dE = dP + dl; (dP < dE); dP++, tP++) {
2270 ulong val = (ulong)*dP * (ulong)*dP + *tP;
2273 *(++tP) += (uint)val;
2274 if (*tP < (uint)val) {
2276 // Account for the first carry
2279 // Keep adding until no carry
2280 while ((*tP3++) == 0)
2289 while (tt [bi.length-1] == 0 && bi.length > 1) bi.length--;
2295 * Never called in BigInteger (and part of a private class)
2296 * public static bool Double (uint [] u, int l)
2302 u [i] = (x << 1) | carry;
2303 carry = x >> (32 - 1);
2306 if (carry != 0) u [l] = carry;
2312 #region Number Theory
2314 public static BigInteger gcd (BigInteger a, BigInteger b)
2321 while (x.length > 1) {
2327 if (x == 0) return g;
2329 // TODO: should we have something here if we can convert to long?
2332 // Now we can just do it with single precision. I am using the binary gcd method,
2333 // as it should be faster.
2336 uint yy = x.data [0];
2341 while (((xx | yy) & 1) == 0) {
2342 xx >>= 1; yy >>= 1; t++;
2345 while ((xx & 1) == 0) xx >>= 1;
2346 while ((yy & 1) == 0) yy >>= 1;
2348 xx = (xx - yy) >> 1;
2350 yy = (yy - xx) >> 1;
2356 public static uint modInverse (BigInteger bi, uint modulus)
2358 uint a = modulus, b = bi % modulus;
2359 uint p0 = 0, p1 = 1;
2379 public static BigInteger modInverse (BigInteger bi, BigInteger modulus)
2381 if (modulus.length == 1) return modInverse (bi, modulus.data [0]);
2383 BigInteger [] p = { 0, 1 };
2384 BigInteger [] q = new BigInteger [2]; // quotients
2385 BigInteger [] r = { 0, 0 }; // remainders
2389 BigInteger a = modulus;
2392 ModulusRing mr = new ModulusRing (modulus);
2398 BigInteger pval = mr.Difference (p [0], p [1] * q [0]);
2399 p [0] = p [1]; p [1] = pval;
2402 BigInteger [] divret = multiByteDivide (a, b);
2404 q [0] = q [1]; q [1] = divret [0];
2405 r [0] = r [1]; r [1] = divret [1];
2413 throw (new ArithmeticException ("No inverse!"));
2415 return mr.Difference (p [0], p [1] * q [0]);