554 lines
12 KiB
C++
554 lines
12 KiB
C++
|
/* bigint_ext - external portion of large integer package
|
|||
|
|
**
|
|||
|
|
** Copyright <EFBFBD> 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 );
|
|||
|
|
}
|
|||
|
|
}
|