/src/rt/bigint/bigint_ext.cpp
C++ | 553 lines | 403 code | 84 blank | 66 comment | 62 complexity | 69ef421e2b423e04b66db38b55b2c8ac MD5 | raw file
1/* bigint_ext - external portion of large integer package 2** 3** Copyright � 2000 by Jef Poskanzer <jef@mail.acme.com>. 4** All rights reserved. 5** 6** Redistribution and use in source and binary forms, with or without 7** modification, are permitted provided that the following conditions 8** are met: 9** 1. Redistributions of source code must retain the above copyright 10** notice, this list of conditions and the following disclaimer. 11** 2. Redistributions in binary form must reproduce the above copyright 12** notice, this list of conditions and the following disclaimer in the 13** documentation and/or other materials provided with the distribution. 14** 15** THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 16** ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 17** IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 18** ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 19** FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 20** DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 21** OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 22** HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 23** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 24** OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 25** SUCH DAMAGE. 26*/ 27 28#include <sys/types.h> 29#include <signal.h> 30#include <stdio.h> 31#include <stdlib.h> 32#include <unistd.h> 33#include <time.h> 34 35#include "bigint.h" 36#include "low_primes.h" 37 38 39bigint bi_0, bi_1, bi_2, bi_10, bi_m1, bi_maxint, bi_minint; 40 41 42/* Forwards. */ 43static void print_pos( FILE* f, bigint bi ); 44 45 46bigint 47str_to_bi( char* str ) 48 { 49 int sign; 50 bigint biR; 51 52 sign = 1; 53 if ( *str == '-' ) 54 { 55 sign = -1; 56 ++str; 57 } 58 for ( biR = bi_0; *str >= '0' && *str <= '9'; ++str ) 59 biR = bi_int_add( bi_int_multiply( biR, 10 ), *str - '0' ); 60 if ( sign == -1 ) 61 biR = bi_negate( biR ); 62 return biR; 63 } 64 65 66void 67bi_print( FILE* f, bigint bi ) 68 { 69 if ( bi_is_negative( bi_copy( bi ) ) ) 70 { 71 putc( '-', f ); 72 bi = bi_negate( bi ); 73 } 74 print_pos( f, bi ); 75 } 76 77 78bigint 79bi_scan( FILE* f ) 80 { 81 int sign; 82 int c; 83 bigint biR; 84 85 sign = 1; 86 c = getc( f ); 87 if ( c == '-' ) 88 sign = -1; 89 else 90 ungetc( c, f ); 91 92 biR = bi_0; 93 for (;;) 94 { 95 c = getc( f ); 96 if ( c < '0' || c > '9' ) 97 break; 98 biR = bi_int_add( bi_int_multiply( biR, 10 ), c - '0' ); 99 } 100 101 if ( sign == -1 ) 102 biR = bi_negate( biR ); 103 return biR; 104 } 105 106 107static void 108print_pos( FILE* f, bigint bi ) 109 { 110 if ( bi_compare( bi_copy( bi ), bi_10 ) >= 0 ) 111 print_pos( f, bi_int_divide( bi_copy( bi ), 10 ) ); 112 putc( bi_int_mod( bi, 10 ) + '0', f ); 113 } 114 115 116int 117bi_int_mod( bigint bi, int m ) 118 { 119 int r; 120 121 if ( m <= 0 ) 122 { 123 (void) fprintf( stderr, "bi_int_mod: zero or negative modulus\n" ); 124 (void) kill( getpid(), SIGFPE ); 125 } 126 r = bi_int_rem( bi, m ); 127 if ( r < 0 ) 128 r += m; 129 return r; 130 } 131 132 133bigint 134bi_rem( bigint bia, bigint bim ) 135 { 136 return bi_subtract( 137 bia, bi_multiply( bi_divide( bi_copy( bia ), bi_copy( bim ) ), bim ) ); 138 } 139 140 141bigint 142bi_mod( bigint bia, bigint bim ) 143 { 144 bigint biR; 145 146 if ( bi_compare( bi_copy( bim ), bi_0 ) <= 0 ) 147 { 148 (void) fprintf( stderr, "bi_mod: zero or negative modulus\n" ); 149 (void) kill( getpid(), SIGFPE ); 150 } 151 biR = bi_rem( bia, bi_copy( bim ) ); 152 if ( bi_is_negative( bi_copy( biR ) ) ) 153 biR = bi_add( biR, bim ); 154 else 155 bi_free( bim ); 156 return biR; 157 } 158 159 160bigint 161bi_square( bigint bi ) 162 { 163 bigint biR; 164 165 biR = bi_multiply( bi_copy( bi ), bi_copy( bi ) ); 166 bi_free( bi ); 167 return biR; 168 } 169 170 171bigint 172bi_power( bigint bi, bigint biexp ) 173 { 174 bigint biR; 175 176 if ( bi_is_negative( bi_copy( biexp ) ) ) 177 { 178 (void) fprintf( stderr, "bi_power: negative exponent\n" ); 179 (void) kill( getpid(), SIGFPE ); 180 } 181 biR = bi_1; 182 for (;;) 183 { 184 if ( bi_is_odd( bi_copy( biexp ) ) ) 185 biR = bi_multiply( biR, bi_copy( bi ) ); 186 biexp = bi_half( biexp ); 187 if ( bi_compare( bi_copy( biexp ), bi_0 ) <= 0 ) 188 break; 189 bi = bi_multiply( bi_copy( bi ), bi ); 190 } 191 bi_free( bi ); 192 bi_free( biexp ); 193 return biR; 194 } 195 196 197bigint 198bi_factorial( bigint bi ) 199 { 200 bigint biR; 201 202 biR = bi_1; 203 while ( bi_compare( bi_copy( bi ), bi_1 ) > 0 ) 204 { 205 biR = bi_multiply( biR, bi_copy( bi ) ); 206 bi = bi_int_subtract( bi, 1 ); 207 } 208 bi_free( bi ); 209 return biR; 210 } 211 212 213int 214bi_is_even( bigint bi ) 215 { 216 return ! bi_is_odd( bi ); 217 } 218 219 220bigint 221bi_mod_power( bigint bi, bigint biexp, bigint bim ) 222 { 223 int invert; 224 bigint biR; 225 226 invert = 0; 227 if ( bi_is_negative( bi_copy( biexp ) ) ) 228 { 229 biexp = bi_negate( biexp ); 230 invert = 1; 231 } 232 233 biR = bi_1; 234 for (;;) 235 { 236 if ( bi_is_odd( bi_copy( biexp ) ) ) 237 biR = bi_mod( bi_multiply( biR, bi_copy( bi ) ), bi_copy( bim ) ); 238 biexp = bi_half( biexp ); 239 if ( bi_compare( bi_copy( biexp ), bi_0 ) <= 0 ) 240 break; 241 bi = bi_mod( bi_multiply( bi_copy( bi ), bi ), bi_copy( bim ) ); 242 } 243 bi_free( bi ); 244 bi_free( biexp ); 245 246 if ( invert ) 247 biR = bi_mod_inverse( biR, bim ); 248 else 249 bi_free( bim ); 250 return biR; 251 } 252 253 254bigint 255bi_mod_inverse( bigint bi, bigint bim ) 256 { 257 bigint gcd, mul0, mul1; 258 259 gcd = bi_egcd( bi_copy( bim ), bi, &mul0, &mul1 ); 260 261 /* Did we get gcd == 1? */ 262 if ( ! bi_is_one( gcd ) ) 263 { 264 (void) fprintf( stderr, "bi_mod_inverse: not relatively prime\n" ); 265 (void) kill( getpid(), SIGFPE ); 266 } 267 268 bi_free( mul0 ); 269 return bi_mod( mul1, bim ); 270 } 271 272 273/* Euclid's algorithm. */ 274bigint 275bi_gcd( bigint bim, bigint bin ) 276 { 277 bigint bit; 278 279 bim = bi_abs( bim ); 280 bin = bi_abs( bin ); 281 while ( ! bi_is_zero( bi_copy( bin ) ) ) 282 { 283 bit = bi_mod( bim, bi_copy( bin ) ); 284 bim = bin; 285 bin = bit; 286 } 287 bi_free( bin ); 288 return bim; 289 } 290 291 292/* Extended Euclidean algorithm. */ 293bigint 294bi_egcd( bigint bim, bigint bin, bigint* bim_mul, bigint* bin_mul ) 295 { 296 bigint a0, b0, c0, a1, b1, c1, q, t; 297 298 if ( bi_is_negative( bi_copy( bim ) ) ) 299 { 300 bigint biR; 301 302 biR = bi_egcd( bi_negate( bim ), bin, &t, bin_mul ); 303 *bim_mul = bi_negate( t ); 304 return biR; 305 } 306 if ( bi_is_negative( bi_copy( bin ) ) ) 307 { 308 bigint biR; 309 310 biR = bi_egcd( bim, bi_negate( bin ), bim_mul, &t ); 311 *bin_mul = bi_negate( t ); 312 return biR; 313 } 314 315 a0 = bi_1; b0 = bi_0; c0 = bim; 316 a1 = bi_0; b1 = bi_1; c1 = bin; 317 318 while ( ! bi_is_zero( bi_copy( c1 ) ) ) 319 { 320 q = bi_divide( bi_copy( c0 ), bi_copy( c1 ) ); 321 t = a0; 322 a0 = bi_copy( a1 ); 323 a1 = bi_subtract( t, bi_multiply( bi_copy( q ), a1 ) ); 324 t = b0; 325 b0 = bi_copy( b1 ); 326 b1 = bi_subtract( t, bi_multiply( bi_copy( q ), b1 ) ); 327 t = c0; 328 c0 = bi_copy( c1 ); 329 c1 = bi_subtract( t, bi_multiply( bi_copy( q ), c1 ) ); 330 bi_free( q ); 331 } 332 333 bi_free( a1 ); 334 bi_free( b1 ); 335 bi_free( c1 ); 336 *bim_mul = a0; 337 *bin_mul = b0; 338 return c0; 339 } 340 341 342bigint 343bi_lcm( bigint bia, bigint bib ) 344 { 345 bigint biR; 346 347 biR = bi_divide( 348 bi_multiply( bi_copy( bia ), bi_copy( bib ) ), 349 bi_gcd( bi_copy( bia ), bi_copy( bib ) ) ); 350 bi_free( bia ); 351 bi_free( bib ); 352 return biR; 353 } 354 355 356/* The Jacobi symbol. */ 357bigint 358bi_jacobi( bigint bia, bigint bib ) 359 { 360 bigint biR; 361 362 if ( bi_is_even( bi_copy( bib ) ) ) 363 { 364 (void) fprintf( stderr, "bi_jacobi: don't know how to compute Jacobi(n, even)\n" ); 365 (void) kill( getpid(), SIGFPE ); 366 } 367 368 if ( bi_compare( bi_copy( bia ), bi_copy( bib ) ) >= 0 ) 369 return bi_jacobi( bi_mod( bia, bi_copy( bib ) ), bib ); 370 371 if ( bi_is_zero( bi_copy( bia ) ) || bi_is_one( bi_copy( bia ) ) ) 372 { 373 bi_free( bib ); 374 return bia; 375 } 376 377 if ( bi_compare( bi_copy( bia ), bi_2 ) == 0 ) 378 { 379 bi_free( bia ); 380 switch ( bi_int_mod( bib, 8 ) ) 381 { 382 case 1: case 7: 383 return bi_1; 384 case 3: case 5: 385 return bi_m1; 386 } 387 } 388 389 if ( bi_is_even( bi_copy( bia ) ) ) 390 { 391 biR = bi_multiply( 392 bi_jacobi( bi_2, bi_copy( bib ) ), 393 bi_jacobi( bi_half( bia ), bi_copy( bib ) ) ); 394 bi_free( bib ); 395 return biR; 396 } 397 398 if ( bi_int_mod( bi_copy( bia ), 4 ) == 3 && 399 bi_int_mod( bi_copy( bib ), 4 ) == 3 ) 400 return bi_negate( bi_jacobi( bib, bia ) ); 401 else 402 return bi_jacobi( bib, bia ); 403 } 404 405 406/* Probabalistic prime checking. */ 407int 408bi_is_probable_prime( bigint bi, int certainty ) 409 { 410 int i, p; 411 bigint bim1; 412 413 /* First do trial division by a list of small primes. This eliminates 414 ** many candidates. 415 */ 416 for ( i = 0; i < sizeof(low_primes)/sizeof(*low_primes); ++i ) 417 { 418 p = low_primes[i]; 419 switch ( bi_compare( int_to_bi( p ), bi_copy( bi ) ) ) 420 { 421 case 0: 422 bi_free( bi ); 423 return 1; 424 case 1: 425 bi_free( bi ); 426 return 0; 427 } 428 if ( bi_int_mod( bi_copy( bi ), p ) == 0 ) 429 { 430 bi_free( bi ); 431 return 0; 432 } 433 } 434 435 /* Now do the probabilistic tests. */ 436 bim1 = bi_int_subtract( bi_copy( bi ), 1 ); 437 for ( i = 0; i < certainty; ++i ) 438 { 439 bigint a, j, jac; 440 441 /* Pick random test number. */ 442 a = bi_random( bi_copy( bi ) ); 443 444 /* Decide whether to run the Fermat test or the Solovay-Strassen 445 ** test. The Fermat test is fast but lets some composite numbers 446 ** through. Solovay-Strassen runs slower but is more certain. 447 ** So the compromise here is we run the Fermat test a couple of 448 ** times to quickly reject most composite numbers, and then do 449 ** the rest of the iterations with Solovay-Strassen so nothing 450 ** slips through. 451 */ 452 if ( i < 2 && certainty >= 5 ) 453 { 454 /* Fermat test. Note that this is not state of the art. There's a 455 ** class of numbers called Carmichael numbers which are composite 456 ** but look prime to this test - it lets them slip through no 457 ** matter how many reps you run. However, it's nice and fast so 458 ** we run it anyway to help quickly reject most of the composites. 459 */ 460 if ( ! bi_is_one( bi_mod_power( bi_copy( a ), bi_copy( bim1 ), bi_copy( bi ) ) ) ) 461 { 462 bi_free( bi ); 463 bi_free( bim1 ); 464 bi_free( a ); 465 return 0; 466 } 467 } 468 else 469 { 470 /* GCD test. This rarely hits, but we need it for Solovay-Strassen. */ 471 if ( ! bi_is_one( bi_gcd( bi_copy( bi ), bi_copy( a ) ) ) ) 472 { 473 bi_free( bi ); 474 bi_free( bim1 ); 475 bi_free( a ); 476 return 0; 477 } 478 479 /* Solovay-Strassen test. First compute pseudo Jacobi. */ 480 j = bi_mod_power( 481 bi_copy( a ), bi_half( bi_copy( bim1 ) ), bi_copy( bi ) ); 482 if ( bi_compare( bi_copy( j ), bi_copy( bim1 ) ) == 0 ) 483 { 484 bi_free( j ); 485 j = bi_m1; 486 } 487 488 /* Now compute real Jacobi. */ 489 jac = bi_jacobi( bi_copy( a ), bi_copy( bi ) ); 490 491 /* If they're not equal, the number is definitely composite. */ 492 if ( bi_compare( j, jac ) != 0 ) 493 { 494 bi_free( bi ); 495 bi_free( bim1 ); 496 bi_free( a ); 497 return 0; 498 } 499 } 500 501 bi_free( a ); 502 } 503 504 bi_free( bim1 ); 505 506 bi_free( bi ); 507 return 1; 508 } 509 510 511bigint 512bi_generate_prime( int bits, int certainty ) 513 { 514 bigint bimo2, bip; 515 int i, inc = 0; 516 517 bimo2 = bi_power( bi_2, int_to_bi( bits - 1 ) ); 518 for (;;) 519 { 520 bip = bi_add( bi_random( bi_copy( bimo2 ) ), bi_copy( bimo2 ) ); 521 /* By shoving the candidate numbers up to the next highest multiple 522 ** of six plus or minus one, we pre-eliminate all multiples of 523 ** two and/or three. 524 */ 525 switch ( bi_int_mod( bi_copy( bip ), 6 ) ) 526 { 527 case 0: inc = 4; bip = bi_int_add( bip, 1 ); break; 528 case 1: inc = 4; break; 529 case 2: inc = 2; bip = bi_int_add( bip, 3 ); break; 530 case 3: inc = 2; bip = bi_int_add( bip, 2 ); break; 531 case 4: inc = 2; bip = bi_int_add( bip, 1 ); break; 532 case 5: inc = 2; break; 533 } 534 /* Starting from the generated random number, check a bunch of 535 ** numbers in sequence. This is just to avoid calls to bi_random(), 536 ** which is more expensive than a simple add. 537 */ 538 for ( i = 0; i < 1000; ++i ) /* arbitrary */ 539 { 540 if ( bi_is_probable_prime( bi_copy( bip ), certainty ) ) 541 { 542 bi_free( bimo2 ); 543 return bip; 544 } 545 bip = bi_int_add( bip, inc ); 546 inc = 6 - inc; 547 } 548 /* We ran through the whole sequence and didn't find a prime. 549 ** Shrug, just try a different random starting point. 550 */ 551 bi_free( bip ); 552 } 553 }