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