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