diff options
Diffstat (limited to 'src/rt/bigint/bigint_ext.cpp')
| -rw-r--r-- | src/rt/bigint/bigint_ext.cpp | 553 |
1 files changed, 553 insertions, 0 deletions
diff --git a/src/rt/bigint/bigint_ext.cpp b/src/rt/bigint/bigint_ext.cpp new file mode 100644 index 00000000000..66d79106f48 --- /dev/null +++ b/src/rt/bigint/bigint_ext.cpp @@ -0,0 +1,553 @@ +/* bigint_ext - external portion of large integer package +** +** Copyright © 2000 by Jef Poskanzer <jef@mail.acme.com>. +** All rights reserved. +** +** Redistribution and use in source and binary forms, with or without +** modification, are permitted provided that the following conditions +** are met: +** 1. Redistributions of source code must retain the above copyright +** notice, this list of conditions and the following disclaimer. +** 2. Redistributions in binary form must reproduce the above copyright +** notice, this list of conditions and the following disclaimer in the +** documentation and/or other materials provided with the distribution. +** +** THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND +** ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +** IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +** ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE +** FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +** DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS +** OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +** HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY +** OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF +** SUCH DAMAGE. +*/ + +#include <sys/types.h> +#include <signal.h> +#include <stdio.h> +#include <stdlib.h> +#include <unistd.h> +#include <time.h> + +#include "bigint.h" +#include "low_primes.h" + + +bigint bi_0, bi_1, bi_2, bi_10, bi_m1, bi_maxint, bi_minint; + + +/* Forwards. */ +static void print_pos( FILE* f, bigint bi ); + + +bigint +str_to_bi( char* str ) + { + int sign; + bigint biR; + + sign = 1; + if ( *str == '-' ) + { + sign = -1; + ++str; + } + for ( biR = bi_0; *str >= '0' && *str <= '9'; ++str ) + biR = bi_int_add( bi_int_multiply( biR, 10 ), *str - '0' ); + if ( sign == -1 ) + biR = bi_negate( biR ); + return biR; + } + + +void +bi_print( FILE* f, bigint bi ) + { + if ( bi_is_negative( bi_copy( bi ) ) ) + { + putc( '-', f ); + bi = bi_negate( bi ); + } + print_pos( f, bi ); + } + + +bigint +bi_scan( FILE* f ) + { + int sign; + int c; + bigint biR; + + sign = 1; + c = getc( f ); + if ( c == '-' ) + sign = -1; + else + ungetc( c, f ); + + biR = bi_0; + for (;;) + { + c = getc( f ); + if ( c < '0' || c > '9' ) + break; + biR = bi_int_add( bi_int_multiply( biR, 10 ), c - '0' ); + } + + if ( sign == -1 ) + biR = bi_negate( biR ); + return biR; + } + + +static void +print_pos( FILE* f, bigint bi ) + { + if ( bi_compare( bi_copy( bi ), bi_10 ) >= 0 ) + print_pos( f, bi_int_divide( bi_copy( bi ), 10 ) ); + putc( bi_int_mod( bi, 10 ) + '0', f ); + } + + +int +bi_int_mod( bigint bi, int m ) + { + int r; + + if ( m <= 0 ) + { + (void) fprintf( stderr, "bi_int_mod: zero or negative modulus\n" ); + (void) kill( getpid(), SIGFPE ); + } + r = bi_int_rem( bi, m ); + if ( r < 0 ) + r += m; + return r; + } + + +bigint +bi_rem( bigint bia, bigint bim ) + { + return bi_subtract( + bia, bi_multiply( bi_divide( bi_copy( bia ), bi_copy( bim ) ), bim ) ); + } + + +bigint +bi_mod( bigint bia, bigint bim ) + { + bigint biR; + + if ( bi_compare( bi_copy( bim ), bi_0 ) <= 0 ) + { + (void) fprintf( stderr, "bi_mod: zero or negative modulus\n" ); + (void) kill( getpid(), SIGFPE ); + } + biR = bi_rem( bia, bi_copy( bim ) ); + if ( bi_is_negative( bi_copy( biR ) ) ) + biR = bi_add( biR, bim ); + else + bi_free( bim ); + return biR; + } + + +bigint +bi_square( bigint bi ) + { + bigint biR; + + biR = bi_multiply( bi_copy( bi ), bi_copy( bi ) ); + bi_free( bi ); + return biR; + } + + +bigint +bi_power( bigint bi, bigint biexp ) + { + bigint biR; + + if ( bi_is_negative( bi_copy( biexp ) ) ) + { + (void) fprintf( stderr, "bi_power: negative exponent\n" ); + (void) kill( getpid(), SIGFPE ); + } + biR = bi_1; + for (;;) + { + if ( bi_is_odd( bi_copy( biexp ) ) ) + biR = bi_multiply( biR, bi_copy( bi ) ); + biexp = bi_half( biexp ); + if ( bi_compare( bi_copy( biexp ), bi_0 ) <= 0 ) + break; + bi = bi_multiply( bi_copy( bi ), bi ); + } + bi_free( bi ); + bi_free( biexp ); + return biR; + } + + +bigint +bi_factorial( bigint bi ) + { + bigint biR; + + biR = bi_1; + while ( bi_compare( bi_copy( bi ), bi_1 ) > 0 ) + { + biR = bi_multiply( biR, bi_copy( bi ) ); + bi = bi_int_subtract( bi, 1 ); + } + bi_free( bi ); + return biR; + } + + +int +bi_is_even( bigint bi ) + { + return ! bi_is_odd( bi ); + } + + +bigint +bi_mod_power( bigint bi, bigint biexp, bigint bim ) + { + int invert; + bigint biR; + + invert = 0; + if ( bi_is_negative( bi_copy( biexp ) ) ) + { + biexp = bi_negate( biexp ); + invert = 1; + } + + biR = bi_1; + for (;;) + { + if ( bi_is_odd( bi_copy( biexp ) ) ) + biR = bi_mod( bi_multiply( biR, bi_copy( bi ) ), bi_copy( bim ) ); + biexp = bi_half( biexp ); + if ( bi_compare( bi_copy( biexp ), bi_0 ) <= 0 ) + break; + bi = bi_mod( bi_multiply( bi_copy( bi ), bi ), bi_copy( bim ) ); + } + bi_free( bi ); + bi_free( biexp ); + + if ( invert ) + biR = bi_mod_inverse( biR, bim ); + else + bi_free( bim ); + return biR; + } + + +bigint +bi_mod_inverse( bigint bi, bigint bim ) + { + bigint gcd, mul0, mul1; + + gcd = bi_egcd( bi_copy( bim ), bi, &mul0, &mul1 ); + + /* Did we get gcd == 1? */ + if ( ! bi_is_one( gcd ) ) + { + (void) fprintf( stderr, "bi_mod_inverse: not relatively prime\n" ); + (void) kill( getpid(), SIGFPE ); + } + + bi_free( mul0 ); + return bi_mod( mul1, bim ); + } + + +/* Euclid's algorithm. */ +bigint +bi_gcd( bigint bim, bigint bin ) + { + bigint bit; + + bim = bi_abs( bim ); + bin = bi_abs( bin ); + while ( ! bi_is_zero( bi_copy( bin ) ) ) + { + bit = bi_mod( bim, bi_copy( bin ) ); + bim = bin; + bin = bit; + } + bi_free( bin ); + return bim; + } + + +/* Extended Euclidean algorithm. */ +bigint +bi_egcd( bigint bim, bigint bin, bigint* bim_mul, bigint* bin_mul ) + { + bigint a0, b0, c0, a1, b1, c1, q, t; + + if ( bi_is_negative( bi_copy( bim ) ) ) + { + bigint biR; + + biR = bi_egcd( bi_negate( bim ), bin, &t, bin_mul ); + *bim_mul = bi_negate( t ); + return biR; + } + if ( bi_is_negative( bi_copy( bin ) ) ) + { + bigint biR; + + biR = bi_egcd( bim, bi_negate( bin ), bim_mul, &t ); + *bin_mul = bi_negate( t ); + return biR; + } + + a0 = bi_1; b0 = bi_0; c0 = bim; + a1 = bi_0; b1 = bi_1; c1 = bin; + + while ( ! bi_is_zero( bi_copy( c1 ) ) ) + { + q = bi_divide( bi_copy( c0 ), bi_copy( c1 ) ); + t = a0; + a0 = bi_copy( a1 ); + a1 = bi_subtract( t, bi_multiply( bi_copy( q ), a1 ) ); + t = b0; + b0 = bi_copy( b1 ); + b1 = bi_subtract( t, bi_multiply( bi_copy( q ), b1 ) ); + t = c0; + c0 = bi_copy( c1 ); + c1 = bi_subtract( t, bi_multiply( bi_copy( q ), c1 ) ); + bi_free( q ); + } + + bi_free( a1 ); + bi_free( b1 ); + bi_free( c1 ); + *bim_mul = a0; + *bin_mul = b0; + return c0; + } + + +bigint +bi_lcm( bigint bia, bigint bib ) + { + bigint biR; + + biR = bi_divide( + bi_multiply( bi_copy( bia ), bi_copy( bib ) ), + bi_gcd( bi_copy( bia ), bi_copy( bib ) ) ); + bi_free( bia ); + bi_free( bib ); + return biR; + } + + +/* The Jacobi symbol. */ +bigint +bi_jacobi( bigint bia, bigint bib ) + { + bigint biR; + + if ( bi_is_even( bi_copy( bib ) ) ) + { + (void) fprintf( stderr, "bi_jacobi: don't know how to compute Jacobi(n, even)\n" ); + (void) kill( getpid(), SIGFPE ); + } + + if ( bi_compare( bi_copy( bia ), bi_copy( bib ) ) >= 0 ) + return bi_jacobi( bi_mod( bia, bi_copy( bib ) ), bib ); + + if ( bi_is_zero( bi_copy( bia ) ) || bi_is_one( bi_copy( bia ) ) ) + { + bi_free( bib ); + return bia; + } + + if ( bi_compare( bi_copy( bia ), bi_2 ) == 0 ) + { + bi_free( bia ); + switch ( bi_int_mod( bib, 8 ) ) + { + case 1: case 7: + return bi_1; + case 3: case 5: + return bi_m1; + } + } + + if ( bi_is_even( bi_copy( bia ) ) ) + { + biR = bi_multiply( + bi_jacobi( bi_2, bi_copy( bib ) ), + bi_jacobi( bi_half( bia ), bi_copy( bib ) ) ); + bi_free( bib ); + return biR; + } + + if ( bi_int_mod( bi_copy( bia ), 4 ) == 3 && + bi_int_mod( bi_copy( bib ), 4 ) == 3 ) + return bi_negate( bi_jacobi( bib, bia ) ); + else + return bi_jacobi( bib, bia ); + } + + +/* Probabalistic prime checking. */ +int +bi_is_probable_prime( bigint bi, int certainty ) + { + int i, p; + bigint bim1; + + /* First do trial division by a list of small primes. This eliminates + ** many candidates. + */ + for ( i = 0; i < sizeof(low_primes)/sizeof(*low_primes); ++i ) + { + p = low_primes[i]; + switch ( bi_compare( int_to_bi( p ), bi_copy( bi ) ) ) + { + case 0: + bi_free( bi ); + return 1; + case 1: + bi_free( bi ); + return 0; + } + if ( bi_int_mod( bi_copy( bi ), p ) == 0 ) + { + bi_free( bi ); + return 0; + } + } + + /* Now do the probabilistic tests. */ + bim1 = bi_int_subtract( bi_copy( bi ), 1 ); + for ( i = 0; i < certainty; ++i ) + { + bigint a, j, jac; + + /* Pick random test number. */ + a = bi_random( bi_copy( bi ) ); + + /* Decide whether to run the Fermat test or the Solovay-Strassen + ** test. The Fermat test is fast but lets some composite numbers + ** through. Solovay-Strassen runs slower but is more certain. + ** So the compromise here is we run the Fermat test a couple of + ** times to quickly reject most composite numbers, and then do + ** the rest of the iterations with Solovay-Strassen so nothing + ** slips through. + */ + if ( i < 2 && certainty >= 5 ) + { + /* Fermat test. Note that this is not state of the art. There's a + ** class of numbers called Carmichael numbers which are composite + ** but look prime to this test - it lets them slip through no + ** matter how many reps you run. However, it's nice and fast so + ** we run it anyway to help quickly reject most of the composites. + */ + if ( ! bi_is_one( bi_mod_power( bi_copy( a ), bi_copy( bim1 ), bi_copy( bi ) ) ) ) + { + bi_free( bi ); + bi_free( bim1 ); + bi_free( a ); + return 0; + } + } + else + { + /* GCD test. This rarely hits, but we need it for Solovay-Strassen. */ + if ( ! bi_is_one( bi_gcd( bi_copy( bi ), bi_copy( a ) ) ) ) + { + bi_free( bi ); + bi_free( bim1 ); + bi_free( a ); + return 0; + } + + /* Solovay-Strassen test. First compute pseudo Jacobi. */ + j = bi_mod_power( + bi_copy( a ), bi_half( bi_copy( bim1 ) ), bi_copy( bi ) ); + if ( bi_compare( bi_copy( j ), bi_copy( bim1 ) ) == 0 ) + { + bi_free( j ); + j = bi_m1; + } + + /* Now compute real Jacobi. */ + jac = bi_jacobi( bi_copy( a ), bi_copy( bi ) ); + + /* If they're not equal, the number is definitely composite. */ + if ( bi_compare( j, jac ) != 0 ) + { + bi_free( bi ); + bi_free( bim1 ); + bi_free( a ); + return 0; + } + } + + bi_free( a ); + } + + bi_free( bim1 ); + + bi_free( bi ); + return 1; + } + + +bigint +bi_generate_prime( int bits, int certainty ) + { + bigint bimo2, bip; + int i, inc = 0; + + bimo2 = bi_power( bi_2, int_to_bi( bits - 1 ) ); + for (;;) + { + bip = bi_add( bi_random( bi_copy( bimo2 ) ), bi_copy( bimo2 ) ); + /* By shoving the candidate numbers up to the next highest multiple + ** of six plus or minus one, we pre-eliminate all multiples of + ** two and/or three. + */ + switch ( bi_int_mod( bi_copy( bip ), 6 ) ) + { + case 0: inc = 4; bip = bi_int_add( bip, 1 ); break; + case 1: inc = 4; break; + case 2: inc = 2; bip = bi_int_add( bip, 3 ); break; + case 3: inc = 2; bip = bi_int_add( bip, 2 ); break; + case 4: inc = 2; bip = bi_int_add( bip, 1 ); break; + case 5: inc = 2; break; + } + /* Starting from the generated random number, check a bunch of + ** numbers in sequence. This is just to avoid calls to bi_random(), + ** which is more expensive than a simple add. + */ + for ( i = 0; i < 1000; ++i ) /* arbitrary */ + { + if ( bi_is_probable_prime( bi_copy( bip ), certainty ) ) + { + bi_free( bimo2 ); + return bip; + } + bip = bi_int_add( bip, inc ); + inc = 6 - inc; + } + /* We ran through the whole sequence and didn't find a prime. + ** Shrug, just try a different random starting point. + */ + bi_free( bip ); + } + } |