?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Merge changes, that were only made in corlib, back into Mono.Security
[mono-project.git] / mcs / class / Mono.Security / Mono.Math / BigInteger.cs
blob0356d0a166487af7ffae6befcd7c59b362d180d0
1 //
2 // BigInteger.cs - Big Integer implementation
3 //
4 // Authors:
5 // Ben Maurer
6 // Chew Keong TAN
7 // Sebastien Pouliot <sebastien@ximian.com>
8 // Pieter Philippaerts <Pieter@mentalis.org>
9 //
10 // Copyright (c) 2003 Ben Maurer
11 // All rights reserved
12 //
13 // Copyright (c) 2002 Chew Keong TAN
14 // All rights reserved.
15 //
16 // Copyright (C) 2004, 2007 Novell, Inc (http://www.novell.com)
17 //
18 // Permission is hereby granted, free of charge, to any person obtaining
19 // a copy of this software and associated documentation files (the
20 // "Software"), to deal in the Software without restriction, including
21 // without limitation the rights to use, copy, modify, merge, publish,
22 // distribute, sublicense, and/or sell copies of the Software, and to
23 // permit persons to whom the Software is furnished to do so, subject to
24 // the following conditions:
25 //
26 // The above copyright notice and this permission notice shall be
27 // included in all copies or substantial portions of the Software.
28 //
29 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
30 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
31 // MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
32 // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
33 // LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
34 // OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
35 // WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
36 //
37
38 using System;
39 using System.Security.Cryptography;
40 using Mono.Math.Prime.Generator;
41 using Mono.Math.Prime;
42
43 namespace Mono.Math {
44
45 #if INSIDE_CORLIB
46 internal
47 #else
48 public
49 #endif
50 class BigInteger {
51
52 #region Data Storage
53
54 /// <summary>
55 /// The Length of this BigInteger
56 /// </summary>
57 uint length = 1;
58
59 /// <summary>
60 /// The data for this BigInteger
61 /// </summary>
62 uint [] data;
63
64 #endregion
65
66 #region Constants
67
68 /// <summary>
69 /// Default length of a BigInteger in bytes
70 /// </summary>
71 const uint DEFAULT_LEN = 20;
72
73 /// <summary>
74 /// Table of primes below 2000.
75 /// </summary>
76 /// <remarks>
77 /// <para>
78 /// This table was generated using Mathematica 4.1 using the following function:
79 /// </para>
80 /// <para>
81 /// <code>
82 /// PrimeTable [x_] := Prime [Range [1, PrimePi [x]]]
83 /// PrimeTable [6000]
84 /// </code>
85 /// </para>
86 /// </remarks>
87 internal static readonly uint [] smallPrimes = {
88 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
89 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
90 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233,
91 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317,
92 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419,
93 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503,
94 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607,
95 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701,
96 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811,
97 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911,
98 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997,
99
100 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087,
101 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181,
102 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279,
103 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373,
104 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471,
105 1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559,
106 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637,
107 1657, 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747,
108 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861, 1867,
109 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973,
110 1979, 1987, 1993, 1997, 1999,
111
112 2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087, 2089,
113 2099, 2111, 2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207,
114 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, 2293, 2297,
115 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377, 2381, 2383, 2389,
116 2393, 2399, 2411, 2417, 2423, 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503,
117 2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, 2621,
118 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707,
119 2711, 2713, 2719, 2729, 2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791, 2797,
120 2801, 2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903,
121 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999,
122
123 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, 3083, 3089, 3109,
124 3119, 3121, 3137, 3163, 3167, 3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221,
125 3229, 3251, 3253, 3257, 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329,
126 3331, 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, 3433, 3449,
127 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533, 3539,
128 3541, 3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631,
129 3637, 3643, 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, 3733,
130 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, 3823, 3833, 3847, 3851,
131 3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923, 3929, 3931, 3943,
132 3947, 3967, 3989,
133
134 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, 4073, 4079, 4091,
135 4093, 4099, 4111, 4127, 4129, 4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211,
136 4217, 4219, 4229, 4231, 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289,
137 4297, 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, 4421, 4423,
138 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, 4507, 4513, 4517, 4519, 4523,
139 4547, 4549, 4561, 4567, 4583, 4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649,
140 4651, 4657, 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, 4759,
141 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, 4861, 4871, 4877, 4889,
142 4903, 4909, 4919, 4931, 4933, 4937, 4943, 4951, 4957, 4967, 4969, 4973, 4987,
143 4993, 4999,
144
145 5003, 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, 5099, 5101,
146 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179, 5189, 5197, 5209, 5227, 5231,
147 5233, 5237, 5261, 5273, 5279, 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351,
148 5381, 5387, 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, 5449,
149 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, 5527, 5531, 5557, 5563,
150 5569, 5573, 5581, 5591, 5623, 5639, 5641, 5647, 5651, 5653, 5657, 5659, 5669,
151 5683, 5689, 5693, 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791,
152 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857, 5861, 5867, 5869,
153 5879, 5881, 5897, 5903, 5923, 5927, 5939, 5953, 5981, 5987
154 };
155
156 public enum Sign : int {
157 Negative = -1,
158 Zero = 0,
159 Positive = 1
160 };
161
162 #region Exception Messages
163 const string WouldReturnNegVal = "Operation would return a negative value";
164 #endregion
165
166 #endregion
167
168 #region Constructors
169
170 public BigInteger ()
171 {
172 data = new uint [DEFAULT_LEN];
173 this.length = DEFAULT_LEN;
174 }
175
176 #if !INSIDE_CORLIB
177 [CLSCompliant (false)]
178 #endif
179 public BigInteger (Sign sign, uint len)
180 {
181 this.data = new uint [len];
182 this.length = len;
183 }
184
185 public BigInteger (BigInteger bi)
186 {
187 this.data = (uint [])bi.data.Clone ();
188 this.length = bi.length;
189 }
190
191 #if !INSIDE_CORLIB
192 [CLSCompliant (false)]
193 #endif
194 public BigInteger (BigInteger bi, uint len)
195 {
196
197 this.data = new uint [len];
198
199 for (uint i = 0; i < bi.length; i++)
200 this.data [i] = bi.data [i];
201
202 this.length = bi.length;
203 }
204
205 #endregion
206
207 #region Conversions
208
209 public BigInteger (byte [] inData)
210 {
211 if (inData.Length == 0)
212 inData = new byte [1];
213 length = (uint)inData.Length >> 2;
214 int leftOver = inData.Length & 0x3;
215
216 // length not multiples of 4
217 if (leftOver != 0) length++;
218
219 data = new uint [length];
220
221 for (int i = inData.Length - 1, j = 0; i >= 3; i -= 4, j++) {
222 data [j] = (uint)(
223 (inData [i-3] << (3*8)) |
224 (inData [i-2] << (2*8)) |
225 (inData [i-1] << (1*8)) |
226 (inData [i])
227 );
228 }
229
230 switch (leftOver) {
231 case 1: data [length-1] = (uint)inData [0]; break;
232 case 2: data [length-1] = (uint)((inData [0] << 8) | inData [1]); break;
233 case 3: data [length-1] = (uint)((inData [0] << 16) | (inData [1] << 8) | inData [2]); break;
234 }
235
236 this.Normalize ();
237 }
238
239 #if !INSIDE_CORLIB
240 [CLSCompliant (false)]
241 #endif
242 public BigInteger (uint [] inData)
243 {
244 if (inData.Length == 0)
245 inData = new uint [1];
246 length = (uint)inData.Length;
247
248 data = new uint [length];
249
250 for (int i = (int)length - 1, j = 0; i >= 0; i--, j++)
251 data [j] = inData [i];
252
253 this.Normalize ();
254 }
255
256 #if !INSIDE_CORLIB
257 [CLSCompliant (false)]
258 #endif
259 public BigInteger (uint ui)
260 {
261 data = new uint [] {ui};
262 }
263
264 #if !INSIDE_CORLIB
265 [CLSCompliant (false)]
266 #endif
267 public BigInteger (ulong ul)
268 {
269 data = new uint [2] { (uint)ul, (uint)(ul >> 32)};
270 length = 2;
271
272 this.Normalize ();
273 }
274
275 #if !INSIDE_CORLIB
276 [CLSCompliant (false)]
277 #endif
278 public static implicit operator BigInteger (uint value)
279 {
280 return (new BigInteger (value));
281 }
282
283 public static implicit operator BigInteger (int value)
284 {
285 if (value < 0) throw new ArgumentOutOfRangeException ("value");
286 return (new BigInteger ((uint)value));
287 }
288
289 #if !INSIDE_CORLIB
290 [CLSCompliant (false)]
291 #endif
292 public static implicit operator BigInteger (ulong value)
293 {
294 return (new BigInteger (value));
295 }
296
297 /* This is the BigInteger.Parse method I use. This method works
298 because BigInteger.ToString returns the input I gave to Parse. */
299 public static BigInteger Parse (string number)
300 {
301 if (number == null)
302 throw new ArgumentNullException ("number");
303
304 int i = 0, len = number.Length;
305 char c;
306 bool digits_seen = false;
307 BigInteger val = new BigInteger (0);
308 if (number [i] == '+') {
309 i++;
310 }
311 else if (number [i] == '-') {
312 throw new FormatException (WouldReturnNegVal);
313 }
314
315 for (; i < len; i++) {
316 c = number [i];
317 if (c == '\0') {
318 i = len;
319 continue;
320 }
321 if (c >= '0' && c <= '9') {
322 val = val * 10 + (c - '0');
323 digits_seen = true;
324 }
325 else {
326 if (Char.IsWhiteSpace (c)) {
327 for (i++; i < len; i++) {
328 if (!Char.IsWhiteSpace (number [i]))
329 throw new FormatException ();
330 }
331 break;
332 }
333 else
334 throw new FormatException ();
335 }
336 }
337 if (!digits_seen)
338 throw new FormatException ();
339 return val;
340 }
341
342 #endregion
343
344 #region Operators
345
346 public static BigInteger operator + (BigInteger bi1, BigInteger bi2)
347 {
348 if (bi1 == 0)
349 return new BigInteger (bi2);
350 else if (bi2 == 0)
351 return new BigInteger (bi1);
352 else
353 return Kernel.AddSameSign (bi1, bi2);
354 }
355
356 public static BigInteger operator - (BigInteger bi1, BigInteger bi2)
357 {
358 if (bi2 == 0)
359 return new BigInteger (bi1);
360
361 if (bi1 == 0)
362 throw new ArithmeticException (WouldReturnNegVal);
363
364 switch (Kernel.Compare (bi1, bi2)) {
365
366 case Sign.Zero:
367 return 0;
368
369 case Sign.Positive:
370 return Kernel.Subtract (bi1, bi2);
371
372 case Sign.Negative:
373 throw new ArithmeticException (WouldReturnNegVal);
374 default:
375 throw new Exception ();
376 }
377 }
378
379 public static int operator % (BigInteger bi, int i)
380 {
381 if (i > 0)
382 return (int)Kernel.DwordMod (bi, (uint)i);
383 else
384 return -(int)Kernel.DwordMod (bi, (uint)-i);
385 }
386
387 #if !INSIDE_CORLIB
388 [CLSCompliant (false)]
389 #endif
390 public static uint operator % (BigInteger bi, uint ui)
391 {
392 return Kernel.DwordMod (bi, (uint)ui);
393 }
394
395 public static BigInteger operator % (BigInteger bi1, BigInteger bi2)
396 {
397 return Kernel.multiByteDivide (bi1, bi2)[1];
398 }
399
400 public static BigInteger operator / (BigInteger bi, int i)
401 {
402 if (i > 0)
403 return Kernel.DwordDiv (bi, (uint)i);
404
405 throw new ArithmeticException (WouldReturnNegVal);
406 }
407
408 public static BigInteger operator / (BigInteger bi1, BigInteger bi2)
409 {
410 return Kernel.multiByteDivide (bi1, bi2)[0];
411 }
412
413 public static BigInteger operator * (BigInteger bi1, BigInteger bi2)
414 {
415 if (bi1 == 0 || bi2 == 0) return 0;
416
417 //
418 // Validate pointers
419 //
420 if (bi1.data.Length < bi1.length) throw new IndexOutOfRangeException ("bi1 out of range");
421 if (bi2.data.Length < bi2.length) throw new IndexOutOfRangeException ("bi2 out of range");
422
423 BigInteger ret = new BigInteger (Sign.Positive, bi1.length + bi2.length);
424
425 Kernel.Multiply (bi1.data, 0, bi1.length, bi2.data, 0, bi2.length, ret.data, 0);
426
427 ret.Normalize ();
428 return ret;
429 }
430
431 public static BigInteger operator * (BigInteger bi, int i)
432 {
433 if (i < 0) throw new ArithmeticException (WouldReturnNegVal);
434 if (i == 0) return 0;
435 if (i == 1) return new BigInteger (bi);
436
437 return Kernel.MultiplyByDword (bi, (uint)i);
438 }
439
440 public static BigInteger operator << (BigInteger bi1, int shiftVal)
441 {
442 return Kernel.LeftShift (bi1, shiftVal);
443 }
444
445 public static BigInteger operator >> (BigInteger bi1, int shiftVal)
446 {
447 return Kernel.RightShift (bi1, shiftVal);
448 }
449
450 #endregion
451
452 #region Friendly names for operators
453
454 // with names suggested by FxCop 1.30
455
456 public static BigInteger Add (BigInteger bi1, BigInteger bi2)
457 {
458 return (bi1 + bi2);
459 }
460
461 public static BigInteger Subtract (BigInteger bi1, BigInteger bi2)
462 {
463 return (bi1 - bi2);
464 }
465
466 public static int Modulus (BigInteger bi, int i)
467 {
468 return (bi % i);
469 }
470
471 #if !INSIDE_CORLIB
472 [CLSCompliant (false)]
473 #endif
474 public static uint Modulus (BigInteger bi, uint ui)
475 {
476 return (bi % ui);
477 }
478
479 public static BigInteger Modulus (BigInteger bi1, BigInteger bi2)
480 {
481 return (bi1 % bi2);
482 }
483
484 public static BigInteger Divid (BigInteger bi, int i)
485 {
486 return (bi / i);
487 }
488
489 public static BigInteger Divid (BigInteger bi1, BigInteger bi2)
490 {
491 return (bi1 / bi2);
492 }
493
494 public static BigInteger Multiply (BigInteger bi1, BigInteger bi2)
495 {
496 return (bi1 * bi2);
497 }
498
499 public static BigInteger Multiply (BigInteger bi, int i)
500 {
501 return (bi * i);
502 }
503
504 #endregion
505
506 #region Random
507 private static RandomNumberGenerator rng;
508 private static RandomNumberGenerator Rng {
509 get {
510 if (rng == null)
511 rng = RandomNumberGenerator.Create ();
512 return rng;
513 }
514 }
515
516 /// <summary>
517 /// Generates a new, random BigInteger of the specified length.
518 /// </summary>
519 /// <param name="bits">The number of bits for the new number.</param>
520 /// <param name="rng">A random number generator to use to obtain the bits.</param>
521 /// <returns>A random number of the specified length.</returns>
522 public static BigInteger GenerateRandom (int bits, RandomNumberGenerator rng)
523 {
524 int dwords = bits >> 5;
525 int remBits = bits & 0x1F;
526
527 if (remBits != 0)
528 dwords++;
529
530 BigInteger ret = new BigInteger (Sign.Positive, (uint)dwords + 1);
531 byte [] random = new byte [dwords << 2];
532
533 rng.GetBytes (random);
534 Buffer.BlockCopy (random, 0, ret.data, 0, (int)dwords << 2);
535
536 if (remBits != 0) {
537 uint mask = (uint)(0x01 << (remBits-1));
538 ret.data [dwords-1] |= mask;
539
540 mask = (uint)(0xFFFFFFFF >> (32 - remBits));
541 ret.data [dwords-1] &= mask;
542 }
543 else
544 ret.data [dwords-1] |= 0x80000000;
545
546 ret.Normalize ();
547 return ret;
548 }
549
550 /// <summary>
551 /// Generates a new, random BigInteger of the specified length using the default RNG crypto service provider.
552 /// </summary>
553 /// <param name="bits">The number of bits for the new number.</param>
554 /// <returns>A random number of the specified length.</returns>
555 public static BigInteger GenerateRandom (int bits)
556 {
557 return GenerateRandom (bits, Rng);
558 }
559
560 /// <summary>
561 /// Randomizes the bits in "this" from the specified RNG.
562 /// </summary>
563 /// <param name="rng">A RNG.</param>
564 public void Randomize (RandomNumberGenerator rng)
565 {
566 if (this == 0)
567 return;
568
569 int bits = this.BitCount ();
570 int dwords = bits >> 5;
571 int remBits = bits & 0x1F;
572
573 if (remBits != 0)
574 dwords++;
575
576 byte [] random = new byte [dwords << 2];
577
578 rng.GetBytes (random);
579 Buffer.BlockCopy (random, 0, data, 0, (int)dwords << 2);
580
581 if (remBits != 0) {
582 uint mask = (uint)(0x01 << (remBits-1));
583 data [dwords-1] |= mask;
584
585 mask = (uint)(0xFFFFFFFF >> (32 - remBits));
586 data [dwords-1] &= mask;
587 }
588
589 else
590 data [dwords-1] |= 0x80000000;
591
592 Normalize ();
593 }
594
595 /// <summary>
596 /// Randomizes the bits in "this" from the default RNG.
597 /// </summary>
598 public void Randomize ()
599 {
600 Randomize (Rng);
601 }
602
603 #endregion
604
605 #region Bitwise
606
607 public int BitCount ()
608 {
609 this.Normalize ();
610
611 uint value = data [length - 1];
612 uint mask = 0x80000000;
613 uint bits = 32;
614
615 while (bits > 0 && (value & mask) == 0) {
616 bits--;
617 mask >>= 1;
618 }
619 bits += ((length - 1) << 5);
620
621 return (int)bits;
622 }
623
624 /// <summary>
625 /// Tests if the specified bit is 1.
626 /// </summary>
627 /// <param name="bitNum">The bit to test. The least significant bit is 0.</param>
628 /// <returns>True if bitNum is set to 1, else false.</returns>
629 #if !INSIDE_CORLIB
630 [CLSCompliant (false)]
631 #endif
632 public bool TestBit (uint bitNum)
633 {
634 uint bytePos = bitNum >> 5; // divide by 32
635 byte bitPos = (byte)(bitNum & 0x1F); // get the lowest 5 bits
636
637 uint mask = (uint)1 << bitPos;
638 return ((this.data [bytePos] & mask) != 0);
639 }
640
641 public bool TestBit (int bitNum)
642 {
643 if (bitNum < 0) throw new IndexOutOfRangeException ("bitNum out of range");
644
645 uint bytePos = (uint)bitNum >> 5; // divide by 32
646 byte bitPos = (byte)(bitNum & 0x1F); // get the lowest 5 bits
647
648 uint mask = (uint)1 << bitPos;
649 return ((this.data [bytePos] | mask) == this.data [bytePos]);
650 }
651
652 #if !INSIDE_CORLIB
653 [CLSCompliant (false)]
654 #endif
655 public void SetBit (uint bitNum)
656 {
657 SetBit (bitNum, true);
658 }
659
660 #if !INSIDE_CORLIB
661 [CLSCompliant (false)]
662 #endif
663 public void ClearBit (uint bitNum)
664 {
665 SetBit (bitNum, false);
666 }
667
668 #if !INSIDE_CORLIB
669 [CLSCompliant (false)]
670 #endif
671 public void SetBit (uint bitNum, bool value)
672 {
673 uint bytePos = bitNum >> 5; // divide by 32
674
675 if (bytePos < this.length) {
676 uint mask = (uint)1 << (int)(bitNum & 0x1F);
677 if (value)
678 this.data [bytePos] |= mask;
679 else
680 this.data [bytePos] &= ~mask;
681 }
682 }
683
684 public int LowestSetBit ()
685 {
686 if (this == 0) return -1;
687 int i = 0;
688 while (!TestBit (i)) i++;
689 return i;
690 }
691
692 public byte[] GetBytes ()
693 {
694 if (this == 0) return new byte [1];
695
696 int numBits = BitCount ();
697 int numBytes = numBits >> 3;
698 if ((numBits & 0x7) != 0)
699 numBytes++;
700
701 byte [] result = new byte [numBytes];
702
703 int numBytesInWord = numBytes & 0x3;
704 if (numBytesInWord == 0) numBytesInWord = 4;
705
706 int pos = 0;
707 for (int i = (int)length - 1; i >= 0; i--) {
708 uint val = data [i];
709 for (int j = numBytesInWord - 1; j >= 0; j--) {
710 result [pos+j] = (byte)(val & 0xFF);
711 val >>= 8;
712 }
713 pos += numBytesInWord;
714 numBytesInWord = 4;
715 }
716 return result;
717 }
718
719 #endregion
720
721 #region Compare
722
723 #if !INSIDE_CORLIB
724 [CLSCompliant (false)]
725 #endif
726 public static bool operator == (BigInteger bi1, uint ui)
727 {
728 if (bi1.length != 1) bi1.Normalize ();
729 return bi1.length == 1 && bi1.data [0] == ui;
730 }
731
732 #if !INSIDE_CORLIB
733 [CLSCompliant (false)]
734 #endif
735 public static bool operator != (BigInteger bi1, uint ui)
736 {
737 if (bi1.length != 1) bi1.Normalize ();
738 return !(bi1.length == 1 && bi1.data [0] == ui);
739 }
740
741 public static bool operator == (BigInteger bi1, BigInteger bi2)
742 {
743 // we need to compare with null
744 if ((bi1 as object) == (bi2 as object))
745 return true;
746 if (null == bi1 || null == bi2)
747 return false;
748 return Kernel.Compare (bi1, bi2) == 0;
749 }
750
751 public static bool operator != (BigInteger bi1, BigInteger bi2)
752 {
753 // we need to compare with null
754 if ((bi1 as object) == (bi2 as object))
755 return false;
756 if (null == bi1 || null == bi2)
757 return true;
758 return Kernel.Compare (bi1, bi2) != 0;
759 }
760
761 public static bool operator > (BigInteger bi1, BigInteger bi2)
762 {
763 return Kernel.Compare (bi1, bi2) > 0;
764 }
765
766 public static bool operator < (BigInteger bi1, BigInteger bi2)
767 {
768 return Kernel.Compare (bi1, bi2) < 0;
769 }
770
771 public static bool operator >= (BigInteger bi1, BigInteger bi2)
772 {
773 return Kernel.Compare (bi1, bi2) >= 0;
774 }
775
776 public static bool operator <= (BigInteger bi1, BigInteger bi2)
777 {
778 return Kernel.Compare (bi1, bi2) <= 0;
779 }
780
781 public Sign Compare (BigInteger bi)
782 {
783 return Kernel.Compare (this, bi);
784 }
785
786 #endregion
787
788 #region Formatting
789
790 #if !INSIDE_CORLIB
791 [CLSCompliant (false)]
792 #endif
793 public string ToString (uint radix)
794 {
795 return ToString (radix, "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ");
796 }
797
798 #if !INSIDE_CORLIB
799 [CLSCompliant (false)]
800 #endif
801 public string ToString (uint radix, string characterSet)
802 {
803 if (characterSet.Length < radix)
804 throw new ArgumentException ("charSet length less than radix", "characterSet");
805 if (radix == 1)
806 throw new ArgumentException ("There is no such thing as radix one notation", "radix");
807
808 if (this == 0) return "0";
809 if (this == 1) return "1";
810
811 string result = "";
812
813 BigInteger a = new BigInteger (this);
814
815 while (a != 0) {
816 uint rem = Kernel.SingleByteDivideInPlace (a, radix);
817 result = characterSet [(int) rem] + result;
818 }
819
820 return result;
821 }
822
823 #endregion
824
825 #region Misc
826
827 /// <summary>
828 /// Normalizes this by setting the length to the actual number of
829 /// uints used in data and by setting the sign to Sign.Zero if the
830 /// value of this is 0.
831 /// </summary>
832 private void Normalize ()
833 {
834 // Normalize length
835 while (length > 0 && data [length-1] == 0) length--;
836
837 // Check for zero
838 if (length == 0)
839 length++;
840 }
841
842 public void Clear ()
843 {
844 for (int i=0; i < length; i++)
845 data [i] = 0x00;
846 }
847
848 #endregion
849
850 #region Object Impl
851
852 public override int GetHashCode ()
853 {
854 uint val = 0;
855
856 for (uint i = 0; i < this.length; i++)
857 val ^= this.data [i];
858
859 return (int)val;
860 }
861
862 public override string ToString ()
863 {
864 return ToString (10);
865 }
866
867 public override bool Equals (object o)
868 {
869 if (o == null)
870 return false;
871 if (o is int)
872 return (int)o >= 0 && this == (uint)o;
873
874 BigInteger bi = o as BigInteger;
875 if (bi == null)
876 return false;
877
878 return Kernel.Compare (this, bi) == 0;
879 }
880
881 #endregion
882
883 #region Number Theory
884
885 public BigInteger GCD (BigInteger bi)
886 {
887 return Kernel.gcd (this, bi);
888 }
889
890 public BigInteger ModInverse (BigInteger modulus)
891 {
892 return Kernel.modInverse (this, modulus);
893 }
894
895 public BigInteger ModPow (BigInteger exp, BigInteger n)
896 {
897 ModulusRing mr = new ModulusRing (n);
898 return mr.Pow (this, exp);
899 }
900
901 #endregion
902
903 #region Prime Testing
904
905 public bool IsProbablePrime ()
906 {
907 // can we use our small-prime table ?
908 if (this <= smallPrimes[smallPrimes.Length - 1]) {
909 for (int p = 0; p < smallPrimes.Length; p++) {
910 if (this == smallPrimes[p])
911 return true;
912 }
913 // the list is complete, so it's not a prime
914 return false;
915 }
916
917 // otherwise check if we can divide by one of the small primes
918 for (int p = 0; p < smallPrimes.Length; p++) {
919 if (this % smallPrimes[p] == 0)
920 return false;
921 }
922 // the last step is to confirm the "large" prime with the SPP or Miller-Rabin test
923 return PrimalityTests.Test (this, Prime.ConfidenceFactor.Medium);
924 }
925
926 #endregion
927
928 #region Prime Number Generation
929
930 /// <summary>
931 /// Generates the smallest prime >= bi
932 /// </summary>
933 /// <param name="bi">A BigInteger</param>
934 /// <returns>The smallest prime >= bi. More mathematically, if bi is prime: bi, else Prime [PrimePi [bi] + 1].</returns>
935 public static BigInteger NextHighestPrime (BigInteger bi)
936 {
937 NextPrimeFinder npf = new NextPrimeFinder ();
938 return npf.GenerateNewPrime (0, bi);
939 }
940
941 public static BigInteger GeneratePseudoPrime (int bits)
942 {
943 SequentialSearchPrimeGeneratorBase sspg = new SequentialSearchPrimeGeneratorBase ();
944 return sspg.GenerateNewPrime (bits);
945 }
946
947 /// <summary>
948 /// Increments this by two
949 /// </summary>
950 public void Incr2 ()
951 {
952 int i = 0;
953
954 data [0] += 2;
955
956 // If there was no carry, nothing to do
957 if (data [0] < 2) {
958
959 // Account for the first carry
960 data [++i]++;
961
962 // Keep adding until no carry
963 while (data [i++] == 0x0)
964 data [i]++;
965
966 // See if we increased the data length
967 if (length == (uint)i)
968 length++;
969 }
970 }
971
972 #endregion
973
974 #if INSIDE_CORLIB
975 internal
976 #else
977 public
978 #endif
979 sealed class ModulusRing {
980
981 BigInteger mod, constant;
982
983 public ModulusRing (BigInteger modulus)
984 {
985 this.mod = modulus;
986
987 // calculate constant = b^ (2k) / m
988 uint i = mod.length << 1;
989
990 constant = new BigInteger (Sign.Positive, i + 1);
991 constant.data [i] = 0x00000001;
992
993 constant = constant / mod;
994 }
995
996 public void BarrettReduction (BigInteger x)
997 {
998 BigInteger n = mod;
999 uint k = n.length,
1000 kPlusOne = k+1,
1001 kMinusOne = k-1;
1002
1003 // x < mod, so nothing to do.
1004 if (x.length < k) return;
1005
1006 BigInteger q3;
1007
1008 //
1009 // Validate pointers
1010 //
1011 if (x.data.Length < x.length) throw new IndexOutOfRangeException ("x out of range");
1012
1013 // q1 = x / b^ (k-1)
1014 // q2 = q1 * constant
1015 // q3 = q2 / b^ (k+1), Needs to be accessed with an offset of kPlusOne
1016
1017 // TODO: We should the method in HAC p 604 to do this (14.45)
1018 q3 = new BigInteger (Sign.Positive, x.length - kMinusOne + constant.length);
1019 Kernel.Multiply (x.data, kMinusOne, x.length - kMinusOne, constant.data, 0, constant.length, q3.data, 0);
1020
1021 // r1 = x mod b^ (k+1)
1022 // i.e. keep the lowest (k+1) words
1023
1024 uint lengthToCopy = (x.length > kPlusOne) ? kPlusOne : x.length;
1025
1026 x.length = lengthToCopy;
1027 x.Normalize ();
1028
1029 // r2 = (q3 * n) mod b^ (k+1)
1030 // partial multiplication of q3 and n
1031
1032 BigInteger r2 = new BigInteger (Sign.Positive, kPlusOne);
1033 Kernel.MultiplyMod2p32pmod (q3.data, (int)kPlusOne, (int)q3.length - (int)kPlusOne, n.data, 0, (int)n.length, r2.data, 0, (int)kPlusOne);
1034
1035 r2.Normalize ();
1036
1037 if (r2 <= x) {
1038 Kernel.MinusEq (x, r2);
1039 } else {
1040 BigInteger val = new BigInteger (Sign.Positive, kPlusOne + 1);
1041 val.data [kPlusOne] = 0x00000001;
1042
1043 Kernel.MinusEq (val, r2);
1044 Kernel.PlusEq (x, val);
1045 }
1046
1047 while (x >= n)
1048 Kernel.MinusEq (x, n);
1049 }
1050
1051 public BigInteger Multiply (BigInteger a, BigInteger b)
1052 {
1053 if (a == 0 || b == 0) return 0;
1054
1055 if (a > mod)
1056 a %= mod;
1057
1058 if (b > mod)
1059 b %= mod;
1060
1061 BigInteger ret = a * b;
1062 BarrettReduction (ret);
1063
1064 return ret;
1065 }
1066
1067 public BigInteger Difference (BigInteger a, BigInteger b)
1068 {
1069 Sign cmp = Kernel.Compare (a, b);
1070 BigInteger diff;
1071
1072 switch (cmp) {
1073 case Sign.Zero:
1074 return 0;
1075 case Sign.Positive:
1076 diff = a - b; break;
1077 case Sign.Negative:
1078 diff = b - a; break;
1079 default:
1080 throw new Exception ();
1081 }
1082
1083 if (diff >= mod) {
1084 if (diff.length >= mod.length << 1)
1085 diff %= mod;
1086 else
1087 BarrettReduction (diff);
1088 }
1089 if (cmp == Sign.Negative)
1090 diff = mod - diff;
1091 return diff;
1092 }
1093 #if true
1094 public BigInteger Pow (BigInteger a, BigInteger k)
1095 {
1096 BigInteger b = new BigInteger (1);
1097 if (k == 0)
1098 return b;
1099
1100 BigInteger A = a;
1101 if (k.TestBit (0))
1102 b = a;
1103
1104 for (int i = 1; i < k.BitCount (); i++) {
1105 A = Multiply (A, A);
1106 if (k.TestBit (i))
1107 b = Multiply (A, b);
1108 }
1109 return b;
1110 }
1111 #else
1112 public BigInteger Pow (BigInteger b, BigInteger exp)
1113 {
1114 if ((mod.data [0] & 1) == 1) return OddPow (b, exp);
1115 else return EvenPow (b, exp);
1116 }
1117
1118 public BigInteger EvenPow (BigInteger b, BigInteger exp)
1119 {
1120 BigInteger resultNum = new BigInteger ((BigInteger)1, mod.length << 1);
1121 BigInteger tempNum = new BigInteger (b % mod, mod.length << 1); // ensures (tempNum * tempNum) < b^ (2k)
1122
1123 uint totalBits = (uint)exp.BitCount ();
1124
1125 uint [] wkspace = new uint [mod.length << 1];
1126
1127 // perform squaring and multiply exponentiation
1128 for (uint pos = 0; pos < totalBits; pos++) {
1129 if (exp.TestBit (pos)) {
1130
1131 Array.Clear (wkspace, 0, wkspace.Length);
1132 Kernel.Multiply (resultNum.data, 0, resultNum.length, tempNum.data, 0, tempNum.length, wkspace, 0);
1133 resultNum.length += tempNum.length;
1134 uint [] t = wkspace;
1135 wkspace = resultNum.data;
1136 resultNum.data = t;
1137
1138 BarrettReduction (resultNum);
1139 }
1140
1141 Kernel.SquarePositive (tempNum, ref wkspace);
1142 BarrettReduction (tempNum);
1143
1144 if (tempNum == 1) {
1145 return resultNum;
1146 }
1147 }
1148
1149 return resultNum;
1150 }
1151
1152 private BigInteger OddPow (BigInteger b, BigInteger exp)
1153 {
1154 BigInteger resultNum = new BigInteger (Montgomery.ToMont (1, mod), mod.length << 1);
1155 BigInteger tempNum = new BigInteger (Montgomery.ToMont (b, mod), mod.length << 1); // ensures (tempNum * tempNum) < b^ (2k)
1156 uint mPrime = Montgomery.Inverse (mod.data [0]);
1157 uint totalBits = (uint)exp.BitCount ();
1158
1159 uint [] wkspace = new uint [mod.length << 1];
1160
1161 // perform squaring and multiply exponentiation
1162 for (uint pos = 0; pos < totalBits; pos++) {
1163 if (exp.TestBit (pos)) {
1164
1165 Array.Clear (wkspace, 0, wkspace.Length);
1166 Kernel.Multiply (resultNum.data, 0, resultNum.length, tempNum.data, 0, tempNum.length, wkspace, 0);
1167 resultNum.length += tempNum.length;
1168 uint [] t = wkspace;
1169 wkspace = resultNum.data;
1170 resultNum.data = t;
1171
1172 Montgomery.Reduce (resultNum, mod, mPrime);
1173 }
1174
1175 // the value of tempNum is required in the last loop
1176 if (pos < totalBits - 1) {
1177 Kernel.SquarePositive (tempNum, ref wkspace);
1178 Montgomery.Reduce (tempNum, mod, mPrime);
1179 }
1180 }
1181
1182 Montgomery.Reduce (resultNum, mod, mPrime);
1183 return resultNum;
1184 }
1185 #endif
1186 #region Pow Small Base
1187
1188 // TODO: Make tests for this, not really needed b/c prime stuff
1189 // checks it, but still would be nice
1190 #if !INSIDE_CORLIB
1191 [CLSCompliant (false)]
1192 #endif
1193 #if true
1194 public BigInteger Pow (uint b, BigInteger exp)
1195 {
1196 return Pow (new BigInteger (b), exp);
1197 }
1198 #else
1199 public BigInteger Pow (uint b, BigInteger exp)
1200 {
1201 // if (b != 2) {
1202 if ((mod.data [0] & 1) == 1)
1203 return OddPow (b, exp);
1204 else
1205 return EvenPow (b, exp);
1206 /* buggy in some cases (like the well tested primes)
1207 } else {
1208 if ((mod.data [0] & 1) == 1)
1209 return OddModTwoPow (exp);
1210 else
1211 return EvenModTwoPow (exp);
1212 }*/
1213 }
1214
1215 private unsafe BigInteger OddPow (uint b, BigInteger exp)
1216 {
1217 exp.Normalize ();
1218 uint [] wkspace = new uint [mod.length << 1 + 1];
1219
1220 BigInteger resultNum = Montgomery.ToMont ((BigInteger)b, this.mod);
1221 resultNum = new BigInteger (resultNum, mod.length << 1 +1);
1222
1223 uint mPrime = Montgomery.Inverse (mod.data [0]);
1224
1225 int bc = exp.BitCount () - 2;
1226 uint pos = (bc > 1 ? (uint) bc : 1);
1227
1228 //
1229 // We know that the first itr will make the val b
1230 //
1231
1232 do {
1233 //
1234 // r = r ^ 2 % m
1235 //
1236 Kernel.SquarePositive (resultNum, ref wkspace);
1237 resultNum = Montgomery.Reduce (resultNum, mod, mPrime);
1238
1239 if (exp.TestBit (pos)) {
1240
1241 //
1242 // r = r * b % m
1243 //
1244
1245 // TODO: Is Unsafe really speeding things up?
1246 fixed (uint* u = resultNum.data) {
1247
1248 uint i = 0;
1249 ulong mc = 0;
1250
1251 do {
1252 mc += (ulong)u [i] * (ulong)b;
1253 u [i] = (uint)mc;
1254 mc >>= 32;
1255 } while (++i < resultNum.length);
1256
1257 if (resultNum.length < mod.length) {
1258 if (mc != 0) {
1259 u [i] = (uint)mc;
1260 resultNum.length++;
1261 while (resultNum >= mod)
1262 Kernel.MinusEq (resultNum, mod);
1263 }
1264 } else if (mc != 0) {
1265
1266 //
1267 // First, we estimate the quotient by dividing
1268 // the first part of each of the numbers. Then
1269 // we correct this, if necessary, with a subtraction.
1270 //
1271
1272 uint cc = (uint)mc;
1273
1274 // We would rather have this estimate overshoot,
1275 // so we add one to the divisor
1276 uint divEstimate;
1277 if (mod.data [mod.length - 1] < UInt32.MaxValue) {
1278 divEstimate = (uint) ((((ulong)cc << 32) | (ulong) u [i -1]) /
1279 (mod.data [mod.length-1] + 1));
1280 }
1281 else {
1282 // guess but don't divide by 0
1283 divEstimate = (uint) ((((ulong)cc << 32) | (ulong) u [i -1]) /
1284 (mod.data [mod.length-1]));
1285 }
1286
1287 uint t;
1288
1289 i = 0;
1290 mc = 0;
1291 do {
1292 mc += (ulong)mod.data [i] * (ulong)divEstimate;
1293 t = u [i];
1294 u [i] -= (uint)mc;
1295 mc >>= 32;
1296 if (u [i] > t) mc++;
1297 i++;
1298 } while (i < resultNum.length);
1299 cc -= (uint)mc;
1300
1301 if (cc != 0) {
1302
1303 uint sc = 0, j = 0;
1304 uint [] s = mod.data;
1305 do {
1306 uint a = s [j];
1307 if (((a += sc) < sc) | ((u [j] -= a) > ~a)) sc = 1;
1308 else sc = 0;
1309 j++;
1310 } while (j < resultNum.length);
1311 cc -= sc;
1312 }
1313 while (resultNum >= mod)
1314 Kernel.MinusEq (resultNum, mod);
1315 } else {
1316 while (resultNum >= mod)
1317 Kernel.MinusEq (resultNum, mod);
1318 }
1319 }
1320 }
1321 } while (pos-- > 0);
1322
1323 resultNum = Montgomery.Reduce (resultNum, mod, mPrime);
1324 return resultNum;
1325
1326 }
1327
1328 private unsafe BigInteger EvenPow (uint b, BigInteger exp)
1329 {
1330 exp.Normalize ();
1331 uint [] wkspace = new uint [mod.length << 1 + 1];
1332 BigInteger resultNum = new BigInteger ((BigInteger)b, mod.length << 1 + 1);
1333
1334 uint pos = (uint)exp.BitCount () - 2;
1335
1336 //
1337 // We know that the first itr will make the val b
1338 //
1339
1340 do {
1341 //
1342 // r = r ^ 2 % m
1343 //
1344 Kernel.SquarePositive (resultNum, ref wkspace);
1345 if (!(resultNum.length < mod.length))
1346 BarrettReduction (resultNum);
1347
1348 if (exp.TestBit (pos)) {
1349
1350 //
1351 // r = r * b % m
1352 //
1353
1354 // TODO: Is Unsafe really speeding things up?
1355 fixed (uint* u = resultNum.data) {
1356
1357 uint i = 0;
1358 ulong mc = 0;
1359
1360 do {
1361 mc += (ulong)u [i] * (ulong)b;
1362 u [i] = (uint)mc;
1363 mc >>= 32;
1364 } while (++i < resultNum.length);
1365
1366 if (resultNum.length < mod.length) {
1367 if (mc != 0) {
1368 u [i] = (uint)mc;
1369 resultNum.length++;
1370 while (resultNum >= mod)
1371 Kernel.MinusEq (resultNum, mod);
1372 }
1373 } else if (mc != 0) {
1374
1375 //
1376 // First, we estimate the quotient by dividing
1377 // the first part of each of the numbers. Then
1378 // we correct this, if necessary, with a subtraction.
1379 //
1380
1381 uint cc = (uint)mc;
1382
1383 // We would rather have this estimate overshoot,
1384 // so we add one to the divisor
1385 uint divEstimate = (uint) ((((ulong)cc << 32) | (ulong) u [i -1]) /
1386 (mod.data [mod.length-1] + 1));
1387
1388 uint t;
1389
1390 i = 0;
1391 mc = 0;
1392 do {
1393 mc += (ulong)mod.data [i] * (ulong)divEstimate;
1394 t = u [i];
1395 u [i] -= (uint)mc;
1396 mc >>= 32;
1397 if (u [i] > t) mc++;
1398 i++;
1399 } while (i < resultNum.length);
1400 cc -= (uint)mc;
1401
1402 if (cc != 0) {
1403
1404 uint sc = 0, j = 0;
1405 uint [] s = mod.data;
1406 do {
1407 uint a = s [j];
1408 if (((a += sc) < sc) | ((u [j] -= a) > ~a)) sc = 1;
1409 else sc = 0;
1410 j++;
1411 } while (j < resultNum.length);
1412 cc -= sc;
1413 }
1414 while (resultNum >= mod)
1415 Kernel.MinusEq (resultNum, mod);
1416 } else {
1417 while (resultNum >= mod)
1418 Kernel.MinusEq (resultNum, mod);
1419 }
1420 }
1421 }
1422 } while (pos-- > 0);
1423
1424 return resultNum;
1425 }
1426 #endif
1427 /* known to be buggy in some cases */
1428 #if false
1429 private unsafe BigInteger EvenModTwoPow (BigInteger exp)
1430 {
1431 exp.Normalize ();
1432 uint [] wkspace = new uint [mod.length << 1 + 1];
1433
1434 BigInteger resultNum = new BigInteger (2, mod.length << 1 +1);
1435
1436 uint value = exp.data [exp.length - 1];
1437 uint mask = 0x80000000;
1438
1439 // Find the first bit of the exponent
1440 while ((value & mask) == 0)
1441 mask >>= 1;
1442
1443 //
1444 // We know that the first itr will make the val 2,
1445 // so eat one bit of the exponent
1446 //
1447 mask >>= 1;
1448
1449 uint wPos = exp.length - 1;
1450
1451 do {
1452 value = exp.data [wPos];
1453 do {
1454 Kernel.SquarePositive (resultNum, ref wkspace);
1455 if (resultNum.length >= mod.length)
1456 BarrettReduction (resultNum);
1457
1458 if ((value & mask) != 0) {
1459 //
1460 // resultNum = (resultNum * 2) % mod
1461 //
1462
1463 fixed (uint* u = resultNum.data) {
1464 //
1465 // Double
1466 //
1467 uint* uu = u;
1468 uint* uuE = u + resultNum.length;
1469 uint x, carry = 0;
1470 while (uu < uuE) {
1471 x = *uu;
1472 *uu = (x << 1) | carry;
1473 carry = x >> (32 - 1);
1474 uu++;
1475 }
1476
1477 // subtraction inlined because we know it is square
1478 if (carry != 0 || resultNum >= mod) {
1479 uu = u;
1480 uint c = 0;
1481 uint [] s = mod.data;
1482 uint i = 0;
1483 do {
1484 uint a = s [i];
1485 if (((a += c) < c) | ((* (uu++) -= a) > ~a))
1486 c = 1;
1487 else
1488 c = 0;
1489 i++;
1490 } while (uu < uuE);
1491 }
1492 }
1493 }
1494 } while ((mask >>= 1) > 0);
1495 mask = 0x80000000;
1496 } while (wPos-- > 0);
1497
1498 return resultNum;
1499 }
1500
1501 private unsafe BigInteger OddModTwoPow (BigInteger exp)
1502 {
1503
1504 uint [] wkspace = new uint [mod.length << 1 + 1];
1505
1506 BigInteger resultNum = Montgomery.ToMont ((BigInteger)2, this.mod);
1507 resultNum = new BigInteger (resultNum, mod.length << 1 +1);
1508
1509 uint mPrime = Montgomery.Inverse (mod.data [0]);
1510
1511 //
1512 // TODO: eat small bits, the ones we can do with no modular reduction
1513 //
1514 uint pos = (uint)exp.BitCount () - 2;
1515
1516 do {
1517 Kernel.SquarePositive (resultNum, ref wkspace);
1518 resultNum = Montgomery.Reduce (resultNum, mod, mPrime);
1519
1520 if (exp.TestBit (pos)) {
1521 //
1522 // resultNum = (resultNum * 2) % mod
1523 //
1524
1525 fixed (uint* u = resultNum.data) {
1526 //
1527 // Double
1528 //
1529 uint* uu = u;
1530 uint* uuE = u + resultNum.length;
1531 uint x, carry = 0;
1532 while (uu < uuE) {
1533 x = *uu;
1534 *uu = (x << 1) | carry;
1535 carry = x >> (32 - 1);
1536 uu++;
1537 }
1538
1539 // subtraction inlined because we know it is square
1540 if (carry != 0 || resultNum >= mod) {
1541 fixed (uint* s = mod.data) {
1542 uu = u;
1543 uint c = 0;
1544 uint* ss = s;
1545 do {
1546 uint a = *ss++;
1547 if (((a += c) < c) | ((* (uu++) -= a) > ~a))
1548 c = 1;
1549 else
1550 c = 0;
1551 } while (uu < uuE);
1552 }
1553 }
1554 }
1555 }
1556 } while (pos-- > 0);
1557
1558 resultNum = Montgomery.Reduce (resultNum, mod, mPrime);
1559 return resultNum;
1560 }
1561 #endif
1562 #endregion
1563 }
1564
1565 /// <summary>
1566 /// Low level functions for the BigInteger
1567 /// </summary>
1568 private sealed class Kernel {
1569
1570 #region Addition/Subtraction
1571
1572 /// <summary>
1573 /// Adds two numbers with the same sign.
1574 /// </summary>
1575 /// <param name="bi1">A BigInteger</param>
1576 /// <param name="bi2">A BigInteger</param>
1577 /// <returns>bi1 + bi2</returns>
1578 public static BigInteger AddSameSign (BigInteger bi1, BigInteger bi2)
1579 {
1580 uint [] x, y;
1581 uint yMax, xMax, i = 0;
1582
1583 // x should be bigger
1584 if (bi1.length < bi2.length) {
1585 x = bi2.data;
1586 xMax = bi2.length;
1587 y = bi1.data;
1588 yMax = bi1.length;
1589 } else {
1590 x = bi1.data;
1591 xMax = bi1.length;
1592 y = bi2.data;
1593 yMax = bi2.length;
1594 }
1595
1596 BigInteger result = new BigInteger (Sign.Positive, xMax + 1);
1597
1598 uint [] r = result.data;
1599
1600 ulong sum = 0;
1601
1602 // Add common parts of both numbers
1603 do {
1604 sum = ((ulong)x [i]) + ((ulong)y [i]) + sum;
1605 r [i] = (uint)sum;
1606 sum >>= 32;
1607 } while (++i < yMax);
1608
1609 // Copy remainder of longer number while carry propagation is required
1610 bool carry = (sum != 0);
1611
1612 if (carry) {
1613
1614 if (i < xMax) {
1615 do
1616 carry = ((r [i] = x [i] + 1) == 0);
1617 while (++i < xMax && carry);
1618 }
1619
1620 if (carry) {
1621 r [i] = 1;
1622 result.length = ++i;
1623 return result;
1624 }
1625 }
1626
1627 // Copy the rest
1628 if (i < xMax) {
1629 do
1630 r [i] = x [i];
1631 while (++i < xMax);
1632 }
1633
1634 result.Normalize ();
1635 return result;
1636 }
1637
1638 public static BigInteger Subtract (BigInteger big, BigInteger small)
1639 {
1640 BigInteger result = new BigInteger (Sign.Positive, big.length);
1641
1642 uint [] r = result.data, b = big.data, s = small.data;
1643 uint i = 0, c = 0;
1644
1645 do {
1646
1647 uint x = s [i];
1648 if (((x += c) < c) | ((r [i] = b [i] - x) > ~x))
1649 c = 1;
1650 else
1651 c = 0;
1652
1653 } while (++i < small.length);
1654
1655 if (i == big.length) goto fixup;
1656
1657 if (c == 1) {
1658 do
1659 r [i] = b [i] - 1;
1660 while (b [i++] == 0 && i < big.length);
1661
1662 if (i == big.length) goto fixup;
1663 }
1664
1665 do
1666 r [i] = b [i];
1667 while (++i < big.length);
1668
1669 fixup:
1670
1671 result.Normalize ();
1672 return result;
1673 }
1674
1675 public static void MinusEq (BigInteger big, BigInteger small)
1676 {
1677 uint [] b = big.data, s = small.data;
1678 uint i = 0, c = 0;
1679
1680 do {
1681 uint x = s [i];
1682 if (((x += c) < c) | ((b [i] -= x) > ~x))
1683 c = 1;
1684 else
1685 c = 0;
1686 } while (++i < small.length);
1687
1688 if (i == big.length) goto fixup;
1689
1690 if (c == 1) {
1691 do
1692 b [i]--;
1693 while (b [i++] == 0 && i < big.length);
1694 }
1695
1696 fixup:
1697
1698 // Normalize length
1699 while (big.length > 0 && big.data [big.length-1] == 0) big.length--;
1700
1701 // Check for zero
1702 if (big.length == 0)
1703 big.length++;
1704
1705 }
1706
1707 public static void PlusEq (BigInteger bi1, BigInteger bi2)
1708 {
1709 uint [] x, y;
1710 uint yMax, xMax, i = 0;
1711 bool flag = false;
1712
1713 // x should be bigger
1714 if (bi1.length < bi2.length){
1715 flag = true;
1716 x = bi2.data;
1717 xMax = bi2.length;
1718 y = bi1.data;
1719 yMax = bi1.length;
1720 } else {
1721 x = bi1.data;
1722 xMax = bi1.length;
1723 y = bi2.data;
1724 yMax = bi2.length;
1725 }
1726
1727 uint [] r = bi1.data;
1728
1729 ulong sum = 0;
1730
1731 // Add common parts of both numbers
1732 do {
1733 sum += ((ulong)x [i]) + ((ulong)y [i]);
1734 r [i] = (uint)sum;
1735 sum >>= 32;
1736 } while (++i < yMax);
1737
1738 // Copy remainder of longer number while carry propagation is required
1739 bool carry = (sum != 0);
1740
1741 if (carry){
1742
1743 if (i < xMax) {
1744 do
1745 carry = ((r [i] = x [i] + 1) == 0);
1746 while (++i < xMax && carry);
1747 }
1748
1749 if (carry) {
1750 r [i] = 1;
1751 bi1.length = ++i;
1752 return;
1753 }
1754 }
1755
1756 // Copy the rest
1757 if (flag && i < xMax - 1) {
1758 do
1759 r [i] = x [i];
1760 while (++i < xMax);
1761 }
1762
1763 bi1.length = xMax + 1;
1764 bi1.Normalize ();
1765 }
1766
1767 #endregion
1768
1769 #region Compare
1770
1771 /// <summary>
1772 /// Compares two BigInteger
1773 /// </summary>
1774 /// <param name="bi1">A BigInteger</param>
1775 /// <param name="bi2">A BigInteger</param>
1776 /// <returns>The sign of bi1 - bi2</returns>
1777 public static Sign Compare (BigInteger bi1, BigInteger bi2)
1778 {
1779 //
1780 // Step 1. Compare the lengths
1781 //
1782 uint l1 = bi1.length, l2 = bi2.length;
1783
1784 while (l1 > 0 && bi1.data [l1-1] == 0) l1--;
1785 while (l2 > 0 && bi2.data [l2-1] == 0) l2--;
1786
1787 if (l1 == 0 && l2 == 0) return Sign.Zero;
1788
1789 // bi1 len < bi2 len
1790 if (l1 < l2) return Sign.Negative;
1791 // bi1 len > bi2 len
1792 else if (l1 > l2) return Sign.Positive;
1793
1794 //
1795 // Step 2. Compare the bits
1796 //
1797
1798 uint pos = l1 - 1;
1799
1800 while (pos != 0 && bi1.data [pos] == bi2.data [pos]) pos--;
1801
1802 if (bi1.data [pos] < bi2.data [pos])
1803 return Sign.Negative;
1804 else if (bi1.data [pos] > bi2.data [pos])
1805 return Sign.Positive;
1806 else
1807 return Sign.Zero;
1808 }
1809
1810 #endregion
1811
1812 #region Division
1813
1814 #region Dword
1815
1816 /// <summary>
1817 /// Performs n / d and n % d in one operation.
1818 /// </summary>
1819 /// <param name="n">A BigInteger, upon exit this will hold n / d</param>
1820 /// <param name="d">The divisor</param>
1821 /// <returns>n % d</returns>
1822 public static uint SingleByteDivideInPlace (BigInteger n, uint d)
1823 {
1824 ulong r = 0;
1825 uint i = n.length;
1826
1827 while (i-- > 0) {
1828 r <<= 32;
1829 r |= n.data [i];
1830 n.data [i] = (uint)(r / d);
1831 r %= d;
1832 }
1833 n.Normalize ();
1834
1835 return (uint)r;
1836 }
1837
1838 public static uint DwordMod (BigInteger n, uint d)
1839 {
1840 ulong r = 0;
1841 uint i = n.length;
1842
1843 while (i-- > 0) {
1844 r <<= 32;
1845 r |= n.data [i];
1846 r %= d;
1847 }
1848
1849 return (uint)r;
1850 }
1851
1852 public static BigInteger DwordDiv (BigInteger n, uint d)
1853 {
1854 BigInteger ret = new BigInteger (Sign.Positive, n.length);
1855
1856 ulong r = 0;
1857 uint i = n.length;
1858
1859 while (i-- > 0) {
1860 r <<= 32;
1861 r |= n.data [i];
1862 ret.data [i] = (uint)(r / d);
1863 r %= d;
1864 }
1865 ret.Normalize ();
1866
1867 return ret;
1868 }
1869
1870 public static BigInteger [] DwordDivMod (BigInteger n, uint d)
1871 {
1872 BigInteger ret = new BigInteger (Sign.Positive , n.length);
1873
1874 ulong r = 0;
1875 uint i = n.length;
1876
1877 while (i-- > 0) {
1878 r <<= 32;
1879 r |= n.data [i];
1880 ret.data [i] = (uint)(r / d);
1881 r %= d;
1882 }
1883 ret.Normalize ();
1884
1885 BigInteger rem = (uint)r;
1886
1887 return new BigInteger [] {ret, rem};
1888 }
1889
1890 #endregion
1891
1892 #region BigNum
1893
1894 public static BigInteger [] multiByteDivide (BigInteger bi1, BigInteger bi2)
1895 {
1896 if (Kernel.Compare (bi1, bi2) == Sign.Negative)
1897 return new BigInteger [2] { 0, new BigInteger (bi1) };
1898
1899 bi1.Normalize (); bi2.Normalize ();
1900
1901 if (bi2.length == 1)
1902 return DwordDivMod (bi1, bi2.data [0]);
1903
1904 uint remainderLen = bi1.length + 1;
1905 int divisorLen = (int)bi2.length + 1;
1906
1907 uint mask = 0x80000000;
1908 uint val = bi2.data [bi2.length - 1];
1909 int shift = 0;
1910 int resultPos = (int)bi1.length - (int)bi2.length;
1911
1912 while (mask != 0 && (val & mask) == 0) {
1913 shift++; mask >>= 1;
1914 }
1915
1916 BigInteger quot = new BigInteger (Sign.Positive, bi1.length - bi2.length + 1);
1917 BigInteger rem = (bi1 << shift);
1918
1919 uint [] remainder = rem.data;
1920
1921 bi2 = bi2 << shift;
1922
1923 int j = (int)(remainderLen - bi2.length);
1924 int pos = (int)remainderLen - 1;
1925
1926 uint firstDivisorByte = bi2.data [bi2.length-1];
1927 ulong secondDivisorByte = bi2.data [bi2.length-2];
1928
1929 while (j > 0) {
1930 ulong dividend = ((ulong)remainder [pos] << 32) + (ulong)remainder [pos-1];
1931
1932 ulong q_hat = dividend / (ulong)firstDivisorByte;
1933 ulong r_hat = dividend % (ulong)firstDivisorByte;
1934
1935 do {
1936
1937 if (q_hat == 0x100000000 ||
1938 (q_hat * secondDivisorByte) > ((r_hat << 32) + remainder [pos-2])) {
1939 q_hat--;
1940 r_hat += (ulong)firstDivisorByte;
1941
1942 if (r_hat < 0x100000000)
1943 continue;
1944 }
1945 break;
1946 } while (true);
1947
1948 //
1949 // At this point, q_hat is either exact, or one too large
1950 // (more likely to be exact) so, we attempt to multiply the
1951 // divisor by q_hat, if we get a borrow, we just subtract
1952 // one from q_hat and add the divisor back.
1953 //
1954
1955 uint t;
1956 uint dPos = 0;
1957 int nPos = pos - divisorLen + 1;
1958 ulong mc = 0;
1959 uint uint_q_hat = (uint)q_hat;
1960 do {
1961 mc += (ulong)bi2.data [dPos] * (ulong)uint_q_hat;
1962 t = remainder [nPos];
1963 remainder [nPos] -= (uint)mc;
1964 mc >>= 32;
1965 if (remainder [nPos] > t) mc++;
1966 dPos++; nPos++;
1967 } while (dPos < divisorLen);
1968
1969 nPos = pos - divisorLen + 1;
1970 dPos = 0;
1971
1972 // Overestimate
1973 if (mc != 0) {
1974 uint_q_hat--;
1975 ulong sum = 0;
1976
1977 do {
1978 sum = ((ulong)remainder [nPos]) + ((ulong)bi2.data [dPos]) + sum;
1979 remainder [nPos] = (uint)sum;
1980 sum >>= 32;
1981 dPos++; nPos++;
1982 } while (dPos < divisorLen);
1983
1984 }
1985
1986 quot.data [resultPos--] = (uint)uint_q_hat;
1987
1988 pos--;
1989 j--;
1990 }
1991
1992 quot.Normalize ();
1993 rem.Normalize ();
1994 BigInteger [] ret = new BigInteger [2] { quot, rem };
1995
1996 if (shift != 0)
1997 ret [1] >>= shift;
1998
1999 return ret;
2000 }
2001
2002 #endregion
2003
2004 #endregion
2005
2006 #region Shift
2007 public static BigInteger LeftShift (BigInteger bi, int n)
2008 {
2009 if (n == 0) return new BigInteger (bi, bi.length + 1);
2010
2011 int w = n >> 5;
2012 n &= ((1 << 5) - 1);
2013
2014 BigInteger ret = new BigInteger (Sign.Positive, bi.length + 1 + (uint)w);
2015
2016 uint i = 0, l = bi.length;
2017 if (n != 0) {
2018 uint x, carry = 0;
2019 while (i < l) {
2020 x = bi.data [i];
2021 ret.data [i + w] = (x << n) | carry;
2022 carry = x >> (32 - n);
2023 i++;
2024 }
2025 ret.data [i + w] = carry;
2026 } else {
2027 while (i < l) {
2028 ret.data [i + w] = bi.data [i];
2029 i++;
2030 }
2031 }
2032
2033 ret.Normalize ();
2034 return ret;
2035 }
2036
2037 public static BigInteger RightShift (BigInteger bi, int n)
2038 {
2039 if (n == 0) return new BigInteger (bi);
2040
2041 int w = n >> 5;
2042 int s = n & ((1 << 5) - 1);
2043
2044 BigInteger ret = new BigInteger (Sign.Positive, bi.length - (uint)w + 1);
2045 uint l = (uint)ret.data.Length - 1;
2046
2047 if (s != 0) {
2048
2049 uint x, carry = 0;
2050
2051 while (l-- > 0) {
2052 x = bi.data [l + w];
2053 ret.data [l] = (x >> n) | carry;
2054 carry = x << (32 - n);
2055 }
2056 } else {
2057 while (l-- > 0)
2058 ret.data [l] = bi.data [l + w];
2059
2060 }
2061 ret.Normalize ();
2062 return ret;
2063 }
2064
2065 #endregion
2066
2067 #region Multiply
2068
2069 public static BigInteger MultiplyByDword (BigInteger n, uint f)
2070 {
2071 BigInteger ret = new BigInteger (Sign.Positive, n.length + 1);
2072
2073 uint i = 0;
2074 ulong c = 0;
2075
2076 do {
2077 c += (ulong)n.data [i] * (ulong)f;
2078 ret.data [i] = (uint)c;
2079 c >>= 32;
2080 } while (++i < n.length);
2081 ret.data [i] = (uint)c;
2082 ret.Normalize ();
2083 return ret;
2084
2085 }
2086
2087 /// <summary>
2088 /// Multiplies the data in x [xOffset:xOffset+xLen] by
2089 /// y [yOffset:yOffset+yLen] and puts it into
2090 /// d [dOffset:dOffset+xLen+yLen].
2091 /// </summary>
2092 /// <remarks>
2093 /// This code is unsafe! It is the caller's responsibility to make
2094 /// sure that it is safe to access x [xOffset:xOffset+xLen],
2095 /// y [yOffset:yOffset+yLen], and d [dOffset:dOffset+xLen+yLen].
2096 /// </remarks>
2097 public static unsafe void Multiply (uint [] x, uint xOffset, uint xLen, uint [] y, uint yOffset, uint yLen, uint [] d, uint dOffset)
2098 {
2099 fixed (uint* xx = x, yy = y, dd = d) {
2100 uint* xP = xx + xOffset,
2101 xE = xP + xLen,
2102 yB = yy + yOffset,
2103 yE = yB + yLen,
2104 dB = dd + dOffset;
2105
2106 for (; xP < xE; xP++, dB++) {
2107
2108 if (*xP == 0) continue;
2109
2110 ulong mcarry = 0;
2111
2112 uint* dP = dB;
2113 for (uint* yP = yB; yP < yE; yP++, dP++) {
2114 mcarry += ((ulong)*xP * (ulong)*yP) + (ulong)*dP;
2115
2116 *dP = (uint)mcarry;
2117 mcarry >>= 32;
2118 }
2119
2120 if (mcarry != 0)
2121 *dP = (uint)mcarry;
2122 }
2123 }
2124 }
2125
2126 /// <summary>
2127 /// Multiplies the data in x [xOffset:xOffset+xLen] by
2128 /// y [yOffset:yOffset+yLen] and puts the low mod words into
2129 /// d [dOffset:dOffset+mod].
2130 /// </summary>
2131 /// <remarks>
2132 /// This code is unsafe! It is the caller's responsibility to make
2133 /// sure that it is safe to access x [xOffset:xOffset+xLen],
2134 /// y [yOffset:yOffset+yLen], and d [dOffset:dOffset+mod].
2135 /// </remarks>
2136 public static unsafe void MultiplyMod2p32pmod (uint [] x, int xOffset, int xLen, uint [] y, int yOffest, int yLen, uint [] d, int dOffset, int mod)
2137 {
2138 fixed (uint* xx = x, yy = y, dd = d) {
2139 uint* xP = xx + xOffset,
2140 xE = xP + xLen,
2141 yB = yy + yOffest,
2142 yE = yB + yLen,
2143 dB = dd + dOffset,
2144 dE = dB + mod;
2145
2146 for (; xP < xE; xP++, dB++) {
2147
2148 if (*xP == 0) continue;
2149
2150 ulong mcarry = 0;
2151 uint* dP = dB;
2152 for (uint* yP = yB; yP < yE && dP < dE; yP++, dP++) {
2153 mcarry += ((ulong)*xP * (ulong)*yP) + (ulong)*dP;
2154
2155 *dP = (uint)mcarry;
2156 mcarry >>= 32;
2157 }
2158
2159 if (mcarry != 0 && dP < dE)
2160 *dP = (uint)mcarry;
2161 }
2162 }
2163 }
2164
2165 public static unsafe void SquarePositive (BigInteger bi, ref uint [] wkSpace)
2166 {
2167 uint [] t = wkSpace;
2168 wkSpace = bi.data;
2169 uint [] d = bi.data;
2170 uint dl = bi.length;
2171 bi.data = t;
2172
2173 fixed (uint* dd = d, tt = t) {
2174
2175 uint* ttE = tt + t.Length;
2176 // Clear the dest
2177 for (uint* ttt = tt; ttt < ttE; ttt++)
2178 *ttt = 0;
2179
2180 uint* dP = dd, tP = tt;
2181
2182 for (uint i = 0; i < dl; i++, dP++) {
2183 if (*dP == 0)
2184 continue;
2185
2186 ulong mcarry = 0;
2187 uint bi1val = *dP;
2188
2189 uint* dP2 = dP + 1, tP2 = tP + 2*i + 1;
2190
2191 for (uint j = i + 1; j < dl; j++, tP2++, dP2++) {
2192 // k = i + j
2193 mcarry += ((ulong)bi1val * (ulong)*dP2) + *tP2;
2194
2195 *tP2 = (uint)mcarry;
2196 mcarry >>= 32;
2197 }
2198
2199 if (mcarry != 0)
2200 *tP2 = (uint)mcarry;
2201 }
2202
2203 // Double t. Inlined for speed.
2204
2205 tP = tt;
2206
2207 uint x, carry = 0;
2208 while (tP < ttE) {
2209 x = *tP;
2210 *tP = (x << 1) | carry;
2211 carry = x >> (32 - 1);
2212 tP++;
2213 }
2214 if (carry != 0) *tP = carry;
2215
2216 // Add in the diagnals
2217
2218 dP = dd;
2219 tP = tt;
2220 for (uint* dE = dP + dl; (dP < dE); dP++, tP++) {
2221 ulong val = (ulong)*dP * (ulong)*dP + *tP;
2222 *tP = (uint)val;
2223 val >>= 32;
2224 *(++tP) += (uint)val;
2225 if (*tP < (uint)val) {
2226 uint* tP3 = tP;
2227 // Account for the first carry
2228 (*++tP3)++;
2229
2230 // Keep adding until no carry
2231 while ((*tP3++) == 0)
2232 (*tP3)++;
2233 }
2234
2235 }
2236
2237 bi.length <<= 1;
2238
2239 // Normalize length
2240 while (tt [bi.length-1] == 0 && bi.length > 1) bi.length--;
2241
2242 }
2243 }
2244
2245 /*
2246 * Never called in BigInteger (and part of a private class)
2247 * public static bool Double (uint [] u, int l)
2248 {
2249 uint x, carry = 0;
2250 uint i = 0;
2251 while (i < l) {
2252 x = u [i];
2253 u [i] = (x << 1) | carry;
2254 carry = x >> (32 - 1);
2255 i++;
2256 }
2257 if (carry != 0) u [l] = carry;
2258 return carry != 0;
2259 }*/
2260
2261 #endregion
2262
2263 #region Number Theory
2264
2265 public static BigInteger gcd (BigInteger a, BigInteger b)
2266 {
2267 BigInteger x = a;
2268 BigInteger y = b;
2269
2270 BigInteger g = y;
2271
2272 while (x.length > 1) {
2273 g = x;
2274 x = y % x;
2275 y = g;
2276
2277 }
2278 if (x == 0) return g;
2279
2280 // TODO: should we have something here if we can convert to long?
2281
2282 //
2283 // Now we can just do it with single precision. I am using the binary gcd method,
2284 // as it should be faster.
2285 //
2286
2287 uint yy = x.data [0];
2288 uint xx = y % yy;
2289
2290 int t = 0;
2291
2292 while (((xx | yy) & 1) == 0) {
2293 xx >>= 1; yy >>= 1; t++;
2294 }
2295 while (xx != 0) {
2296 while ((xx & 1) == 0) xx >>= 1;
2297 while ((yy & 1) == 0) yy >>= 1;
2298 if (xx >= yy)
2299 xx = (xx - yy) >> 1;
2300 else
2301 yy = (yy - xx) >> 1;
2302 }
2303
2304 return yy << t;
2305 }
2306
2307 public static uint modInverse (BigInteger bi, uint modulus)
2308 {
2309 uint a = modulus, b = bi % modulus;
2310 uint p0 = 0, p1 = 1;
2311
2312 while (b != 0) {
2313 if (b == 1)
2314 return p1;
2315 p0 += (a / b) * p1;
2316 a %= b;
2317
2318 if (a == 0)
2319 break;
2320 if (a == 1)
2321 return modulus-p0;
2322
2323 p1 += (b / a) * p0;
2324 b %= a;
2325
2326 }
2327 return 0;
2328 }
2329
2330 public static BigInteger modInverse (BigInteger bi, BigInteger modulus)
2331 {
2332 if (modulus.length == 1) return modInverse (bi, modulus.data [0]);
2333
2334 BigInteger [] p = { 0, 1 };
2335 BigInteger [] q = new BigInteger [2]; // quotients
2336 BigInteger [] r = { 0, 0 }; // remainders
2337
2338 int step = 0;
2339
2340 BigInteger a = modulus;
2341 BigInteger b = bi;
2342
2343 ModulusRing mr = new ModulusRing (modulus);
2344
2345 while (b != 0) {
2346
2347 if (step > 1) {
2348
2349 BigInteger pval = mr.Difference (p [0], p [1] * q [0]);
2350 p [0] = p [1]; p [1] = pval;
2351 }
2352
2353 BigInteger [] divret = multiByteDivide (a, b);
2354
2355 q [0] = q [1]; q [1] = divret [0];
2356 r [0] = r [1]; r [1] = divret [1];
2357 a = b;
2358 b = divret [1];
2359
2360 step++;
2361 }
2362
2363 if (r [0] != 1)
2364 throw (new ArithmeticException ("No inverse!"));
2365
2366 return mr.Difference (p [0], p [1] * q [0]);
2367
2368 }
2369 #endregion
2370 }
2371 }
2372 }
Mirror of the official Mono project git repository