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