source: trunk/CrypPluginBase/Miscellaneous/BigInteger.cs @ 1155

Last change on this file since 1155 was 1155, checked in by Arno Wacker, 12 years ago

P2P Bruteforce key searching

  • Fixed serialization bug (JobId > 255)
  • Removed special case for JobId == 0 due to BigInt-bug; BigInt bug fixed, see bellow
  • Cleanup and renaming of the P2P-samples for Worker and Manager

BigInt:

  • getBytes() now returns also a single zeroed byte if value of BigInt is zero
  • added a warning to BigInt.dataLength - in some cases the values returned are wrong ==> yet another bug

AnotherEditor:

  • Invoke -> BeginInvoke for Plugin-Icon-changes and Quickwatch-updates
  • New binaries (CrypWin/AnotherEditor)
File size: 127.4 KB
Line 
1//source: http://www.codeproject.com/KB/cs/biginteger.aspx
2//************************************************************************************
3// BigInteger Class Version 1.03
4//
5// Copyright (c) 2002 Chew Keong TAN
6// All rights reserved.
7//
8// Permission is hereby granted, free of charge, to any person obtaining a
9// copy of this software and associated documentation files (the
10// "Software"), to deal in the Software without restriction, including
11// without limitation the rights to use, copy, modify, merge, publish,
12// distribute, and/or sell copies of the Software, and to permit persons
13// to whom the Software is furnished to do so, provided that the above
14// copyright notice(s) and this permission notice appear in all copies of
15// the Software and that both the above copyright notice(s) and this
16// permission notice appear in supporting documentation.
17//
18// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
19// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
20// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT
21// OF THIRD PARTY RIGHTS. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
22// HOLDERS INCLUDED IN THIS NOTICE BE LIABLE FOR ANY CLAIM, OR ANY SPECIAL
23// INDIRECT OR CONSEQUENTIAL DAMAGES, OR ANY DAMAGES WHATSOEVER RESULTING
24// FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
25// NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION
26// WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
27//
28//
29// Disclaimer
30// ----------
31// Although reasonable care has been taken to ensure the correctness of this
32// implementation, this code should never be used in any application without
33// proper verification and testing.  I disclaim all liability and responsibility
34// to any person or entity with respect to any loss or damage caused, or alleged
35// to be caused, directly or indirectly, by the use of this BigInteger class.
36//
37// Comments, bugs and suggestions to
38// (http://www.codeproject.com/csharp/biginteger.asp)
39//
40//
41// Overloaded Operators +, -, *, /, %, >>, <<, ==, !=, >, <, >=, <=, &, |, ^, ++, --, ~
42//
43// Features
44// --------
45// 1) Arithmetic operations involving large signed integers (2's complement).
46// 2) Primality test using Fermat little theorm, Rabin Miller's method,
47//    Solovay Strassen's method and Lucas strong pseudoprime.
48// 3) Modulo exponential with Barrett's reduction.
49// 4) Inverse modulo.
50// 5) Pseudo prime generation.
51// 6) Co-prime generation.
52//
53//
54// Known Problem
55// -------------
56// This pseudoprime passes my implementation of
57// primality test but failed in JDK's isProbablePrime test.
58//
59//       byte[] pseudoPrime1 = { (byte)0x00,
60//             (byte)0x85, (byte)0x84, (byte)0x64, (byte)0xFD, (byte)0x70, (byte)0x6A,
61//             (byte)0x9F, (byte)0xF0, (byte)0x94, (byte)0x0C, (byte)0x3E, (byte)0x2C,
62//             (byte)0x74, (byte)0x34, (byte)0x05, (byte)0xC9, (byte)0x55, (byte)0xB3,
63//             (byte)0x85, (byte)0x32, (byte)0x98, (byte)0x71, (byte)0xF9, (byte)0x41,
64//             (byte)0x21, (byte)0x5F, (byte)0x02, (byte)0x9E, (byte)0xEA, (byte)0x56,
65//             (byte)0x8D, (byte)0x8C, (byte)0x44, (byte)0xCC, (byte)0xEE, (byte)0xEE,
66//             (byte)0x3D, (byte)0x2C, (byte)0x9D, (byte)0x2C, (byte)0x12, (byte)0x41,
67//             (byte)0x1E, (byte)0xF1, (byte)0xC5, (byte)0x32, (byte)0xC3, (byte)0xAA,
68//             (byte)0x31, (byte)0x4A, (byte)0x52, (byte)0xD8, (byte)0xE8, (byte)0xAF,
69//             (byte)0x42, (byte)0xF4, (byte)0x72, (byte)0xA1, (byte)0x2A, (byte)0x0D,
70//             (byte)0x97, (byte)0xB1, (byte)0x31, (byte)0xB3,
71//       };
72//
73//
74// Change Log
75// ----------
76// 1) September 23, 2002 (Version 1.03)
77//    - Fixed operator- to give correct data length.
78//    - Added Lucas sequence generation.
79//    - Added Strong Lucas Primality test.
80//    - Added integer square root method.
81//    - Added setBit/unsetBit methods.
82//    - New isProbablePrime() method which do not require the
83//      confident parameter.
84//
85// 2) August 29, 2002 (Version 1.02)
86//    - Fixed bug in the exponentiation of negative numbers.
87//    - Faster modular exponentiation using Barrett reduction.
88//    - Added getBytes() method.
89//    - Fixed bug in ToHexString method.
90//    - Added overloading of ^ operator.
91//    - Faster computation of Jacobi symbol.
92//
93// 3) August 19, 2002 (Version 1.01)
94//    - Big integer is stored and manipulated as unsigned integers (4 bytes) instead of
95//      individual bytes this gives significant performance improvement.
96//    - Updated Fermat's Little Theorem test to use a^(p-1) mod p = 1
97//    - Added isProbablePrime method.
98//    - Updated documentation.
99//
100// 4) August 9, 2002 (Version 1.0)
101//    - Initial Release.
102//
103//
104// References
105// [1] D. E. Knuth, "Seminumerical Algorithms", The Art of Computer Programming Vol. 2,
106//     3rd Edition, Addison-Wesley, 1998.
107//
108// [2] K. H. Rosen, "Elementary Number Theory and Its Applications", 3rd Ed,
109//     Addison-Wesley, 1993.
110//
111// [3] B. Schneier, "Applied Cryptography", 2nd Ed, John Wiley & Sons, 1996.
112//
113// [4] A. Menezes, P. van Oorschot, and S. Vanstone, "Handbook of Applied Cryptography",
114//     CRC Press, 1996, www.cacr.math.uwaterloo.ca/hac
115//
116// [5] A. Bosselaers, R. Govaerts, and J. Vandewalle, "Comparison of Three Modular
117//     Reduction Functions," Proc. CRYPTO'93, pp.175-186.
118//
119// [6] R. Baillie and S. S. Wagstaff Jr, "Lucas Pseudoprimes", Mathematics of Computation,
120//     Vol. 35, No. 152, Oct 1980, pp. 1391-1417.
121//
122// [7] H. C. Williams, "Édouard Lucas and Primality Testing", Canadian Mathematical
123//     Society Series of Monographs and Advance Texts, vol. 22, John Wiley & Sons, New York,
124//     NY, 1998.
125//
126// [8] P. Ribenboim, "The new book of prime number records", 3rd edition, Springer-Verlag,
127//     New York, NY, 1995.
128//
129// [9] M. Joye and J.-J. Quisquater, "Efficient computation of full Lucas sequences",
130//     Electronics Letters, 32(6), 1996, pp 537-538.
131//
132//************************************************************************************
133
134using System;
135using System.Collections;
136using System.Collections.Generic;
137
138namespace Cryptool.PluginBase.Miscellaneous
139{
140    public class BigInteger : IComparable
141    {
142        // maximum length of the BigInteger in uint (4 bytes)
143        // change this to suit the required level of precision.
144
145        private const int maxLength = 1024;
146
147        // primes smaller than 2000 to test the generated prime number
148
149        public static readonly int[] primesBelow2000 = {
150            2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97,
151            101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199,
152            211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293,
153            307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397,
154            401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499,
155            503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599,
156            601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691,
157            701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797,
158            809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887,
159            907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997,
160            1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097,
161            1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193,
162            1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, 1297,
163            1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399,
164            1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499,
165            1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597,
166            1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699,
167            1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789,
168            1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889,
169            1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999 };
170
171
172        private uint[] data = null;             // stores bytes from the Big Integer
173       
174        /// <summary>
175        /// This returns the number of the actual chars used
176        /// WARNIG: this sees bugg, e.g. a vlue of 256 lead to dataLength = 1, which should be 2
177        /// This happens if 256 is created with the ++ operator/ other cases not checked)
178        /// </summary>
179        public int dataLength;
180
181
182        //***********************************************************************
183        // Constructor (Default value for BigInteger is 0
184        //***********************************************************************
185
186        public BigInteger()
187        {
188            data = new uint[maxLength];
189            dataLength = 1;
190        }
191
192
193        //***********************************************************************
194        // Constructor (Default value provided by long)
195        //***********************************************************************
196
197        public BigInteger(long value)
198        {
199            data = new uint[maxLength];
200            long tempVal = value;
201
202            // copy bytes from long to BigInteger without any assumption of
203            // the length of the long datatype
204
205            dataLength = 0;
206            while (value != 0 && dataLength < maxLength)
207            {
208                data[dataLength] = (uint)(value & 0xFFFFFFFF);
209                value >>= 32;
210                dataLength++;
211            }
212
213            if (tempVal > 0)         // overflow check for +ve value
214            {
215                if (value != 0 || (data[maxLength - 1] & 0x80000000) != 0)
216                    throw (new ArithmeticException("Positive overflow in constructor."));
217            }
218            else if (tempVal < 0)    // underflow check for -ve value
219            {
220                if (value != -1 || (data[dataLength - 1] & 0x80000000) == 0)
221                    throw (new ArithmeticException("Negative underflow in constructor."));
222            }
223
224            if (dataLength == 0)
225                dataLength = 1;
226        }
227
228
229        //***********************************************************************
230        // Constructor (Default value provided by ulong)
231        //***********************************************************************
232
233        public BigInteger(ulong value)
234        {
235            data = new uint[maxLength];
236
237            // copy bytes from ulong to BigInteger without any assumption of
238            // the length of the ulong datatype
239
240            dataLength = 0;
241            while (value != 0 && dataLength < maxLength)
242            {
243                data[dataLength] = (uint)(value & 0xFFFFFFFF);
244                value >>= 32;
245                dataLength++;
246            }
247
248            if (value != 0 || (data[maxLength - 1] & 0x80000000) != 0)
249                throw (new ArithmeticException("Positive overflow in constructor."));
250
251            if (dataLength == 0)
252                dataLength = 1;
253        }
254
255
256
257        //***********************************************************************
258        // Constructor (Default value provided by BigInteger)
259        //***********************************************************************
260
261        public BigInteger(BigInteger bi)
262        {
263            data = new uint[maxLength];
264
265            dataLength = bi.dataLength;
266
267            for (int i = 0; i < dataLength; i++)
268                data[i] = bi.data[i];
269        }
270
271
272        //***********************************************************************
273        // Constructor (Default value provided by a string of digits of the
274        //              specified base)
275        //
276        // Example (base 10)
277        // -----------------
278        // To initialize "a" with the default value of 1234 in base 10
279        //      BigInteger a = new BigInteger("1234", 10)
280        //
281        // To initialize "a" with the default value of -1234
282        //      BigInteger a = new BigInteger("-1234", 10)
283        //
284        // Example (base 16)
285        // -----------------
286        // To initialize "a" with the default value of 0x1D4F in base 16
287        //      BigInteger a = new BigInteger("1D4F", 16)
288        //
289        // To initialize "a" with the default value of -0x1D4F
290        //      BigInteger a = new BigInteger("-1D4F", 16)
291        //
292        // Note that string values are specified in the <sign><magnitude>
293        // format.
294        //
295        //***********************************************************************
296
297        public BigInteger(string value, int radix)
298        {
299            BigInteger multiplier = new BigInteger(1);
300            BigInteger result = new BigInteger();
301            value = (value.ToUpper()).Trim();
302            int limit = 0;
303
304            if (value[0] == '-')
305                limit = 1;
306
307            for (int i = value.Length - 1; i >= limit; i--)
308            {
309                int posVal = (int)value[i];
310
311                if (posVal >= '0' && posVal <= '9')
312                    posVal -= '0';
313                else if (posVal >= 'A' && posVal <= 'Z')
314                    posVal = (posVal - 'A') + 10;
315                else
316                    posVal = 9999999;       // arbitrary large
317
318
319                if (posVal >= radix)
320                    throw (new ArithmeticException("Invalid string in constructor."));
321                else
322                {
323                    if (value[0] == '-')
324                        posVal = -posVal;
325
326                    result = result + (multiplier * posVal);
327
328                    if ((i - 1) >= limit)
329                        multiplier = multiplier * radix;
330                }
331            }
332
333            if (value[0] == '-')     // negative values
334            {
335                if ((result.data[maxLength - 1] & 0x80000000) == 0)
336                    throw (new ArithmeticException("Negative underflow in constructor."));
337            }
338            else    // positive values
339            {
340                if ((result.data[maxLength - 1] & 0x80000000) != 0)
341                    throw (new ArithmeticException("Positive overflow in constructor."));
342            }
343
344            data = new uint[maxLength];
345            for (int i = 0; i < result.dataLength; i++)
346                data[i] = result.data[i];
347
348            dataLength = result.dataLength;
349        }
350
351
352        //***********************************************************************
353        // Constructor (Default value provided by an array of bytes)
354        //
355        // The lowest index of the input byte array (i.e [0]) should contain the
356        // most significant byte of the number, and the highest index should
357        // contain the least significant byte.
358        //
359        // E.g.
360        // To initialize "a" with the default value of 0x1D4F in base 16
361        //      byte[] temp = { 0x1D, 0x4F };
362        //      BigInteger a = new BigInteger(temp)
363        //
364        // Note that this method of initialization does not allow the
365        // sign to be specified.
366        //
367        //***********************************************************************
368
369        public BigInteger(byte[] inData)
370        {
371            setBytes(inData);
372        }
373
374        /**         
375         * Set the internal byte Array of this BigInteger to the given One
376         * acts like the Constructor
377         */
378        public void setBytes(byte[] inData)
379        {
380            dataLength = inData.Length >> 2;
381
382            int leftOver = inData.Length & 0x3;
383            if (leftOver != 0)         // length not multiples of 4
384                dataLength++;
385
386
387            if (dataLength > maxLength)
388                throw (new ArithmeticException("Byte overflow in setBytes method."));
389
390            data = new uint[maxLength];
391
392            for (int i = inData.Length - 1, j = 0; i >= 3; i -= 4, j++)
393            {
394                data[j] = (uint)((inData[i - 3] << 24) + (inData[i - 2] << 16) +
395                                 (inData[i - 1] << 8) + inData[i]);
396            }
397
398            if (leftOver == 1)
399                data[dataLength - 1] = (uint)inData[0];
400            else if (leftOver == 2)
401                data[dataLength - 1] = (uint)((inData[0] << 8) + inData[1]);
402            else if (leftOver == 3)
403                data[dataLength - 1] = (uint)((inData[0] << 16) + (inData[1] << 8) + inData[2]);
404
405
406            while (dataLength > 1 && data[dataLength - 1] == 0)
407                dataLength--;
408        }
409
410        //***********************************************************************
411        // Constructor (Default value provided by an array of bytes of the
412        // specified length.)
413        //***********************************************************************
414
415        public BigInteger(byte[] inData, int inLen)
416        {
417            dataLength = inLen >> 2;
418
419            int leftOver = inLen & 0x3;
420            if (leftOver != 0)         // length not multiples of 4
421                dataLength++;
422
423            if (dataLength > maxLength || inLen > inData.Length)
424                throw (new ArithmeticException("Byte overflow in constructor."));
425
426
427            data = new uint[maxLength];
428
429            for (int i = inLen - 1, j = 0; i >= 3; i -= 4, j++)
430            {
431                data[j] = (uint)((inData[i - 3] << 24) + (inData[i - 2] << 16) +
432                                 (inData[i - 1] << 8) + inData[i]);
433            }
434
435            if (leftOver == 1)
436                data[dataLength - 1] = (uint)inData[0];
437            else if (leftOver == 2)
438                data[dataLength - 1] = (uint)((inData[0] << 8) + inData[1]);
439            else if (leftOver == 3)
440                data[dataLength - 1] = (uint)((inData[0] << 16) + (inData[1] << 8) + inData[2]);
441
442
443            if (dataLength == 0)
444                dataLength = 1;
445
446            while (dataLength > 1 && data[dataLength - 1] == 0)
447                dataLength--;
448
449            //Console.WriteLine("Len = " + dataLength);
450        }
451
452
453        //***********************************************************************
454        // Constructor (Default value provided by an array of unsigned integers)
455        //*********************************************************************
456
457        public BigInteger(uint[] inData)
458        {
459            dataLength = inData.Length;
460
461            if (dataLength > maxLength)
462                throw (new ArithmeticException("Byte overflow in constructor."));
463
464            data = new uint[maxLength];
465
466            for (int i = dataLength - 1, j = 0; i >= 0; i--, j++)
467                data[j] = inData[i];
468
469            while (dataLength > 1 && data[dataLength - 1] == 0)
470                dataLength--;
471
472            //Console.WriteLine("Len = " + dataLength);
473        }
474
475
476        //***********************************************************************
477        // Overloading of the typecast operator.
478        // For BigInteger bi = 10;
479        //***********************************************************************
480
481        public static implicit operator BigInteger(long value)
482        {
483            return (new BigInteger(value));
484        }
485
486        public static implicit operator BigInteger(ulong value)
487        {
488            return (new BigInteger(value));
489        }
490
491        public static implicit operator BigInteger(int value)
492        {
493            return (new BigInteger((long)value));
494        }
495
496        public static implicit operator BigInteger(uint value)
497        {
498            return (new BigInteger((ulong)value));
499        }
500
501
502        //***********************************************************************
503        // Overloading of addition operator
504        //***********************************************************************
505
506        public static BigInteger operator +(BigInteger bi1, BigInteger bi2)
507        {
508            BigInteger result = new BigInteger();
509
510            result.dataLength = (bi1.dataLength > bi2.dataLength) ? bi1.dataLength : bi2.dataLength;
511
512            long carry = 0;
513            for (int i = 0; i < result.dataLength; i++)
514            {
515                long sum = (long)bi1.data[i] + (long)bi2.data[i] + carry;
516                carry = sum >> 32;
517                result.data[i] = (uint)(sum & 0xFFFFFFFF);
518            }
519
520            if (carry != 0 && result.dataLength < maxLength)
521            {
522                result.data[result.dataLength] = (uint)(carry);
523                result.dataLength++;
524            }
525
526            while (result.dataLength > 1 && result.data[result.dataLength - 1] == 0)
527                result.dataLength--;
528
529
530            // overflow check
531            int lastPos = maxLength - 1;
532            if ((bi1.data[lastPos] & 0x80000000) == (bi2.data[lastPos] & 0x80000000) &&
533               (result.data[lastPos] & 0x80000000) != (bi1.data[lastPos] & 0x80000000))
534            {
535                throw (new ArithmeticException());
536            }
537
538            return result;
539        }
540
541
542        //***********************************************************************
543        // Overloading of the unary ++ operator
544        //***********************************************************************
545
546        public static BigInteger operator ++(BigInteger bi1)
547        {
548            BigInteger result = new BigInteger(bi1);
549
550            long val, carry = 1;
551            int index = 0;
552
553            while (carry != 0 && index < maxLength)
554            {
555                val = (long)(result.data[index]);
556                val++;
557
558                result.data[index] = (uint)(val & 0xFFFFFFFF);
559                carry = val >> 32;
560
561                index++;
562            }
563
564            if (index > result.dataLength)
565                result.dataLength = index;
566            else
567            {
568                while (result.dataLength > 1 && result.data[result.dataLength - 1] == 0)
569                    result.dataLength--;
570            }
571
572            // overflow check
573            int lastPos = maxLength - 1;
574
575            // overflow if initial value was +ve but ++ caused a sign
576            // change to negative.
577
578            if ((bi1.data[lastPos] & 0x80000000) == 0 &&
579               (result.data[lastPos] & 0x80000000) != (bi1.data[lastPos] & 0x80000000))
580            {
581                throw (new ArithmeticException("Overflow in ++."));
582            }
583            return result;
584        }
585
586
587        //***********************************************************************
588        // Overloading of subtraction operator
589        //***********************************************************************
590
591        public static BigInteger operator -(BigInteger bi1, BigInteger bi2)
592        {
593            BigInteger result = new BigInteger();
594
595            result.dataLength = (bi1.dataLength > bi2.dataLength) ? bi1.dataLength : bi2.dataLength;
596
597            long carryIn = 0;
598            for (int i = 0; i < result.dataLength; i++)
599            {
600                long diff;
601
602                diff = (long)bi1.data[i] - (long)bi2.data[i] - carryIn;
603                result.data[i] = (uint)(diff & 0xFFFFFFFF);
604
605                if (diff < 0)
606                    carryIn = 1;
607                else
608                    carryIn = 0;
609            }
610
611            // roll over to negative
612            if (carryIn != 0)
613            {
614                for (int i = result.dataLength; i < maxLength; i++)
615                    result.data[i] = 0xFFFFFFFF;
616                result.dataLength = maxLength;
617            }
618
619            // fixed in v1.03 to give correct datalength for a - (-b)
620            while (result.dataLength > 1 && result.data[result.dataLength - 1] == 0)
621                result.dataLength--;
622
623            // overflow check
624
625            int lastPos = maxLength - 1;
626            if ((bi1.data[lastPos] & 0x80000000) != (bi2.data[lastPos] & 0x80000000) &&
627               (result.data[lastPos] & 0x80000000) != (bi1.data[lastPos] & 0x80000000))
628            {
629                throw (new ArithmeticException());
630            }
631
632            return result;
633        }
634
635
636        //***********************************************************************
637        // Overloading of the unary -- operator
638        //***********************************************************************
639
640        public static BigInteger operator --(BigInteger bi1)
641        {
642            BigInteger result = new BigInteger(bi1);
643
644            long val;
645            bool carryIn = true;
646            int index = 0;
647
648            while (carryIn && index < maxLength)
649            {
650                val = (long)(result.data[index]);
651                val--;
652
653                result.data[index] = (uint)(val & 0xFFFFFFFF);
654
655                if (val >= 0)
656                    carryIn = false;
657
658                index++;
659            }
660
661            if (index > result.dataLength)
662                result.dataLength = index;
663
664            while (result.dataLength > 1 && result.data[result.dataLength - 1] == 0)
665                result.dataLength--;
666
667            // overflow check
668            int lastPos = maxLength - 1;
669
670            // overflow if initial value was -ve but -- caused a sign
671            // change to positive.
672
673            if ((bi1.data[lastPos] & 0x80000000) != 0 &&
674               (result.data[lastPos] & 0x80000000) != (bi1.data[lastPos] & 0x80000000))
675            {
676                throw (new ArithmeticException("Underflow in --."));
677            }
678
679            return result;
680        }
681
682
683        //***********************************************************************
684        // Overloading of multiplication operator
685        //***********************************************************************
686
687        public static BigInteger operator *(BigInteger bi1, BigInteger bi2)
688        {
689            int lastPos = maxLength - 1;
690            bool bi1Neg = false, bi2Neg = false;
691
692            // take the absolute value of the inputs
693            try
694            {
695                if ((bi1.data[lastPos] & 0x80000000) != 0)     // bi1 negative
696                {
697                    bi1Neg = true; bi1 = -bi1;
698                }
699                if ((bi2.data[lastPos] & 0x80000000) != 0)     // bi2 negative
700                {
701                    bi2Neg = true; bi2 = -bi2;
702                }
703            }
704            catch (Exception) { }
705
706            BigInteger result = new BigInteger();
707
708            // multiply the absolute values
709            try
710            {
711                for (int i = 0; i < bi1.dataLength; i++)
712                {
713                    if (bi1.data[i] == 0) continue;
714
715                    ulong mcarry = 0;
716                    for (int j = 0, k = i; j < bi2.dataLength; j++, k++)
717                    {
718                        // k = i + j
719                        ulong val = ((ulong)bi1.data[i] * (ulong)bi2.data[j]) +
720                                     (ulong)result.data[k] + mcarry;
721
722                        result.data[k] = (uint)(val & 0xFFFFFFFF);
723                        mcarry = (val >> 32);
724                    }
725
726                    if (mcarry != 0)
727                        result.data[i + bi2.dataLength] = (uint)mcarry;
728                }
729            }
730            catch (Exception)
731            {
732                throw (new ArithmeticException("Multiplication overflow."));
733            }
734
735
736            result.dataLength = bi1.dataLength + bi2.dataLength;
737            if (result.dataLength > maxLength)
738                result.dataLength = maxLength;
739
740            while (result.dataLength > 1 && result.data[result.dataLength - 1] == 0)
741                result.dataLength--;
742
743            // overflow check (result is -ve)
744            if ((result.data[lastPos] & 0x80000000) != 0)
745            {
746                if (bi1Neg != bi2Neg && result.data[lastPos] == 0x80000000)    // different sign
747                {
748                    // handle the special case where multiplication produces
749                    // a max negative number in 2's complement.
750
751                    if (result.dataLength == 1)
752                        return result;
753                    else
754                    {
755                        bool isMaxNeg = true;
756                        for (int i = 0; i < result.dataLength - 1 && isMaxNeg; i++)
757                        {
758                            if (result.data[i] != 0)
759                                isMaxNeg = false;
760                        }
761
762                        if (isMaxNeg)
763                            return result;
764                    }
765                }
766
767                throw (new ArithmeticException("Multiplication overflow."));
768            }
769
770            // if input has different signs, then result is -ve
771            if (bi1Neg != bi2Neg)
772                return -result;
773
774            return result;
775        }
776
777
778
779        //***********************************************************************
780        // Overloading of unary << operators
781        //***********************************************************************
782
783        public static BigInteger operator <<(BigInteger bi1, int shiftVal)
784        {
785            BigInteger result = new BigInteger(bi1);
786            result.dataLength = shiftLeft(result.data, shiftVal);
787
788            return result;
789        }
790
791
792        // least significant bits at lower part of buffer
793
794        private static int shiftLeft(uint[] buffer, int shiftVal)
795        {
796            int shiftAmount = 32;
797            int bufLen = buffer.Length;
798
799            while (bufLen > 1 && buffer[bufLen - 1] == 0)
800                bufLen--;
801
802            for (int count = shiftVal; count > 0; )
803            {
804                if (count < shiftAmount)
805                    shiftAmount = count;
806
807                //Console.WriteLine("shiftAmount = {0}", shiftAmount);
808
809                ulong carry = 0;
810                for (int i = 0; i < bufLen; i++)
811                {
812                    ulong val = ((ulong)buffer[i]) << shiftAmount;
813                    val |= carry;
814
815                    buffer[i] = (uint)(val & 0xFFFFFFFF);
816                    carry = val >> 32;
817                }
818
819                if (carry != 0)
820                {
821                    if (bufLen + 1 <= buffer.Length)
822                    {
823                        buffer[bufLen] = (uint)carry;
824                        bufLen++;
825                    }
826                }
827                count -= shiftAmount;
828            }
829            return bufLen;
830        }
831
832
833        //***********************************************************************
834        // Overloading of unary >> operators
835        //***********************************************************************
836
837        public static BigInteger operator >>(BigInteger bi1, int shiftVal)
838        {
839            BigInteger result = new BigInteger(bi1);
840            result.dataLength = shiftRight(result.data, shiftVal);
841
842
843            if ((bi1.data[maxLength - 1] & 0x80000000) != 0) // negative
844            {
845                for (int i = maxLength - 1; i >= result.dataLength; i--)
846                    result.data[i] = 0xFFFFFFFF;
847
848                uint mask = 0x80000000;
849                for (int i = 0; i < 32; i++)
850                {
851                    if ((result.data[result.dataLength - 1] & mask) != 0)
852                        break;
853
854                    result.data[result.dataLength - 1] |= mask;
855                    mask >>= 1;
856                }
857                result.dataLength = maxLength;
858            }
859
860            return result;
861        }
862
863
864        private static int shiftRight(uint[] buffer, int shiftVal)
865        {
866            int shiftAmount = 32;
867            int invShift = 0;
868            int bufLen = buffer.Length;
869
870            while (bufLen > 1 && buffer[bufLen - 1] == 0)
871                bufLen--;
872
873            //Console.WriteLine("bufLen = " + bufLen + " buffer.Length = " + buffer.Length);
874
875            for (int count = shiftVal; count > 0; )
876            {
877                if (count < shiftAmount)
878                {
879                    shiftAmount = count;
880                    invShift = 32 - shiftAmount;
881                }
882
883                //Console.WriteLine("shiftAmount = {0}", shiftAmount);
884
885                ulong carry = 0;
886                for (int i = bufLen - 1; i >= 0; i--)
887                {
888                    ulong val = ((ulong)buffer[i]) >> shiftAmount;
889                    val |= carry;
890
891                    carry = ((ulong)buffer[i]) << invShift;
892                    buffer[i] = (uint)(val);
893                }
894
895                count -= shiftAmount;
896            }
897
898            while (bufLen > 1 && buffer[bufLen - 1] == 0)
899                bufLen--;
900
901            return bufLen;
902        }
903
904
905        //***********************************************************************
906        // Overloading of the NOT operator (1's complement)
907        //***********************************************************************
908
909        public static BigInteger operator ~(BigInteger bi1)
910        {
911            BigInteger result = new BigInteger(bi1);
912
913            for (int i = 0; i < maxLength; i++)
914                result.data[i] = (uint)(~(bi1.data[i]));
915
916            result.dataLength = maxLength;
917
918            while (result.dataLength > 1 && result.data[result.dataLength - 1] == 0)
919                result.dataLength--;
920
921            return result;
922        }
923
924
925        //***********************************************************************
926        // Overloading of the NEGATE operator (2's complement)
927        //***********************************************************************
928
929        public static BigInteger operator -(BigInteger bi1)
930        {
931            // handle neg of zero separately since it'll cause an overflow
932            // if we proceed.
933
934            if (bi1.dataLength == 1 && bi1.data[0] == 0)
935                return (new BigInteger());
936
937            BigInteger result = new BigInteger(bi1);
938
939            // 1's complement
940            for (int i = 0; i < maxLength; i++)
941                result.data[i] = (uint)(~(bi1.data[i]));
942
943            // add one to result of 1's complement
944            long val, carry = 1;
945            int index = 0;
946
947            while (carry != 0 && index < maxLength)
948            {
949                val = (long)(result.data[index]);
950                val++;
951
952                result.data[index] = (uint)(val & 0xFFFFFFFF);
953                carry = val >> 32;
954
955                index++;
956            }
957
958            if ((bi1.data[maxLength - 1] & 0x80000000) == (result.data[maxLength - 1] & 0x80000000))
959                throw (new ArithmeticException("Overflow in negation.\n"));
960
961            result.dataLength = maxLength;
962
963            while (result.dataLength > 1 && result.data[result.dataLength - 1] == 0)
964                result.dataLength--;
965            return result;
966        }
967
968
969        //***********************************************************************
970        // Overloading of equality operator
971        //***********************************************************************
972
973        public static bool operator ==(BigInteger bi1, BigInteger bi2)
974        {
975            // same instance?
976            if (System.Object.ReferenceEquals(bi1, bi2))
977                return true;
978
979            // null?
980            if ((object)bi1 == null || (object)bi2 == null)
981                return false;
982
983            return bi1.Equals(bi2);
984        }
985
986
987        public static bool operator !=(BigInteger bi1, BigInteger bi2)
988        {
989            if (System.Object.ReferenceEquals(bi1, bi2))
990                return false;
991
992            if ((object)bi1 == null || (object)bi2 == null)
993                return true;
994
995            return !(bi1.Equals(bi2));
996        }
997
998
999        public override bool Equals(object o)
1000        {
1001            BigInteger bi = o as BigInteger;
1002            if (bi == null)
1003                return false;
1004
1005            if (this.dataLength != bi.dataLength)
1006                return false;
1007
1008            for (int i = 0; i < this.dataLength; i++)
1009            {
1010                if (this.data[i] != bi.data[i])
1011                    return false;
1012            }
1013            return true;
1014        }
1015
1016
1017        public override int GetHashCode()
1018        {
1019            return this.ToString().GetHashCode();
1020        }
1021
1022
1023        //***********************************************************************
1024        // Overloading of inequality operator
1025        //***********************************************************************
1026
1027        public static bool operator >(BigInteger bi1, BigInteger bi2)
1028        {
1029            int pos = maxLength - 1;
1030
1031            // bi1 is negative, bi2 is positive
1032            if ((bi1.data[pos] & 0x80000000) != 0 && (bi2.data[pos] & 0x80000000) == 0)
1033                return false;
1034
1035                // bi1 is positive, bi2 is negative
1036            else if ((bi1.data[pos] & 0x80000000) == 0 && (bi2.data[pos] & 0x80000000) != 0)
1037                return true;
1038
1039            // same sign
1040            int len = (bi1.dataLength > bi2.dataLength) ? bi1.dataLength : bi2.dataLength;
1041            for (pos = len - 1; pos >= 0 && bi1.data[pos] == bi2.data[pos]; pos--) ;
1042
1043            if (pos >= 0)
1044            {
1045                if (bi1.data[pos] > bi2.data[pos])
1046                    return true;
1047                return false;
1048            }
1049            return false;
1050        }
1051
1052
1053        public static bool operator <(BigInteger bi1, BigInteger bi2)
1054        {
1055            int pos = maxLength - 1;
1056
1057            // bi1 is negative, bi2 is positive
1058            if ((bi1.data[pos] & 0x80000000) != 0 && (bi2.data[pos] & 0x80000000) == 0)
1059                return true;
1060
1061                // bi1 is positive, bi2 is negative
1062            else if ((bi1.data[pos] & 0x80000000) == 0 && (bi2.data[pos] & 0x80000000) != 0)
1063                return false;
1064
1065            // same sign
1066            int len = (bi1.dataLength > bi2.dataLength) ? bi1.dataLength : bi2.dataLength;
1067            for (pos = len - 1; pos >= 0 && bi1.data[pos] == bi2.data[pos]; pos--) ;
1068
1069            if (pos >= 0)
1070            {
1071                if (bi1.data[pos] < bi2.data[pos])
1072                    return true;
1073                return false;
1074            }
1075            return false;
1076        }
1077
1078
1079        public static bool operator >=(BigInteger bi1, BigInteger bi2)
1080        {
1081            return (bi1 == bi2 || bi1 > bi2);
1082        }
1083
1084
1085        public static bool operator <=(BigInteger bi1, BigInteger bi2)
1086        {
1087            return (bi1 == bi2 || bi1 < bi2);
1088        }
1089
1090
1091        //***********************************************************************
1092        // Private function that supports the division of two numbers with
1093        // a divisor that has more than 1 digit.
1094        //
1095        // Algorithm taken from [1]
1096        //***********************************************************************
1097
1098        private static void multiByteDivide(BigInteger bi1, BigInteger bi2,
1099                                            BigInteger outQuotient, BigInteger outRemainder)
1100        {
1101            uint[] result = new uint[maxLength];
1102
1103            int remainderLen = bi1.dataLength + 1;
1104            uint[] remainder = new uint[remainderLen];
1105
1106            uint mask = 0x80000000;
1107            uint val = bi2.data[bi2.dataLength - 1];
1108            int shift = 0, resultPos = 0;
1109
1110            while (mask != 0 && (val & mask) == 0)
1111            {
1112                shift++; mask >>= 1;
1113            }
1114
1115            //Console.WriteLine("shift = {0}", shift);
1116            //Console.WriteLine("Before bi1 Len = {0}, bi2 Len = {1}", bi1.dataLength, bi2.dataLength);
1117
1118            for (int i = 0; i < bi1.dataLength; i++)
1119                remainder[i] = bi1.data[i];
1120            shiftLeft(remainder, shift);
1121            bi2 = bi2 << shift;
1122
1123            /*
1124            Console.WriteLine("bi1 Len = {0}, bi2 Len = {1}", bi1.dataLength, bi2.dataLength);
1125            Console.WriteLine("dividend = " + bi1 + "\ndivisor = " + bi2);
1126            for(int q = remainderLen - 1; q >= 0; q--)
1127                    Console.Write("{0:x2}", remainder[q]);
1128            Console.WriteLine();
1129            */
1130
1131            int j = remainderLen - bi2.dataLength;
1132            int pos = remainderLen - 1;
1133
1134            ulong firstDivisorByte = bi2.data[bi2.dataLength - 1];
1135            ulong secondDivisorByte = bi2.data[bi2.dataLength - 2];
1136
1137            int divisorLen = bi2.dataLength + 1;
1138            uint[] dividendPart = new uint[divisorLen];
1139
1140            while (j > 0)
1141            {
1142                ulong dividend = ((ulong)remainder[pos] << 32) + (ulong)remainder[pos - 1];
1143                //Console.WriteLine("dividend = {0}", dividend);
1144
1145                ulong q_hat = dividend / firstDivisorByte;
1146                ulong r_hat = dividend % firstDivisorByte;
1147
1148                //Console.WriteLine("q_hat = {0:X}, r_hat = {1:X}", q_hat, r_hat);
1149
1150                bool done = false;
1151                while (!done)
1152                {
1153                    done = true;
1154
1155                    if (q_hat == 0x100000000 ||
1156                       (q_hat * secondDivisorByte) > ((r_hat << 32) + remainder[pos - 2]))
1157                    {
1158                        q_hat--;
1159                        r_hat += firstDivisorByte;
1160
1161                        if (r_hat < 0x100000000)
1162                            done = false;
1163                    }
1164                }
1165
1166                for (int h = 0; h < divisorLen; h++)
1167                    dividendPart[h] = remainder[pos - h];
1168
1169                BigInteger kk = new BigInteger(dividendPart);
1170                BigInteger ss = bi2 * (long)q_hat;
1171
1172                //Console.WriteLine("ss before = " + ss);
1173                while (ss > kk)
1174                {
1175                    q_hat--;
1176                    ss -= bi2;
1177                    //Console.WriteLine(ss);
1178                }
1179                BigInteger yy = kk - ss;
1180
1181                //Console.WriteLine("ss = " + ss);
1182                //Console.WriteLine("kk = " + kk);
1183                //Console.WriteLine("yy = " + yy);
1184
1185                for (int h = 0; h < divisorLen; h++)
1186                    remainder[pos - h] = yy.data[bi2.dataLength - h];
1187
1188                /*
1189                Console.WriteLine("dividend = ");
1190                for(int q = remainderLen - 1; q >= 0; q--)
1191                        Console.Write("{0:x2}", remainder[q]);
1192                Console.WriteLine("\n************ q_hat = {0:X}\n", q_hat);
1193                */
1194
1195                result[resultPos++] = (uint)q_hat;
1196
1197                pos--;
1198                j--;
1199            }
1200
1201            outQuotient.dataLength = resultPos;
1202            int y = 0;
1203            for (int x = outQuotient.dataLength - 1; x >= 0; x--, y++)
1204                outQuotient.data[y] = result[x];
1205            for (; y < maxLength; y++)
1206                outQuotient.data[y] = 0;
1207
1208            while (outQuotient.dataLength > 1 && outQuotient.data[outQuotient.dataLength - 1] == 0)
1209                outQuotient.dataLength--;
1210
1211            if (outQuotient.dataLength == 0)
1212                outQuotient.dataLength = 1;
1213
1214            outRemainder.dataLength = shiftRight(remainder, shift);
1215
1216            for (y = 0; y < outRemainder.dataLength; y++)
1217                outRemainder.data[y] = remainder[y];
1218            for (; y < maxLength; y++)
1219                outRemainder.data[y] = 0;
1220        }
1221
1222
1223        //***********************************************************************
1224        // Private function that supports the division of two numbers with
1225        // a divisor that has only 1 digit.
1226        //***********************************************************************
1227
1228        private static void singleByteDivide(BigInteger bi1, BigInteger bi2,
1229                                             BigInteger outQuotient, BigInteger outRemainder)
1230        {
1231            uint[] result = new uint[maxLength];
1232            int resultPos = 0;
1233
1234            // copy dividend to reminder
1235            for (int i = 0; i < maxLength; i++)
1236                outRemainder.data[i] = bi1.data[i];
1237            outRemainder.dataLength = bi1.dataLength;
1238
1239            while (outRemainder.dataLength > 1 && outRemainder.data[outRemainder.dataLength - 1] == 0)
1240                outRemainder.dataLength--;
1241
1242            ulong divisor = (ulong)bi2.data[0];
1243            int pos = outRemainder.dataLength - 1;
1244            ulong dividend = (ulong)outRemainder.data[pos];
1245
1246            //Console.WriteLine("divisor = " + divisor + " dividend = " + dividend);
1247            //Console.WriteLine("divisor = " + bi2 + "\ndividend = " + bi1);
1248
1249            if (dividend >= divisor)
1250            {
1251                ulong quotient = dividend / divisor;
1252                result[resultPos++] = (uint)quotient;
1253
1254                outRemainder.data[pos] = (uint)(dividend % divisor);
1255            }
1256            pos--;
1257
1258            while (pos >= 0)
1259            {
1260                //Console.WriteLine(pos);
1261
1262                dividend = ((ulong)outRemainder.data[pos + 1] << 32) + (ulong)outRemainder.data[pos];
1263                ulong quotient = dividend / divisor;
1264                result[resultPos++] = (uint)quotient;
1265
1266                outRemainder.data[pos + 1] = 0;
1267                outRemainder.data[pos--] = (uint)(dividend % divisor);
1268                //Console.WriteLine(">>>> " + bi1);
1269            }
1270
1271            outQuotient.dataLength = resultPos;
1272            int j = 0;
1273            for (int i = outQuotient.dataLength - 1; i >= 0; i--, j++)
1274                outQuotient.data[j] = result[i];
1275            for (; j < maxLength; j++)
1276                outQuotient.data[j] = 0;
1277
1278            while (outQuotient.dataLength > 1 && outQuotient.data[outQuotient.dataLength - 1] == 0)
1279                outQuotient.dataLength--;
1280
1281            if (outQuotient.dataLength == 0)
1282                outQuotient.dataLength = 1;
1283
1284            while (outRemainder.dataLength > 1 && outRemainder.data[outRemainder.dataLength - 1] == 0)
1285                outRemainder.dataLength--;
1286        }
1287
1288
1289        //***********************************************************************
1290        // Overloading of division operator
1291        //***********************************************************************
1292
1293        public static BigInteger operator /(BigInteger bi1, BigInteger bi2)
1294        {
1295            BigInteger quotient = new BigInteger();
1296            BigInteger remainder = new BigInteger();
1297
1298            int lastPos = maxLength - 1;
1299            bool divisorNeg = false, dividendNeg = false;
1300
1301            if ((bi1.data[lastPos] & 0x80000000) != 0)     // bi1 negative
1302            {
1303                bi1 = -bi1;
1304                dividendNeg = true;
1305            }
1306            if ((bi2.data[lastPos] & 0x80000000) != 0)     // bi2 negative
1307            {
1308                bi2 = -bi2;
1309                divisorNeg = true;
1310            }
1311
1312            if (bi1 < bi2)
1313            {
1314                return quotient;
1315            }
1316
1317            else
1318            {
1319                if (bi2.dataLength == 1)
1320                    singleByteDivide(bi1, bi2, quotient, remainder);
1321                else
1322                    multiByteDivide(bi1, bi2, quotient, remainder);
1323
1324                if (dividendNeg != divisorNeg)
1325                    return -quotient;
1326
1327                return quotient;
1328            }
1329        }
1330
1331
1332        //***********************************************************************
1333        // Overloading of modulus operator
1334        //***********************************************************************
1335
1336        public static BigInteger operator %(BigInteger bi1, BigInteger bi2)
1337        {
1338            BigInteger quotient = new BigInteger();
1339            BigInteger remainder = new BigInteger(bi1);
1340
1341            int lastPos = maxLength - 1;
1342            bool dividendNeg = false;
1343
1344            if ((bi1.data[lastPos] & 0x80000000) != 0)     // bi1 negative
1345            {
1346                bi1 = -bi1;
1347                dividendNeg = true;
1348            }
1349            if ((bi2.data[lastPos] & 0x80000000) != 0)     // bi2 negative
1350                bi2 = -bi2;
1351
1352            if (bi1 < bi2)
1353            {
1354                return remainder;
1355            }
1356
1357            else
1358            {
1359                if (bi2.dataLength == 1)
1360                    singleByteDivide(bi1, bi2, quotient, remainder);
1361                else
1362                    multiByteDivide(bi1, bi2, quotient, remainder);
1363
1364                if (dividendNeg)
1365                    return -remainder;
1366
1367                return remainder;
1368            }
1369        }
1370
1371
1372        //***********************************************************************
1373        // Overloading of bitwise AND operator
1374        //***********************************************************************
1375
1376        public static BigInteger operator &(BigInteger bi1, BigInteger bi2)
1377        {
1378            BigInteger result = new BigInteger();
1379
1380            int len = (bi1.dataLength > bi2.dataLength) ? bi1.dataLength : bi2.dataLength;
1381
1382            for (int i = 0; i < len; i++)
1383            {
1384                uint sum = (uint)(bi1.data[i] & bi2.data[i]);
1385                result.data[i] = sum;
1386            }
1387
1388            result.dataLength = maxLength;
1389
1390            while (result.dataLength > 1 && result.data[result.dataLength - 1] == 0)
1391                result.dataLength--;
1392
1393            return result;
1394        }
1395
1396
1397        //***********************************************************************
1398        // Overloading of bitwise OR operator
1399        //***********************************************************************
1400
1401        public static BigInteger operator |(BigInteger bi1, BigInteger bi2)
1402        {
1403            BigInteger result = new BigInteger();
1404
1405            int len = (bi1.dataLength > bi2.dataLength) ? bi1.dataLength : bi2.dataLength;
1406
1407            for (int i = 0; i < len; i++)
1408            {
1409                uint sum = (uint)(bi1.data[i] | bi2.data[i]);
1410                result.data[i] = sum;
1411            }
1412
1413            result.dataLength = maxLength;
1414
1415            while (result.dataLength > 1 && result.data[result.dataLength - 1] == 0)
1416                result.dataLength--;
1417
1418            return result;
1419        }
1420
1421
1422        //***********************************************************************
1423        // Overloading of bitwise XOR operator
1424        //***********************************************************************
1425
1426        public static BigInteger operator ^(BigInteger bi1, BigInteger bi2)
1427        {
1428            BigInteger result = new BigInteger();
1429
1430            int len = (bi1.dataLength > bi2.dataLength) ? bi1.dataLength : bi2.dataLength;
1431
1432            for (int i = 0; i < len; i++)
1433            {
1434                uint sum = (uint)(bi1.data[i] ^ bi2.data[i]);
1435                result.data[i] = sum;
1436            }
1437
1438            result.dataLength = maxLength;
1439
1440            while (result.dataLength > 1 && result.data[result.dataLength - 1] == 0)
1441                result.dataLength--;
1442
1443            return result;
1444        }
1445
1446
1447        //***********************************************************************
1448        // Returns max(this, bi)
1449        //***********************************************************************
1450
1451        public BigInteger max(BigInteger bi)
1452        {
1453            if (this > bi)
1454                return (new BigInteger(this));
1455            else
1456                return (new BigInteger(bi));
1457        }
1458
1459
1460        //***********************************************************************
1461        // Returns min(this, bi)
1462        //***********************************************************************
1463
1464        public BigInteger min(BigInteger bi)
1465        {
1466            if (this < bi)
1467                return (new BigInteger(this));
1468            else
1469                return (new BigInteger(bi));
1470
1471        }
1472
1473
1474        //***********************************************************************
1475        // Returns the absolute value
1476        //***********************************************************************
1477
1478        public BigInteger abs()
1479        {
1480            if ((this.data[maxLength - 1] & 0x80000000) != 0)
1481                return (-this);
1482            else
1483                return (new BigInteger(this));
1484        }
1485
1486
1487        //***********************************************************************
1488        // Returns a string representing the BigInteger in base 10.
1489        //***********************************************************************
1490
1491        public override string ToString()
1492        {
1493            return ToString(10);
1494        }
1495
1496
1497        //***********************************************************************
1498        // Returns a string representing the BigInteger in sign-and-magnitude
1499        // format in the specified radix.
1500        //
1501        // Example
1502        // -------
1503        // If the value of BigInteger is -255 in base 10, then
1504        // ToString(16) returns "-FF"
1505        //
1506        //***********************************************************************
1507
1508        public string ToString(int radix)
1509        {
1510            if (radix < 2 || radix > 36)
1511                throw (new ArgumentException("Radix must be >= 2 and <= 36"));
1512
1513            string charSet = "ABCDEFGHIJKLMNOPQRSTUVWXYZ";
1514            string result = "";
1515
1516            BigInteger a = this;
1517
1518            bool negative = false;
1519            if ((a.data[maxLength - 1] & 0x80000000) != 0)
1520            {
1521                negative = true;
1522                try
1523                {
1524                    a = -a;
1525                }
1526                catch (Exception) { }
1527            }
1528
1529            BigInteger quotient = new BigInteger();
1530            BigInteger remainder = new BigInteger();
1531            BigInteger biRadix = new BigInteger(radix);
1532
1533            if (a.dataLength == 1 && a.data[0] == 0)
1534                result = "0";
1535            else
1536            {
1537                while (a.dataLength > 1 || (a.dataLength == 1 && a.data[0] != 0))
1538                {
1539                    singleByteDivide(a, biRadix, quotient, remainder);
1540
1541                    if (remainder.data[0] < 10)
1542                        result = remainder.data[0] + result;
1543                    else
1544                        result = charSet[(int)remainder.data[0] - 10] + result;
1545
1546                    a = quotient;
1547                }
1548                if (negative)
1549                    result = "-" + result;
1550            }
1551
1552            return result;
1553        }
1554
1555
1556        //***********************************************************************
1557        // Returns a hex string showing the contains of the BigInteger
1558        //
1559        // Examples
1560        // -------
1561        // 1) If the value of BigInteger is 255 in base 10, then
1562        //    ToHexString() returns "FF"
1563        //
1564        // 2) If the value of BigInteger is -255 in base 10, then
1565        //    ToHexString() returns ".....FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF01",
1566        //    which is the 2's complement representation of -255.
1567        //
1568        //***********************************************************************
1569
1570        public string ToHexString()
1571        {
1572            string result = data[dataLength - 1].ToString("X");
1573
1574            for (int i = dataLength - 2; i >= 0; i--)
1575            {
1576                result += data[i].ToString("X8");
1577            }
1578
1579            return result;
1580        }
1581
1582
1583
1584        //***********************************************************************
1585        // Modulo Exponentiation
1586        //***********************************************************************
1587
1588        public BigInteger modPow(BigInteger exp, BigInteger n)
1589        {
1590            if ((exp.data[maxLength - 1] & 0x80000000) != 0)
1591                throw (new ArithmeticException("Positive exponents only."));
1592
1593            BigInteger resultNum = 1;
1594            BigInteger tempNum;
1595            bool thisNegative = false;
1596
1597            if ((this.data[maxLength - 1] & 0x80000000) != 0)   // negative this
1598            {
1599                tempNum = -this % n;
1600                thisNegative = true;
1601            }
1602            else
1603                tempNum = this % n;  // ensures (tempNum * tempNum) < b^(2k)
1604
1605            if ((n.data[maxLength - 1] & 0x80000000) != 0)   // negative n
1606                n = -n;
1607
1608            // calculate constant = b^(2k) / m
1609            BigInteger constant = new BigInteger();
1610
1611            int i = n.dataLength << 1;
1612            constant.data[i] = 0x00000001;
1613            constant.dataLength = i + 1;
1614
1615            constant = constant / n;
1616            int totalBits = exp.bitCount();
1617            int count = 0;
1618
1619            // perform squaring and multiply exponentiation
1620            for (int pos = 0; pos < exp.dataLength; pos++)
1621            {
1622                uint mask = 0x01;
1623                //Console.WriteLine("pos = " + pos);
1624
1625                for (int index = 0; index < 32; index++)
1626                {
1627                    if ((exp.data[pos] & mask) != 0)
1628                        resultNum = BarrettReduction(resultNum * tempNum, n, constant);
1629
1630                    mask <<= 1;
1631
1632                    tempNum = BarrettReduction(tempNum * tempNum, n, constant);
1633
1634
1635                    if (tempNum.dataLength == 1 && tempNum.data[0] == 1)
1636                    {
1637                        if (thisNegative && (exp.data[0] & 0x1) != 0)    //odd exp
1638                            return -resultNum;
1639                        return resultNum;
1640                    }
1641                    count++;
1642                    if (count == totalBits)
1643                        break;
1644                }
1645            }
1646
1647            if (thisNegative && (exp.data[0] & 0x1) != 0)    //odd exp
1648                return -resultNum;
1649
1650            return resultNum;
1651        }
1652
1653
1654
1655        //***********************************************************************
1656        // Fast calculation of modular reduction using Barrett's reduction.
1657        // Requires x < b^(2k), where b is the base.  In this case, base is
1658        // 2^32 (uint).
1659        //
1660        // Reference [4]
1661        //***********************************************************************
1662
1663        private BigInteger BarrettReduction(BigInteger x, BigInteger n, BigInteger constant)
1664        {
1665            int k = n.dataLength,
1666                kPlusOne = k + 1,
1667                kMinusOne = k - 1;
1668
1669            BigInteger q1 = new BigInteger();
1670
1671            // q1 = x / b^(k-1)
1672            for (int i = kMinusOne, j = 0; i < x.dataLength; i++, j++)
1673                q1.data[j] = x.data[i];
1674            q1.dataLength = x.dataLength - kMinusOne;
1675            if (q1.dataLength <= 0)
1676                q1.dataLength = 1;
1677
1678
1679            BigInteger q2 = q1 * constant;
1680            BigInteger q3 = new BigInteger();
1681
1682            // q3 = q2 / b^(k+1)
1683            for (int i = kPlusOne, j = 0; i < q2.dataLength; i++, j++)
1684                q3.data[j] = q2.data[i];
1685            q3.dataLength = q2.dataLength - kPlusOne;
1686            if (q3.dataLength <= 0)
1687                q3.dataLength = 1;
1688
1689
1690            // r1 = x mod b^(k+1)
1691            // i.e. keep the lowest (k+1) words
1692            BigInteger r1 = new BigInteger();
1693            int lengthToCopy = (x.dataLength > kPlusOne) ? kPlusOne : x.dataLength;
1694            for (int i = 0; i < lengthToCopy; i++)
1695                r1.data[i] = x.data[i];
1696            r1.dataLength = lengthToCopy;
1697
1698
1699            // r2 = (q3 * n) mod b^(k+1)
1700            // partial multiplication of q3 and n
1701
1702            BigInteger r2 = new BigInteger();
1703            for (int i = 0; i < q3.dataLength; i++)
1704            {
1705                if (q3.data[i] == 0) continue;
1706
1707                ulong mcarry = 0;
1708                int t = i;
1709                for (int j = 0; j < n.dataLength && t < kPlusOne; j++, t++)
1710                {
1711                    // t = i + j
1712                    ulong val = ((ulong)q3.data[i] * (ulong)n.data[j]) +
1713                                 (ulong)r2.data[t] + mcarry;
1714
1715                    r2.data[t] = (uint)(val & 0xFFFFFFFF);
1716                    mcarry = (val >> 32);
1717                }
1718
1719                if (t < kPlusOne)
1720                    r2.data[t] = (uint)mcarry;
1721            }
1722            r2.dataLength = kPlusOne;
1723            while (r2.dataLength > 1 && r2.data[r2.dataLength - 1] == 0)
1724                r2.dataLength--;
1725
1726            r1 -= r2;
1727            if ((r1.data[maxLength - 1] & 0x80000000) != 0)        // negative
1728            {
1729                BigInteger val = new BigInteger();
1730                val.data[kPlusOne] = 0x00000001;
1731                val.dataLength = kPlusOne + 1;
1732                r1 += val;
1733            }
1734
1735            while (r1 >= n)
1736                r1 -= n;
1737
1738            return r1;
1739        }
1740
1741
1742        //***********************************************************************
1743        // Returns gcd(this, bi)
1744        //***********************************************************************
1745
1746        public BigInteger gcd(BigInteger bi)
1747        {
1748            BigInteger x;
1749            BigInteger y;
1750
1751            if ((data[maxLength - 1] & 0x80000000) != 0)     // negative
1752                x = -this;
1753            else
1754                x = this;
1755
1756            if ((bi.data[maxLength - 1] & 0x80000000) != 0)     // negative
1757                y = -bi;
1758            else
1759                y = bi;
1760
1761            BigInteger g = y;
1762
1763            while (x.dataLength > 1 || (x.dataLength == 1 && x.data[0] != 0))
1764            {
1765                g = x;
1766                x = y % x;
1767                y = g;
1768            }
1769
1770            return g;
1771        }
1772
1773
1774        //***********************************************************************
1775        // Populates "this" with the specified amount of random bits
1776        //***********************************************************************
1777
1778        public void genRandomBits(int bits, Random rand)
1779        {
1780            int dwords = bits >> 5;
1781            int remBits = bits & 0x1F;
1782
1783            if (remBits != 0)
1784                dwords++;
1785
1786            if (dwords > maxLength)
1787                throw (new ArithmeticException("Number of required bits > maxLength."));
1788
1789            for (int i = 0; i < dwords; i++)
1790                data[i] = (uint)(rand.NextDouble() * 0x100000000);
1791
1792            for (int i = dwords; i < maxLength; i++)
1793                data[i] = 0;
1794
1795            if (remBits != 0)
1796            {
1797                uint mask = (uint)(0x01 << (remBits - 1));
1798                data[dwords - 1] |= mask;
1799
1800                mask = (uint)(0xFFFFFFFF >> (32 - remBits));
1801                data[dwords - 1] &= mask;
1802            }
1803            else
1804                data[dwords - 1] |= 0x80000000;
1805
1806            dataLength = dwords;
1807
1808            if (dataLength == 0)
1809                dataLength = 1;
1810        }
1811
1812
1813        //***********************************************************************
1814        // Returns the position of the most significant bit in the BigInteger.
1815        //
1816        // Eg.  The result is 0, if the value of BigInteger is 0...0000 0000
1817        //      The result is 1, if the value of BigInteger is 0...0000 0001
1818        //      The result is 2, if the value of BigInteger is 0...0000 0010
1819        //      The result is 2, if the value of BigInteger is 0...0000 0011
1820        //
1821        //***********************************************************************
1822
1823        public int bitCount()
1824        {
1825            while (dataLength > 1 && data[dataLength - 1] == 0)
1826                dataLength--;
1827
1828            uint value = data[dataLength - 1];
1829            uint mask = 0x80000000;
1830            int bits = 32;
1831
1832            while (bits > 0 && (value & mask) == 0)
1833            {
1834                bits--;
1835                mask >>= 1;
1836            }
1837            bits += ((dataLength - 1) << 5);
1838
1839            return bits;
1840        }
1841
1842
1843        //***********************************************************************
1844        // Probabilistic prime test based on Fermat's little theorem
1845        //
1846        // for any a < p (p does not divide a) if
1847        //      a^(p-1) mod p != 1 then p is not prime.
1848        //
1849        // Otherwise, p is probably prime (pseudoprime to the chosen base).
1850        //
1851        // Returns
1852        // -------
1853        // True if "this" is a pseudoprime to randomly chosen
1854        // bases.  The number of chosen bases is given by the "confidence"
1855        // parameter.
1856        //
1857        // False if "this" is definitely NOT prime.
1858        //
1859        // Note - this method is fast but fails for Carmichael numbers except
1860        // when the randomly chosen base is a factor of the number.
1861        //
1862        //***********************************************************************
1863
1864        public bool FermatLittleTest(int confidence)
1865        {
1866            BigInteger thisVal;
1867            if ((this.data[maxLength - 1] & 0x80000000) != 0)        // negative
1868                thisVal = -this;
1869            else
1870                thisVal = this;
1871
1872            if (thisVal.dataLength == 1)
1873            {
1874                // test small numbers
1875                if (thisVal.data[0] == 0 || thisVal.data[0] == 1)
1876                    return false;
1877                else if (thisVal.data[0] == 2 || thisVal.data[0] == 3)
1878                    return true;
1879            }
1880
1881            if ((thisVal.data[0] & 0x1) == 0)     // even numbers
1882                return false;
1883
1884            int bits = thisVal.bitCount();
1885            BigInteger a = new BigInteger();
1886            BigInteger p_sub1 = thisVal - (new BigInteger(1));
1887            Random rand = new Random();
1888
1889            for (int round = 0; round < confidence; round++)
1890            {
1891                bool done = false;
1892
1893                while (!done)           // generate a < n
1894                {
1895                    int testBits = 0;
1896
1897                    // make sure "a" has at least 2 bits
1898                    while (testBits < 2)
1899                        testBits = (int)(rand.NextDouble() * bits);
1900
1901                    a.genRandomBits(testBits, rand);
1902
1903                    int byteLen = a.dataLength;
1904
1905                    // make sure "a" is not 0
1906                    if (byteLen > 1 || (byteLen == 1 && a.data[0] != 1))
1907                        done = true;
1908                }
1909
1910                // check whether a factor exists (fix for version 1.03)
1911                BigInteger gcdTest = a.gcd(thisVal);
1912                if (gcdTest.dataLength == 1 && gcdTest.data[0] != 1)
1913                    return false;
1914
1915                // calculate a^(p-1) mod p
1916                BigInteger expResult = a.modPow(p_sub1, thisVal);
1917
1918                int resultLen = expResult.dataLength;
1919
1920                // is NOT prime is a^(p-1) mod p != 1
1921
1922                if (resultLen > 1 || (resultLen == 1 && expResult.data[0] != 1))
1923                {
1924                    //Console.WriteLine("a = " + a.ToString());
1925                    return false;
1926                }
1927            }
1928
1929            return true;
1930        }
1931
1932
1933        //***********************************************************************
1934        // Probabilistic prime test based on Rabin-Miller's
1935        //
1936        // for any p > 0 with p - 1 = 2^s * t
1937        //
1938        // p is probably prime (strong pseudoprime) if for any a < p,
1939        // 1) a^t mod p = 1 or
1940        // 2) a^((2^j)*t) mod p = p-1 for some 0 <= j <= s-1
1941        //
1942        // Otherwise, p is composite.
1943        //
1944        // Returns
1945        // -------
1946        // True if "this" is a strong pseudoprime to randomly chosen
1947        // bases.  The number of chosen bases is given by the "confidence"
1948        // parameter.
1949        //
1950        // False if "this" is definitely NOT prime.
1951        //
1952        //***********************************************************************
1953
1954        public bool RabinMillerTest(int confidence)
1955        {
1956            BigInteger thisVal;
1957            if ((this.data[maxLength - 1] & 0x80000000) != 0)        // negative
1958                thisVal = -this;
1959            else
1960                thisVal = this;
1961
1962            if (thisVal.dataLength == 1)
1963            {
1964                // test small numbers
1965                if (thisVal.data[0] == 0 || thisVal.data[0] == 1)
1966                    return false;
1967                else if (thisVal.data[0] == 2 || thisVal.data[0] == 3)
1968                    return true;
1969            }
1970
1971            if ((thisVal.data[0] & 0x1) == 0)     // even numbers
1972                return false;
1973
1974
1975            // calculate values of s and t
1976            BigInteger p_sub1 = thisVal - (new BigInteger(1));
1977            int s = 0;
1978
1979            for (int index = 0; index < p_sub1.dataLength; index++)
1980            {
1981                uint mask = 0x01;
1982
1983                for (int i = 0; i < 32; i++)
1984                {
1985                    if ((p_sub1.data[index] & mask) != 0)
1986                    {
1987                        index = p_sub1.dataLength;      // to break the outer loop
1988                        break;
1989                    }
1990                    mask <<= 1;
1991                    s++;
1992                }
1993            }
1994
1995            BigInteger t = p_sub1 >> s;
1996
1997            int bits = thisVal.bitCount();
1998            BigInteger a = new BigInteger();
1999            Random rand = new Random();
2000
2001            for (int round = 0; round < confidence; round++)
2002            {
2003                bool done = false;
2004
2005                while (!done)           // generate a < n
2006                {
2007                    int testBits = 0;
2008
2009                    // make sure "a" has at least 2 bits
2010                    while (testBits < 2)
2011                        testBits = (int)(rand.NextDouble() * bits);
2012
2013                    a.genRandomBits(testBits, rand);
2014
2015                    int byteLen = a.dataLength;
2016
2017                    // make sure "a" is not 0
2018                    if (byteLen > 1 || (byteLen == 1 && a.data[0] != 1))
2019                        done = true;
2020                }
2021
2022                // check whether a factor exists (fix for version 1.03)
2023                BigInteger gcdTest = a.gcd(thisVal);
2024                if (gcdTest.dataLength == 1 && gcdTest.data[0] != 1)
2025                    return false;
2026
2027                BigInteger b = a.modPow(t, thisVal);
2028
2029                /*
2030                Console.WriteLine("a = " + a.ToString(10));
2031                Console.WriteLine("b = " + b.ToString(10));
2032                Console.WriteLine("t = " + t.ToString(10));
2033                Console.WriteLine("s = " + s);
2034                */
2035
2036                bool result = false;
2037
2038                if (b.dataLength == 1 && b.data[0] == 1)         // a^t mod p = 1
2039                    result = true;
2040
2041                for (int j = 0; result == false && j < s; j++)
2042                {
2043                    if (b == p_sub1)         // a^((2^j)*t) mod p = p-1 for some 0 <= j <= s-1
2044                    {
2045                        result = true;
2046                        break;
2047                    }
2048
2049                    b = (b * b) % thisVal;
2050                }
2051
2052                if (result == false)
2053                    return false;
2054            }
2055            return true;
2056        }
2057
2058
2059        //***********************************************************************
2060        // Probabilistic prime test based on Solovay-Strassen (Euler Criterion)
2061        //
2062        // p is probably prime if for any a < p (a is not multiple of p),
2063        // a^((p-1)/2) mod p = J(a, p)
2064        //
2065        // where J is the Jacobi symbol.
2066        //
2067        // Otherwise, p is composite.
2068        //
2069        // Returns
2070        // -------
2071        // True if "this" is a Euler pseudoprime to randomly chosen
2072        // bases.  The number of chosen bases is given by the "confidence"
2073        // parameter.
2074        //
2075        // False if "this" is definitely NOT prime.
2076        //
2077        //***********************************************************************
2078
2079        public bool SolovayStrassenTest(int confidence)
2080        {
2081            BigInteger thisVal;
2082            if ((this.data[maxLength - 1] & 0x80000000) != 0)        // negative
2083                thisVal = -this;
2084            else
2085                thisVal = this;
2086
2087            if (thisVal.dataLength == 1)
2088            {
2089                // test small numbers
2090                if (thisVal.data[0] == 0 || thisVal.data[0] == 1)
2091                    return false;
2092                else if (thisVal.data[0] == 2 || thisVal.data[0] == 3)
2093                    return true;
2094            }
2095
2096            if ((thisVal.data[0] & 0x1) == 0)     // even numbers
2097                return false;
2098
2099
2100            int bits = thisVal.bitCount();
2101            BigInteger a = new BigInteger();
2102            BigInteger p_sub1 = thisVal - 1;
2103            BigInteger p_sub1_shift = p_sub1 >> 1;
2104
2105            Random rand = new Random();
2106
2107            for (int round = 0; round < confidence; round++)
2108            {
2109                bool done = false;
2110
2111                while (!done)           // generate a < n
2112                {
2113                    int testBits = 0;
2114
2115                    // make sure "a" has at least 2 bits
2116                    while (testBits < 2)
2117                        testBits = (int)(rand.NextDouble() * bits);
2118
2119                    a.genRandomBits(testBits, rand);
2120
2121                    int byteLen = a.dataLength;
2122
2123                    // make sure "a" is not 0
2124                    if (byteLen > 1 || (byteLen == 1 && a.data[0] != 1))
2125                        done = true;
2126                }
2127
2128                // check whether a factor exists (fix for version 1.03)
2129                BigInteger gcdTest = a.gcd(thisVal);
2130                if (gcdTest.dataLength == 1 && gcdTest.data[0] != 1)
2131                    return false;
2132
2133                // calculate a^((p-1)/2) mod p
2134
2135                BigInteger expResult = a.modPow(p_sub1_shift, thisVal);
2136                if (expResult == p_sub1)
2137                    expResult = -1;
2138
2139                // calculate Jacobi symbol
2140                BigInteger jacob = Jacobi(a, thisVal);
2141
2142                //Console.WriteLine("a = " + a.ToString(10) + " b = " + thisVal.ToString(10));
2143                //Console.WriteLine("expResult = " + expResult.ToString(10) + " Jacob = " + jacob.ToString(10));
2144
2145                // if they are different then it is not prime
2146                if (expResult != jacob)
2147                    return false;
2148            }
2149
2150            return true;
2151        }
2152
2153
2154        //***********************************************************************
2155        // Implementation of the Lucas Strong Pseudo Prime test.
2156        //
2157        // Let n be an odd number with gcd(n,D) = 1, and n - J(D, n) = 2^s * d
2158        // with d odd and s >= 0.
2159        //
2160        // If Ud mod n = 0 or V2^r*d mod n = 0 for some 0 <= r < s, then n
2161        // is a strong Lucas pseudoprime with parameters (P, Q).  We select
2162        // P and Q based on Selfridge.
2163        //
2164        // Returns True if number is a strong Lucus pseudo prime.
2165        // Otherwise, returns False indicating that number is composite.
2166        //***********************************************************************
2167
2168        public bool LucasStrongTest()
2169        {
2170            BigInteger thisVal;
2171            if ((this.data[maxLength - 1] & 0x80000000) != 0)        // negative
2172                thisVal = -this;
2173            else
2174                thisVal = this;
2175
2176            if (thisVal.dataLength == 1)
2177            {
2178                // test small numbers
2179                if (thisVal.data[0] == 0 || thisVal.data[0] == 1)
2180                    return false;
2181                else if (thisVal.data[0] == 2 || thisVal.data[0] == 3)
2182                    return true;
2183            }
2184
2185            if ((thisVal.data[0] & 0x1) == 0)     // even numbers
2186                return false;
2187
2188            return LucasStrongTestHelper(thisVal);
2189        }
2190
2191
2192        private bool LucasStrongTestHelper(BigInteger thisVal)
2193        {
2194            // Do the test (selects D based on Selfridge)
2195            // Let D be the first element of the sequence
2196            // 5, -7, 9, -11, 13, ... for which J(D,n) = -1
2197            // Let P = 1, Q = (1-D) / 4
2198
2199            long D = 5, sign = -1, dCount = 0;
2200            bool done = false;
2201
2202            while (!done)
2203            {
2204                int Jresult = BigInteger.Jacobi(D, thisVal);
2205
2206                if (Jresult == -1)
2207                    done = true;    // J(D, this) = 1
2208                else
2209                {
2210                    if (Jresult == 0 && Math.Abs(D) < thisVal)       // divisor found
2211                        return false;
2212
2213                    if (dCount == 20)
2214                    {
2215                        // check for square
2216                        BigInteger root = thisVal.sqrt();
2217                        if (root * root == thisVal)
2218                            return false;
2219                    }
2220
2221                    //Console.WriteLine(D);
2222                    D = (Math.Abs(D) + 2) * sign;
2223                    sign = -sign;
2224                }
2225                dCount++;
2226            }
2227
2228            long Q = (1 - D) >> 2;
2229
2230            /*
2231            Console.WriteLine("D = " + D);
2232            Console.WriteLine("Q = " + Q);
2233            Console.WriteLine("(n,D) = " + thisVal.gcd(D));
2234            Console.WriteLine("(n,Q) = " + thisVal.gcd(Q));
2235            Console.WriteLine("J(D|n) = " + BigInteger.Jacobi(D, thisVal));
2236            */
2237
2238            BigInteger p_add1 = thisVal + 1;
2239            int s = 0;
2240
2241            for (int index = 0; index < p_add1.dataLength; index++)
2242            {
2243                uint mask = 0x01;
2244
2245                for (int i = 0; i < 32; i++)
2246                {
2247                    if ((p_add1.data[index] & mask) != 0)
2248                    {
2249                        index = p_add1.dataLength;      // to break the outer loop
2250                        break;
2251                    }
2252                    mask <<= 1;
2253                    s++;
2254                }
2255            }
2256
2257            BigInteger t = p_add1 >> s;
2258
2259            // calculate constant = b^(2k) / m
2260            // for Barrett Reduction
2261            BigInteger constant = new BigInteger();
2262
2263            int nLen = thisVal.dataLength << 1;
2264            constant.data[nLen] = 0x00000001;
2265            constant.dataLength = nLen + 1;
2266
2267            constant = constant / thisVal;
2268
2269            BigInteger[] lucas = LucasSequenceHelper(1, Q, t, thisVal, constant, 0);
2270            bool isPrime = false;
2271
2272            if ((lucas[0].dataLength == 1 && lucas[0].data[0] == 0) ||
2273               (lucas[1].dataLength == 1 && lucas[1].data[0] == 0))
2274            {
2275                // u(t) = 0 or V(t) = 0
2276                isPrime = true;
2277            }
2278
2279            for (int i = 1; i < s; i++)
2280            {
2281                if (!isPrime)
2282                {
2283                    // doubling of index
2284                    lucas[1] = thisVal.BarrettReduction(lucas[1] * lucas[1], thisVal, constant);
2285                    lucas[1] = (lucas[1] - (lucas[2] << 1)) % thisVal;
2286
2287                    //lucas[1] = ((lucas[1] * lucas[1]) - (lucas[2] << 1)) % thisVal;
2288
2289                    if ((lucas[1].dataLength == 1 && lucas[1].data[0] == 0))
2290                        isPrime = true;
2291                }
2292
2293                lucas[2] = thisVal.BarrettReduction(lucas[2] * lucas[2], thisVal, constant);     //Q^k
2294            }
2295
2296
2297            if (isPrime)     // additional checks for composite numbers
2298            {
2299                // If n is prime and gcd(n, Q) == 1, then
2300                // Q^((n+1)/2) = Q * Q^((n-1)/2) is congruent to (Q * J(Q, n)) mod n
2301
2302                BigInteger g = thisVal.gcd(Q);
2303                if (g.dataLength == 1 && g.data[0] == 1)         // gcd(this, Q) == 1
2304                {
2305                    if ((lucas[2].data[maxLength - 1] & 0x80000000) != 0)
2306                        lucas[2] += thisVal;
2307
2308                    BigInteger temp = (Q * BigInteger.Jacobi(Q, thisVal)) % thisVal;
2309                    if ((temp.data[maxLength - 1] & 0x80000000) != 0)
2310                        temp += thisVal;
2311
2312                    if (lucas[2] != temp)
2313                        isPrime = false;
2314                }
2315            }
2316
2317            return isPrime;
2318        }
2319
2320
2321        //***********************************************************************
2322        // Determines whether a number is probably prime, using the Rabin-Miller's
2323        // test.  Before applying the test, the number is tested for divisibility
2324        // by primes < 2000
2325        //
2326        // Returns true if number is probably prime.
2327        //***********************************************************************
2328
2329        public bool isProbablePrime(int confidence)
2330        {
2331            BigInteger thisVal;
2332            if ((this.data[maxLength - 1] & 0x80000000) != 0)        // negative
2333                thisVal = -this;
2334            else
2335                thisVal = this;
2336
2337
2338            // test for divisibility by primes < 2000
2339            for (int p = 0; p < primesBelow2000.Length; p++)
2340            {
2341                BigInteger divisor = primesBelow2000[p];
2342
2343                if (divisor >= thisVal)
2344                    break;
2345
2346                BigInteger resultNum = thisVal % divisor;
2347                if (resultNum.IntValue() == 0)
2348                {
2349                    /*
2350    Console.WriteLine("Not prime!  Divisible by {0}\n",
2351                                      primesBelow2000[p]);
2352                    */
2353                    return false;
2354                }
2355            }
2356
2357            if (thisVal.RabinMillerTest(confidence))
2358                return true;
2359            else
2360            {
2361                //Console.WriteLine("Not prime!  Failed primality test\n");
2362                return false;
2363            }
2364        }
2365
2366
2367        //***********************************************************************
2368        // Determines whether this BigInteger is probably prime using a
2369        // combination of base 2 strong pseudoprime test and Lucas strong
2370        // pseudoprime test.
2371        //
2372        // The sequence of the primality test is as follows,
2373        //
2374        // 1) Trial divisions are carried out using prime numbers below 2000.
2375        //    if any of the primes divides this BigInteger, then it is not prime.
2376        //
2377        // 2) Perform base 2 strong pseudoprime test.  If this BigInteger is a
2378        //    base 2 strong pseudoprime, proceed on to the next step.
2379        //
2380        // 3) Perform strong Lucas pseudoprime test.
2381        //
2382        // Returns True if this BigInteger is both a base 2 strong pseudoprime
2383        // and a strong Lucas pseudoprime.
2384        //
2385        // For a detailed discussion of this primality test, see [6].
2386        //
2387        //***********************************************************************
2388
2389        public bool isProbablePrime()
2390        {
2391            BigInteger thisVal;
2392            if ((this.data[maxLength - 1] & 0x80000000) != 0)        // negative
2393                thisVal = -this;
2394            else
2395                thisVal = this;
2396
2397            if (thisVal.dataLength == 1)
2398            {
2399                // test small numbers
2400                if (thisVal.data[0] == 0 || thisVal.data[0] == 1)
2401                    return false;
2402                else if (thisVal.data[0] == 2 || thisVal.data[0] == 3)
2403                    return true;
2404            }
2405
2406            if ((thisVal.data[0] & 0x1) == 0)     // even numbers
2407                return false;
2408
2409
2410            // test for divisibility by primes < 2000
2411            for (int p = 0; p < primesBelow2000.Length; p++)
2412            {
2413                BigInteger divisor = primesBelow2000[p];
2414
2415                if (divisor >= thisVal)
2416                    break;
2417
2418                BigInteger resultNum = thisVal % divisor;
2419                if (resultNum.IntValue() == 0)
2420                {
2421                    //Console.WriteLine("Not prime!  Divisible by {0}\n",
2422                    //                  primesBelow2000[p]);
2423
2424                    return false;
2425                }
2426            }
2427
2428            // Perform BASE 2 Rabin-Miller Test
2429
2430            // calculate values of s and t
2431            BigInteger p_sub1 = thisVal - (new BigInteger(1));
2432            int s = 0;
2433
2434            for (int index = 0; index < p_sub1.dataLength; index++)
2435            {
2436                uint mask = 0x01;
2437
2438                for (int i = 0; i < 32; i++)
2439                {
2440                    if ((p_sub1.data[index] & mask) != 0)
2441                    {
2442                        index = p_sub1.dataLength;      // to break the outer loop
2443                        break;
2444                    }
2445                    mask <<= 1;
2446                    s++;
2447                }
2448            }
2449
2450            BigInteger t = p_sub1 >> s;
2451
2452            int bits = thisVal.bitCount();
2453            BigInteger a = 2;
2454
2455            // b = a^t mod p
2456            BigInteger b = a.modPow(t, thisVal);
2457            bool result = false;
2458
2459            if (b.dataLength == 1 && b.data[0] == 1)         // a^t mod p = 1
2460                result = true;
2461
2462            for (int j = 0; result == false && j < s; j++)
2463            {
2464                if (b == p_sub1)         // a^((2^j)*t) mod p = p-1 for some 0 <= j <= s-1
2465                {
2466                    result = true;
2467                    break;
2468                }
2469
2470                b = (b * b) % thisVal;
2471            }
2472
2473            // if number is strong pseudoprime to base 2, then do a strong lucas test
2474            if (result)
2475                result = LucasStrongTestHelper(thisVal);
2476
2477            return result;
2478        }
2479
2480
2481
2482        //***********************************************************************
2483        // Returns the lowest 4 bytes of the BigInteger as an int.
2484        //***********************************************************************
2485
2486        public int IntValue()
2487        {
2488            return (int)data[0];
2489        }
2490
2491
2492        //***********************************************************************
2493        // Returns the lowest 8 bytes of the BigInteger as a long.
2494        //***********************************************************************
2495
2496        public long LongValue()
2497        {
2498            long val = 0;
2499
2500            val = (long)data[0];
2501            try
2502            {       // exception if maxLength = 1
2503                val |= (long)data[1] << 32;
2504            }
2505            catch (Exception)
2506            {
2507                if ((data[0] & 0x80000000) != 0) // negative
2508                    val = (int)data[0];
2509            }
2510
2511            return val;
2512        }
2513
2514
2515        //***********************************************************************
2516        // Computes the Jacobi Symbol for a and b.
2517        // Algorithm adapted from [3] and [4] with some optimizations
2518        //***********************************************************************
2519
2520        public static int Jacobi(BigInteger a, BigInteger b)
2521        {
2522            // Jacobi defined only for odd integers
2523            if ((b.data[0] & 0x1) == 0)
2524                throw (new ArgumentException("Jacobi defined only for odd integers."));
2525
2526            if (a >= b) a %= b;
2527            if (a.dataLength == 1 && a.data[0] == 0) return 0;  // a == 0
2528            if (a.dataLength == 1 && a.data[0] == 1) return 1;  // a == 1
2529
2530            if (a < 0)
2531            {
2532                if ((((b - 1).data[0]) & 0x2) == 0)       //if( (((b-1) >> 1).data[0] & 0x1) == 0)
2533                    return Jacobi(-a, b);
2534                else
2535                    return -Jacobi(-a, b);
2536            }
2537
2538            int e = 0;
2539            for (int index = 0; index < a.dataLength; index++)
2540            {
2541                uint mask = 0x01;
2542
2543                for (int i = 0; i < 32; i++)
2544                {
2545                    if ((a.data[index] & mask) != 0)
2546                    {
2547                        index = a.dataLength;      // to break the outer loop
2548                        break;
2549                    }
2550                    mask <<= 1;
2551                    e++;
2552                }
2553            }
2554
2555            BigInteger a1 = a >> e;
2556
2557            int s = 1;
2558            if ((e & 0x1) != 0 && ((b.data[0] & 0x7) == 3 || (b.data[0] & 0x7) == 5))
2559                s = -1;
2560
2561            if ((b.data[0] & 0x3) == 3 && (a1.data[0] & 0x3) == 3)
2562                s = -s;
2563
2564            if (a1.dataLength == 1 && a1.data[0] == 1)
2565                return s;
2566            else
2567                return (s * Jacobi(b % a1, a1));
2568        }
2569
2570
2571
2572        //***********************************************************************
2573        // Generates a positive BigInteger that is probably prime.
2574        //***********************************************************************
2575
2576        public static BigInteger genPseudoPrime(int bits, int confidence, Random rand)
2577        {
2578            BigInteger result = new BigInteger();
2579            bool done = false;
2580
2581            while (!done)
2582            {
2583                result.genRandomBits(bits, rand);
2584                result.data[0] |= 0x01;         // make it odd
2585
2586                // prime test
2587                done = result.isProbablePrime(confidence);
2588            }
2589            return result;
2590        }
2591
2592
2593        //***********************************************************************
2594        // Generates a random number with the specified number of bits such
2595        // that gcd(number, this) = 1
2596        //***********************************************************************
2597
2598        public BigInteger genCoPrime(int bits, Random rand)
2599        {
2600            bool done = false;
2601            BigInteger result = new BigInteger();
2602
2603            while (!done)
2604            {
2605                result.genRandomBits(bits, rand);
2606                //Console.WriteLine(result.ToString(16));
2607
2608                // gcd test
2609                BigInteger g = result.gcd(this);
2610                if (g.dataLength == 1 && g.data[0] == 1)
2611                    done = true;
2612            }
2613
2614            return result;
2615        }
2616
2617
2618        //***********************************************************************
2619        // Returns the modulo inverse of this.  Throws ArithmeticException if
2620        // the inverse does not exist.  (i.e. gcd(this, modulus) != 1)
2621        //***********************************************************************
2622
2623        public BigInteger modInverse(BigInteger modulus)
2624        {
2625            BigInteger[] p = { 0, 1 };
2626            BigInteger[] q = new BigInteger[2];    // quotients
2627            BigInteger[] r = { 0, 0 };             // remainders
2628
2629            int step = 0;
2630
2631            BigInteger a = modulus;
2632            BigInteger b = this;
2633
2634            while (b.dataLength > 1 || (b.dataLength == 1 && b.data[0] != 0))
2635            {
2636                BigInteger quotient = new BigInteger();
2637                BigInteger remainder = new BigInteger();
2638
2639                if (step > 1)
2640                {
2641                    BigInteger pval = (p[0] - (p[1] * q[0])) % modulus;
2642                    p[0] = p[1];
2643                    p[1] = pval;
2644                }
2645
2646                if (b.dataLength == 1)
2647                    singleByteDivide(a, b, quotient, remainder);
2648                else
2649                    multiByteDivide(a, b, quotient, remainder);
2650
2651                /*
2652                Console.WriteLine(quotient.dataLength);
2653                Console.WriteLine("{0} = {1}({2}) + {3}  p = {4}", a.ToString(10),
2654                                  b.ToString(10), quotient.ToString(10), remainder.ToString(10),
2655                                  p[1].ToString(10));
2656                */
2657
2658                q[0] = q[1];
2659                r[0] = r[1];
2660                q[1] = quotient; r[1] = remainder;
2661
2662                a = b;
2663                b = remainder;
2664
2665                step++;
2666            }
2667
2668            if (r[0].dataLength > 1 || (r[0].dataLength == 1 && r[0].data[0] != 1))
2669                throw (new ArithmeticException("No inverse!"));
2670
2671            BigInteger result = ((p[0] - (p[1] * q[0])) % modulus);
2672
2673            if ((result.data[maxLength - 1] & 0x80000000) != 0)
2674                result += modulus;  // get the least positive modulus
2675
2676            return result;
2677        }
2678
2679
2680        //***********************************************************************
2681        // Returns the value of the BigInteger as a byte array.  The lowest
2682        // index contains the MSB.
2683        //***********************************************************************
2684
2685        public byte[] getBytes()
2686        {
2687            int numBits = bitCount();
2688
2689            int numBytes = numBits >> 3;
2690            if ((numBits == 0) || ((numBits & 0x7) != 0)) // if numbits == 0, we still need a single byte to represent zero!
2691                numBytes++;
2692
2693            byte[] result = new byte[numBytes];
2694
2695            //Console.WriteLine(result.Length);
2696
2697            int pos = 0;
2698            uint tempVal, val = data[dataLength - 1];
2699
2700            bool jumped = false;
2701            if ((tempVal = (val >> 24 & 0xFF)) != 0)
2702            {
2703                result[pos++] = (byte)tempVal;
2704                jumped = true;
2705            }
2706            if ((tempVal = (val >> 16 & 0xFF)) != 0 || jumped)
2707            {
2708                result[pos++] = (byte)tempVal;
2709                jumped = true;
2710            }
2711            if ((tempVal = (val >> 8 & 0xFF)) != 0 || jumped)
2712            {
2713                result[pos++] = (byte)tempVal;
2714                jumped = true;
2715            }
2716            if ((tempVal = (val & 0xFF)) != 0 || jumped)
2717            {
2718                result[pos++] = (byte)tempVal;
2719                jumped = true;
2720            }
2721
2722            for (int i = dataLength - 2; i >= 0; i--, pos += 4)
2723            {
2724                val = data[i];
2725                result[pos + 3] = (byte)(val & 0xFF);
2726                val >>= 8;
2727                result[pos + 2] = (byte)(val & 0xFF);
2728                val >>= 8;
2729                result[pos + 1] = (byte)(val & 0xFF);
2730                val >>= 8;
2731                result[pos] = (byte)(val & 0xFF);
2732            }
2733
2734            return result;
2735        }
2736
2737
2738        //***********************************************************************
2739        // Sets the value of the specified bit to 1
2740        // The Least Significant Bit position is 0.
2741        //***********************************************************************
2742
2743        public void setBit(uint bitNum)
2744        {
2745            uint bytePos = bitNum >> 5;             // divide by 32
2746            byte bitPos = (byte)(bitNum & 0x1F);    // get the lowest 5 bits
2747
2748            uint mask = (uint)1 << bitPos;
2749            this.data[bytePos] |= mask;
2750
2751            if (bytePos >= this.dataLength)
2752                this.dataLength = (int)bytePos + 1;
2753        }
2754
2755
2756        //***********************************************************************
2757        // Sets the value of the specified bit to 0
2758        // The Least Significant Bit position is 0.
2759        //***********************************************************************
2760
2761        public void unsetBit(uint bitNum)
2762        {
2763            uint bytePos = bitNum >> 5;
2764
2765            if (bytePos < this.dataLength)
2766            {
2767                byte bitPos = (byte)(bitNum & 0x1F);
2768
2769                uint mask = (uint)1 << bitPos;
2770                uint mask2 = 0xFFFFFFFF ^ mask;
2771
2772                this.data[bytePos] &= mask2;
2773
2774                if (this.dataLength > 1 && this.data[this.dataLength - 1] == 0)
2775                    this.dataLength--;
2776            }
2777        }
2778
2779
2780        //***********************************************************************
2781        // Returns a value that is equivalent to the integer square root
2782        // of the BigInteger.
2783        //
2784        // The integer square root of "this" is defined as the largest integer n
2785        // such that (n * n) <= this
2786        //
2787        //***********************************************************************
2788
2789        public BigInteger sqrt()
2790        {
2791            uint numBits = (uint)this.bitCount();
2792
2793            if ((numBits & 0x1) != 0)        // odd number of bits
2794                numBits = (numBits >> 1) + 1;
2795            else
2796                numBits = (numBits >> 1);
2797
2798            uint bytePos = numBits >> 5;
2799            byte bitPos = (byte)(numBits & 0x1F);
2800
2801            uint mask;
2802
2803            BigInteger result = new BigInteger();
2804            if (bitPos == 0)
2805                mask = 0x80000000;
2806            else
2807            {
2808                mask = (uint)1 << bitPos;
2809                bytePos++;
2810            }
2811            result.dataLength = (int)bytePos;
2812
2813            for (int i = (int)bytePos - 1; i >= 0; i--)
2814            {
2815                while (mask != 0)
2816                {
2817                    // guess
2818                    result.data[i] ^= mask;
2819
2820                    // undo the guess if its square is larger than this
2821                    if ((result * result) > this)
2822                        result.data[i] ^= mask;
2823
2824                    mask >>= 1;
2825                }
2826                mask = 0x80000000;
2827            }
2828            return result;
2829        }
2830
2831
2832        //***********************************************************************
2833        // Returns the k_th number in the Lucas Sequence reduced modulo n.
2834        //
2835        // Uses index doubling to speed up the process.  For example, to calculate V(k),
2836        // we maintain two numbers in the sequence V(n) and V(n+1).
2837        //
2838        // To obtain V(2n), we use the identity
2839        //      V(2n) = (V(n) * V(n)) - (2 * Q^n)
2840        // To obtain V(2n+1), we first write it as
2841        //      V(2n+1) = V((n+1) + n)
2842        // and use the identity
2843        //      V(m+n) = V(m) * V(n) - Q * V(m-n)
2844        // Hence,
2845        //      V((n+1) + n) = V(n+1) * V(n) - Q^n * V((n+1) - n)
2846        //                   = V(n+1) * V(n) - Q^n * V(1)
2847        //                   = V(n+1) * V(n) - Q^n * P
2848        //
2849        // We use k in its binary expansion and perform index doubling for each
2850        // bit position.  For each bit position that is set, we perform an
2851        // index doubling followed by an index addition.  This means that for V(n),
2852        // we need to update it to V(2n+1).  For V(n+1), we need to update it to
2853        // V((2n+1)+1) = V(2*(n+1))
2854        //
2855        // This function returns
2856        // [0] = U(k)
2857        // [1] = V(k)
2858        // [2] = Q^n
2859        //
2860        // Where U(0) = 0 % n, U(1) = 1 % n
2861        //       V(0) = 2 % n, V(1) = P % n
2862        //***********************************************************************
2863
2864        public static BigInteger[] LucasSequence(BigInteger P, BigInteger Q,
2865                                                 BigInteger k, BigInteger n)
2866        {
2867            if (k.dataLength == 1 && k.data[0] == 0)
2868            {
2869                BigInteger[] result = new BigInteger[3];
2870
2871                result[0] = 0; result[1] = 2 % n; result[2] = 1 % n;
2872                return result;
2873            }
2874
2875            // calculate constant = b^(2k) / m
2876            // for Barrett Reduction
2877            BigInteger constant = new BigInteger();
2878
2879            int nLen = n.dataLength << 1;
2880            constant.data[nLen] = 0x00000001;
2881            constant.dataLength = nLen + 1;
2882
2883            constant = constant / n;
2884
2885            // calculate values of s and t
2886            int s = 0;
2887
2888            for (int index = 0; index < k.dataLength; index++)
2889            {
2890                uint mask = 0x01;
2891
2892                for (int i = 0; i < 32; i++)
2893                {
2894                    if ((k.data[index] & mask) != 0)
2895                    {
2896                        index = k.dataLength;      // to break the outer loop
2897                        break;
2898                    }
2899                    mask <<= 1;
2900                    s++;
2901                }
2902            }
2903
2904            BigInteger t = k >> s;
2905
2906            //Console.WriteLine("s = " + s + " t = " + t);
2907            return LucasSequenceHelper(P, Q, t, n, constant, s);
2908        }
2909
2910
2911        //***********************************************************************
2912        // Performs the calculation of the kth term in the Lucas Sequence.
2913        // For details of the algorithm, see reference [9].
2914        //
2915        // k must be odd.  i.e LSB == 1
2916        //***********************************************************************
2917
2918        private static BigInteger[] LucasSequenceHelper(BigInteger P, BigInteger Q,
2919                                                        BigInteger k, BigInteger n,
2920                                                        BigInteger constant, int s)
2921        {
2922            BigInteger[] result = new BigInteger[3];
2923
2924            if ((k.data[0] & 0x00000001) == 0)
2925                throw (new ArgumentException("Argument k must be odd."));
2926
2927            int numbits = k.bitCount();
2928            uint mask = (uint)0x1 << ((numbits & 0x1F) - 1);
2929
2930            // v = v0, v1 = v1, u1 = u1, Q_k = Q^0
2931
2932            BigInteger v = 2 % n, Q_k = 1 % n,
2933                       v1 = P % n, u1 = Q_k;
2934            bool flag = true;
2935
2936            for (int i = k.dataLength - 1; i >= 0; i--)     // iterate on the binary expansion of k
2937            {
2938                //Console.WriteLine("round");
2939                while (mask != 0)
2940                {
2941                    if (i == 0 && mask == 0x00000001)        // last bit
2942                        break;
2943
2944                    if ((k.data[i] & mask) != 0)             // bit is set
2945                    {
2946                        // index doubling with addition
2947
2948                        u1 = (u1 * v1) % n;
2949
2950                        v = ((v * v1) - (P * Q_k)) % n;
2951                        v1 = n.BarrettReduction(v1 * v1, n, constant);
2952                        v1 = (v1 - ((Q_k * Q) << 1)) % n;
2953
2954                        if (flag)
2955                            flag = false;
2956                        else
2957                            Q_k = n.BarrettReduction(Q_k * Q_k, n, constant);
2958
2959                        Q_k = (Q_k * Q) % n;
2960                    }
2961                    else
2962                    {
2963                        // index doubling
2964                        u1 = ((u1 * v) - Q_k) % n;
2965
2966                        v1 = ((v * v1) - (P * Q_k)) % n;
2967                        v = n.BarrettReduction(v * v, n, constant);
2968                        v = (v - (Q_k << 1)) % n;
2969
2970                        if (flag)
2971                        {
2972                            Q_k = Q % n;
2973                            flag = false;
2974                        }
2975                        else
2976                            Q_k = n.BarrettReduction(Q_k * Q_k, n, constant);
2977                    }
2978
2979                    mask >>= 1;
2980                }
2981                mask = 0x80000000;
2982            }
2983
2984            // at this point u1 = u(n+1) and v = v(n)
2985            // since the last bit always 1, we need to transform u1 to u(2n+1) and v to v(2n+1)
2986
2987            u1 = ((u1 * v) - Q_k) % n;
2988            v = ((v * v1) - (P * Q_k)) % n;
2989            if (flag)
2990                flag = false;
2991            else
2992                Q_k = n.BarrettReduction(Q_k * Q_k, n, constant);
2993
2994            Q_k = (Q_k * Q) % n;
2995
2996
2997            for (int i = 0; i < s; i++)
2998            {
2999                // index doubling
3000                u1 = (u1 * v) % n;
3001                v = ((v * v) - (Q_k << 1)) % n;
3002
3003                if (flag)
3004                {
3005                    Q_k = Q % n;
3006                    flag = false;
3007                }
3008                else
3009                    Q_k = n.BarrettReduction(Q_k * Q_k, n, constant);
3010            }
3011
3012            result[0] = u1;
3013            result[1] = v;
3014            result[2] = Q_k;
3015
3016            return result;
3017        }
3018
3019
3020        //***********************************************************************
3021        // Tests the correct implementation of the /, %, * and + operators
3022        //***********************************************************************
3023
3024        public static void MulDivTest(int rounds)
3025        {
3026            Random rand = new Random();
3027            byte[] val = new byte[64];
3028            byte[] val2 = new byte[64];
3029
3030            for (int count = 0; count < rounds; count++)
3031            {
3032                // generate 2 numbers of random length
3033                int t1 = 0;
3034                while (t1 == 0)
3035                    t1 = (int)(rand.NextDouble() * 65);
3036
3037                int t2 = 0;
3038                while (t2 == 0)
3039                    t2 = (int)(rand.NextDouble() * 65);
3040
3041                bool done = false;
3042                while (!done)
3043                {
3044                    for (int i = 0; i < 64; i++)
3045                    {
3046                        if (i < t1)
3047                            val[i] = (byte)(rand.NextDouble() * 256);
3048                        else
3049                            val[i] = 0;
3050
3051                        if (val[i] != 0)
3052                            done = true;
3053                    }
3054                }
3055
3056                done = false;
3057                while (!done)
3058                {
3059                    for (int i = 0; i < 64; i++)
3060                    {
3061                        if (i < t2)
3062                            val2[i] = (byte)(rand.NextDouble() * 256);
3063                        else
3064                            val2[i] = 0;
3065
3066                        if (val2[i] != 0)
3067                            done = true;
3068                    }
3069                }
3070
3071                while (val[0] == 0)
3072                    val[0] = (byte)(rand.NextDouble() * 256);
3073                while (val2[0] == 0)
3074                    val2[0] = (byte)(rand.NextDouble() * 256);
3075
3076                Console.WriteLine(count);
3077                BigInteger bn1 = new BigInteger(val, t1);
3078                BigInteger bn2 = new BigInteger(val2, t2);
3079
3080
3081                // Determine the quotient and remainder by dividing
3082                // the first number by the second.
3083
3084                BigInteger bn3 = bn1 / bn2;
3085                BigInteger bn4 = bn1 % bn2;
3086
3087                // Recalculate the number
3088                BigInteger bn5 = (bn3 * bn2) + bn4;
3089
3090                // Make sure they're the same
3091                if (bn5 != bn1)
3092                {
3093                    Console.WriteLine("Error at " + count);
3094                    Console.WriteLine(bn1 + "\n");
3095                    Console.WriteLine(bn2 + "\n");
3096                    Console.WriteLine(bn3 + "\n");
3097                    Console.WriteLine(bn4 + "\n");
3098                    Console.WriteLine(bn5 + "\n");
3099                    return;
3100                }
3101            }
3102        }
3103
3104
3105        //***********************************************************************
3106        // Tests the correct implementation of the modulo exponential function
3107        // using RSA encryption and decryption (using pre-computed encryption and
3108        // decryption keys).
3109        //***********************************************************************
3110
3111        public static void RSATest(int rounds)
3112        {
3113            Random rand = new Random(1);
3114            byte[] val = new byte[64];
3115
3116            // private and public key
3117            BigInteger bi_e = new BigInteger("a932b948feed4fb2b692609bd22164fc9edb59fae7880cc1eaff7b3c9626b7e5b241c27a974833b2622ebe09beb451917663d47232488f23a117fc97720f1e7", 16);
3118            BigInteger bi_d = new BigInteger("4adf2f7a89da93248509347d2ae506d683dd3a16357e859a980c4f77a4e2f7a01fae289f13a851df6e9db5adaa60bfd2b162bbbe31f7c8f828261a6839311929d2cef4f864dde65e556ce43c89bbbf9f1ac5511315847ce9cc8dc92470a747b8792d6a83b0092d2e5ebaf852c85cacf34278efa99160f2f8aa7ee7214de07b7", 16);
3119            BigInteger bi_n = new BigInteger("e8e77781f36a7b3188d711c2190b560f205a52391b3479cdb99fa010745cbeba5f2adc08e1de6bf38398a0487c4a73610d94ec36f17f3f46ad75e17bc1adfec99839589f45f95ccc94cb2a5c500b477eb3323d8cfab0c8458c96f0147a45d27e45a4d11d54d77684f65d48f15fafcc1ba208e71e921b9bd9017c16a5231af7f", 16);
3120
3121            Console.WriteLine("e =\n" + bi_e.ToString(10));
3122            Console.WriteLine("\nd =\n" + bi_d.ToString(10));
3123            Console.WriteLine("\nn =\n" + bi_n.ToString(10) + "\n");
3124
3125            for (int count = 0; count < rounds; count++)
3126            {
3127                // generate data of random length
3128                int t1 = 0;
3129                while (t1 == 0)
3130                    t1 = (int)(rand.NextDouble() * 65);
3131
3132                bool done = false;
3133                while (!done)
3134                {
3135                    for (int i = 0; i < 64; i++)
3136                    {
3137                        if (i < t1)
3138                            val[i] = (byte)(rand.NextDouble() * 256);
3139                        else
3140                            val[i] = 0;
3141
3142                        if (val[i] != 0)
3143                            done = true;
3144                    }
3145                }
3146
3147                while (val[0] == 0)
3148                    val[0] = (byte)(rand.NextDouble() * 256);
3149
3150                Console.Write("Round = " + count);
3151
3152                // encrypt and decrypt data
3153                BigInteger bi_data = new BigInteger(val, t1);
3154                BigInteger bi_encrypted = bi_data.modPow(bi_e, bi_n);
3155                BigInteger bi_decrypted = bi_encrypted.modPow(bi_d, bi_n);
3156
3157                // compare
3158                if (bi_decrypted != bi_data)
3159                {
3160                    Console.WriteLine("\nError at round " + count);
3161                    Console.WriteLine(bi_data + "\n");
3162                    return;
3163                }
3164                Console.WriteLine(" <PASSED>.");
3165            }
3166
3167        }
3168
3169
3170        //***********************************************************************
3171        // Tests the correct implementation of the modulo exponential and
3172        // inverse modulo functions using RSA encryption and decryption.  The two
3173        // pseudoprimes p and q are fixed, but the two RSA keys are generated
3174        // for each round of testing.
3175        //***********************************************************************
3176
3177        public static void RSATest2(int rounds)
3178        {
3179            Random rand = new Random();
3180            byte[] val = new byte[64];
3181
3182            byte[] pseudoPrime1 = {
3183                            (byte)0x85, (byte)0x84, (byte)0x64, (byte)0xFD, (byte)0x70, (byte)0x6A,
3184                            (byte)0x9F, (byte)0xF0, (byte)0x94, (byte)0x0C, (byte)0x3E, (byte)0x2C,
3185                            (byte)0x74, (byte)0x34, (byte)0x05, (byte)0xC9, (byte)0x55, (byte)0xB3,
3186                            (byte)0x85, (byte)0x32, (byte)0x98, (byte)0x71, (byte)0xF9, (byte)0x41,
3187                            (byte)0x21, (byte)0x5F, (byte)0x02, (byte)0x9E, (byte)0xEA, (byte)0x56,
3188                            (byte)0x8D, (byte)0x8C, (byte)0x44, (byte)0xCC, (byte)0xEE, (byte)0xEE,
3189                            (byte)0x3D, (byte)0x2C, (byte)0x9D, (byte)0x2C, (byte)0x12, (byte)0x41,
3190                            (byte)0x1E, (byte)0xF1, (byte)0xC5, (byte)0x32, (byte)0xC3, (byte)0xAA,
3191                            (byte)0x31, (byte)0x4A, (byte)0x52, (byte)0xD8, (byte)0xE8, (byte)0xAF,
3192                            (byte)0x42, (byte)0xF4, (byte)0x72, (byte)0xA1, (byte)0x2A, (byte)0x0D,
3193                            (byte)0x97, (byte)0xB1, (byte)0x31, (byte)0xB3,
3194                    };
3195
3196            byte[] pseudoPrime2 = {
3197                            (byte)0x99, (byte)0x98, (byte)0xCA, (byte)0xB8, (byte)0x5E, (byte)0xD7,
3198                            (byte)0xE5, (byte)0xDC, (byte)0x28, (byte)0x5C, (byte)0x6F, (byte)0x0E,
3199                            (byte)0x15, (byte)0x09, (byte)0x59, (byte)0x6E, (byte)0x84, (byte)0xF3,
3200                            (byte)0x81, (byte)0xCD, (byte)0xDE, (byte)0x42, (byte)0xDC, (byte)0x93,
3201                            (byte)0xC2, (byte)0x7A, (byte)0x62, (byte)0xAC, (byte)0x6C, (byte)0xAF,
3202                            (byte)0xDE, (byte)0x74, (byte)0xE3, (byte)0xCB, (byte)0x60, (byte)0x20,
3203                            (byte)0x38, (byte)0x9C, (byte)0x21, (byte)0xC3, (byte)0xDC, (byte)0xC8,
3204                            (byte)0xA2, (byte)0x4D, (byte)0xC6, (byte)0x2A, (byte)0x35, (byte)0x7F,
3205                            (byte)0xF3, (byte)0xA9, (byte)0xE8, (byte)0x1D, (byte)0x7B, (byte)0x2C,
3206                            (byte)0x78, (byte)0xFA, (byte)0xB8, (byte)0x02, (byte)0x55, (byte)0x80,
3207                            (byte)0x9B, (byte)0xC2, (byte)0xA5, (byte)0xCB,
3208                    };
3209
3210
3211            BigInteger bi_p = new BigInteger(pseudoPrime1);
3212            BigInteger bi_q = new BigInteger(pseudoPrime2);
3213            BigInteger bi_pq = (bi_p - 1) * (bi_q - 1);
3214            BigInteger bi_n = bi_p * bi_q;
3215
3216            for (int count = 0; count < rounds; count++)
3217            {
3218                // generate private and public key
3219                BigInteger bi_e = bi_pq.genCoPrime(512, rand);
3220                BigInteger bi_d = bi_e.modInverse(bi_pq);
3221
3222                Console.WriteLine("\ne =\n" + bi_e.ToString(10));
3223                Console.WriteLine("\nd =\n" + bi_d.ToString(10));
3224                Console.WriteLine("\nn =\n" + bi_n.ToString(10) + "\n");
3225
3226                // generate data of random length
3227                int t1 = 0;
3228                while (t1 == 0)
3229                    t1 = (int)(rand.NextDouble() * 65);
3230
3231                bool done = false;
3232                while (!done)
3233                {
3234                    for (int i = 0; i < 64; i++)
3235                    {
3236                        if (i < t1)
3237                            val[i] = (byte)(rand.NextDouble() * 256);
3238                        else
3239                            val[i] = 0;
3240
3241                        if (val[i] != 0)
3242                            done = true;
3243                    }
3244                }
3245
3246                while (val[0] == 0)
3247                    val[0] = (byte)(rand.NextDouble() * 256);
3248
3249                Console.Write("Round = " + count);
3250
3251                // encrypt and decrypt data
3252                BigInteger bi_data = new BigInteger(val, t1);
3253                BigInteger bi_encrypted = bi_data.modPow(bi_e, bi_n);
3254                BigInteger bi_decrypted = bi_encrypted.modPow(bi_d, bi_n);
3255
3256                // compare
3257                if (bi_decrypted != bi_data)
3258                {
3259                    Console.WriteLine("\nError at round " + count);
3260                    Console.WriteLine(bi_data + "\n");
3261                    return;
3262                }
3263                Console.WriteLine(" <PASSED>.");
3264            }
3265
3266        }
3267
3268
3269        //***********************************************************************
3270        // Tests the correct implementation of sqrt() method.
3271        //***********************************************************************
3272
3273        public static void SqrtTest(int rounds)
3274        {
3275            Random rand = new Random();
3276            for (int count = 0; count < rounds; count++)
3277            {
3278                // generate data of random length
3279                int t1 = 0;
3280                while (t1 == 0)
3281                    t1 = (int)(rand.NextDouble() * 1024);
3282
3283                Console.Write("Round = " + count);
3284
3285                BigInteger a = new BigInteger();
3286                a.genRandomBits(t1, rand);
3287
3288                BigInteger b = a.sqrt();
3289                BigInteger c = (b + 1) * (b + 1);
3290
3291                // check that b is the largest integer such that b*b <= a
3292                if (c <= a)
3293                {
3294                    Console.WriteLine("\nError at round " + count);
3295                    Console.WriteLine(a + "\n");
3296                    return;
3297                }
3298                Console.WriteLine(" <PASSED>.");
3299            }
3300        }
3301
3302
3303
3304        public static void Main(string[] args)
3305        {
3306            // Known problem -> these two pseudoprimes passes my implementation of
3307            // primality test but failed in JDK's isProbablePrime test.
3308
3309            byte[] pseudoPrime1 = { (byte)0x00,
3310                            (byte)0x85, (byte)0x84, (byte)0x64, (byte)0xFD, (byte)0x70, (byte)0x6A,
3311                            (byte)0x9F, (byte)0xF0, (byte)0x94, (byte)0x0C, (byte)0x3E, (byte)0x2C,
3312                            (byte)0x74, (byte)0x34, (byte)0x05, (byte)0xC9, (byte)0x55, (byte)0xB3,
3313                            (byte)0x85, (byte)0x32, (byte)0x98, (byte)0x71, (byte)0xF9, (byte)0x41,
3314                            (byte)0x21, (byte)0x5F, (byte)0x02, (byte)0x9E, (byte)0xEA, (byte)0x56,
3315                            (byte)0x8D, (byte)0x8C, (byte)0x44, (byte)0xCC, (byte)0xEE, (byte)0xEE,
3316                            (byte)0x3D, (byte)0x2C, (byte)0x9D, (byte)0x2C, (byte)0x12, (byte)0x41,
3317                            (byte)0x1E, (byte)0xF1, (byte)0xC5, (byte)0x32, (byte)0xC3, (byte)0xAA,
3318                            (byte)0x31, (byte)0x4A, (byte)0x52, (byte)0xD8, (byte)0xE8, (byte)0xAF,
3319                            (byte)0x42, (byte)0xF4, (byte)0x72, (byte)0xA1, (byte)0x2A, (byte)0x0D,
3320                            (byte)0x97, (byte)0xB1, (byte)0x31, (byte)0xB3,
3321                    };
3322
3323            byte[] pseudoPrime2 = { (byte)0x00,
3324                            (byte)0x99, (byte)0x98, (byte)0xCA, (byte)0xB8, (byte)0x5E, (byte)0xD7,
3325                            (byte)0xE5, (byte)0xDC, (byte)0x28, (byte)0x5C, (byte)0x6F, (byte)0x0E,
3326                            (byte)0x15, (byte)0x09, (byte)0x59, (byte)0x6E, (byte)0x84, (byte)0xF3,
3327                            (byte)0x81, (byte)0xCD, (byte)0xDE, (byte)0x42, (byte)0xDC, (byte)0x93,
3328                            (byte)0xC2, (byte)0x7A, (byte)0x62, (byte)0xAC, (byte)0x6C, (byte)0xAF,
3329                            (byte)0xDE, (byte)0x74, (byte)0xE3, (byte)0xCB, (byte)0x60, (byte)0x20,
3330                            (byte)0x38, (byte)0x9C, (byte)0x21, (byte)0xC3, (byte)0xDC, (byte)0xC8,
3331                            (byte)0xA2, (byte)0x4D, (byte)0xC6, (byte)0x2A, (byte)0x35, (byte)0x7F,
3332                            (byte)0xF3, (byte)0xA9, (byte)0xE8, (byte)0x1D, (byte)0x7B, (byte)0x2C,
3333                            (byte)0x78, (byte)0xFA, (byte)0xB8, (byte)0x02, (byte)0x55, (byte)0x80,
3334                            (byte)0x9B, (byte)0xC2, (byte)0xA5, (byte)0xCB,
3335                    };
3336
3337            Console.WriteLine("List of primes < 2000\n---------------------");
3338            int limit = 100, count = 0;
3339            for (int i = 0; i < 2000; i++)
3340            {
3341                if (i >= limit)
3342                {
3343                    Console.WriteLine();
3344                    limit += 100;
3345                }
3346
3347                BigInteger p = new BigInteger(-i);
3348
3349                if (p.isProbablePrime())
3350                {
3351                    Console.Write(i + ", ");
3352                    count++;
3353                }
3354            }
3355            Console.WriteLine("\nCount = " + count);
3356
3357
3358            BigInteger bi1 = new BigInteger(pseudoPrime1);
3359            Console.WriteLine("\n\nPrimality testing for\n" + bi1.ToString() + "\n");
3360            Console.WriteLine("SolovayStrassenTest(5) = " + bi1.SolovayStrassenTest(5));
3361            Console.WriteLine("RabinMillerTest(5) = " + bi1.RabinMillerTest(5));
3362            Console.WriteLine("FermatLittleTest(5) = " + bi1.FermatLittleTest(5));
3363            Console.WriteLine("isProbablePrime() = " + bi1.isProbablePrime());
3364
3365            Console.Write("\nGenerating 512-bits random pseudoprime. . .");
3366            Random rand = new Random();
3367            BigInteger prime = BigInteger.genPseudoPrime(512, 5, rand);
3368            Console.WriteLine("\n" + prime);
3369
3370            //int dwStart = System.Environment.TickCount;
3371            //BigInteger.MulDivTest(100000);
3372            //BigInteger.RSATest(10);
3373            //BigInteger.RSATest2(10);
3374            //Console.WriteLine(System.Environment.TickCount - dwStart);
3375
3376        }
3377       
3378        /**
3379         * Calculates the log bas of this BigInteger         
3380         */
3381        public double log(double bas) {
3382            int b = this.bitCount() - 1;
3383            double c = 0;
3384            double d = 1;
3385            for (int i = b; i >= 0; --i) {
3386                if ((this.data[i/32] & (1 << (i % 32))) != 0)
3387                    c += d;
3388                d *= 0.5;
3389            }
3390            return (Math.Log(c) + Math.Log(2) * b) / Math.Log(bas);
3391        }
3392
3393        /**         
3394         * Calculates the power with exponent exp of this BigInteger
3395         */
3396        public BigInteger pow(BigInteger exp)
3397        {
3398            if (exp >= 0)
3399            {
3400                BigInteger tmp = 1;
3401                for (BigInteger i = 0; i < exp; i++)
3402                    tmp *= this;
3403                return tmp;
3404            }
3405            else
3406                throw new Exception("Pow without mod not possible with negative exponent!");
3407        }
3408
3409        #region IComparable Members
3410
3411        /*
3412         * Implements the CompareTo method of IComparable
3413         */
3414        int IComparable.CompareTo(object obj)
3415        {
3416            BigInteger other = null;
3417            if (obj is BigInteger)
3418                other = (BigInteger)obj;
3419            else if (obj is String)
3420                other = new BigInteger(obj as String, 10);
3421            else if (obj is Int32)
3422                other = new BigInteger((int)obj);
3423            else
3424                throw new Exception("Invalid comparison with BigInteger");
3425
3426            if (this == other)
3427                return 0;
3428            else if (this < other)
3429                return -1;
3430            else
3431                return 1;
3432        }
3433
3434        #endregion
3435
3436        #region internal stuff of expression parser
3437
3438        private struct TOKEN
3439        {
3440            public enum Ttype { MULTIPLY, DIVIDE, PLUS, MINUS, POW, BRACKETOPEN, BRACKETCLOSE, INTEGER };
3441            public Ttype ttype;
3442            public BigInteger integer;
3443        }
3444
3445        private static Stack<TOKEN> scan(string expr)
3446        {
3447            TOKEN t = new TOKEN();
3448            int startIndex = 0;
3449            if (expr == "")
3450                return new Stack<TOKEN>();
3451            switch (expr[0])
3452            {
3453                case ' ':
3454                    return scan(expr.Substring(1));
3455                case '(':                   
3456                    t.ttype = TOKEN.Ttype.BRACKETOPEN;
3457                    startIndex = 1;
3458                    break;
3459                case ')':                   
3460                    t.ttype = TOKEN.Ttype.BRACKETCLOSE;
3461                    startIndex = 1;
3462                    break;
3463                case '+':
3464                    t.ttype = TOKEN.Ttype.PLUS;
3465                    startIndex = 1;
3466                    break;
3467                case '-':
3468                    t.ttype = TOKEN.Ttype.MINUS;
3469                    startIndex = 1;
3470                    break;
3471                case '*':
3472                    t.ttype = TOKEN.Ttype.MULTIPLY;
3473                    startIndex = 1;
3474                    break;
3475                case '/':
3476                    t.ttype = TOKEN.Ttype.DIVIDE;
3477                    startIndex = 1;
3478                    break;
3479                case '^':
3480                    t.ttype = TOKEN.Ttype.POW;
3481                    startIndex = 1;
3482                    break;
3483                case '0':
3484                case '1':
3485                case '2':
3486                case '3':
3487                case '4':
3488                case '5':
3489                case '6':
3490                case '7':
3491                case '8':
3492                case '9':
3493                    int length = 1;
3494                    for (; length < expr.Length; length++)
3495                        if (!(expr[length] >= '0' && expr[length] <= '9'))
3496                            break;
3497                    t.integer = new BigInteger(expr.Substring(0, length), 10);
3498                    t.ttype = TOKEN.Ttype.INTEGER;
3499                    startIndex = length;
3500                    break;
3501                default:
3502                    throw new Exception("Expression parsing failed at character " + expr[0]);                   
3503            }
3504            Stack<TOKEN> st = scan(expr.Substring(startIndex));
3505            st.Push(t);
3506            return st;           
3507        }
3508
3509        private enum Priority { ALL, POW, MULTDIV, ADDSUB };
3510
3511        private static BigInteger parse(Stack<TOKEN> stack, Priority priority, bool endbracket)
3512        {
3513            if (stack.Count == 0)
3514                throw new Exception("Expression Parsing Error.");
3515            int minus = 1;
3516            BigInteger v = 0;
3517            TOKEN t = stack.Pop();  //get -, +, integer or bracket
3518
3519            if (t.ttype == TOKEN.Ttype.MINUS)
3520            {
3521                minus = -1;
3522                t = stack.Pop();    //get integer or bracket
3523            }
3524            else if (t.ttype == TOKEN.Ttype.PLUS)
3525            {
3526                minus = 1;
3527                t = stack.Pop();    //get integer or bracket
3528            }
3529
3530            if (t.ttype == TOKEN.Ttype.INTEGER)
3531            {
3532                v = minus * t.integer;
3533            }
3534            else if (t.ttype == TOKEN.Ttype.BRACKETOPEN)
3535            {
3536                v = minus * parse(stack, Priority.ALL, true);
3537                stack.Pop();    //pop the closing bracket
3538            }
3539
3540            while (stack.Count != 0)
3541            {       
3542                switch (stack.Peek().ttype) //next operator
3543                {
3544                    case TOKEN.Ttype.PLUS:
3545                        if (priority == Priority.MULTDIV || priority == Priority.POW)
3546                            return v;
3547                        stack.Pop();
3548                        v = v + parse(stack, Priority.ADDSUB, endbracket);
3549                        break;
3550                    case TOKEN.Ttype.MINUS:
3551                        if (priority == Priority.MULTDIV || priority == Priority.POW)
3552                            return v;
3553                        stack.Pop();
3554                        v = v - parse(stack, Priority.ADDSUB, endbracket);
3555                        break;
3556                    case TOKEN.Ttype.MULTIPLY:
3557                        if (priority == Priority.POW)
3558                            return v;
3559                        stack.Pop();
3560                        v = v * parse(stack, Priority.MULTDIV, endbracket);
3561                        break;
3562                    case TOKEN.Ttype.DIVIDE:
3563                        if (priority == Priority.POW)
3564                            return v;
3565                        stack.Pop();
3566                        v = v / parse(stack, Priority.MULTDIV, endbracket);
3567                        break;
3568                    case TOKEN.Ttype.POW:
3569                        stack.Pop();
3570                        v = v.pow(parse(stack, Priority.POW, endbracket));
3571                        break;
3572                    case TOKEN.Ttype.BRACKETCLOSE:
3573                        if (endbracket)
3574                            return v;
3575                        else
3576                            throw new Exception("Expression Parsing Error (closing bracket misplaced).");
3577                    default:
3578                        throw new Exception("Expression Parsing Error.");
3579                }
3580            }
3581            if (endbracket)
3582                throw new Exception("Expression Parsing Error (closing bracket missing).");
3583
3584            return v;
3585        }
3586
3587        #endregion
3588
3589        /*         
3590         * Parses a math expression (example: (2+2)^(17-5) )
3591         * and returns a BigInteger based on this expression
3592         *
3593         * throws an exception when expression is not valid or the Number gets too big
3594         */
3595        public static BigInteger parseExpression(string expr) {
3596            Stack<TOKEN> stack = scan(expr);
3597            BigInteger i = parse(stack, Priority.ALL, false);
3598            return i;
3599        }
3600    }
3601}
Note: See TracBrowser for help on using the repository browser.