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