[3897] in cryptography@c2.net mail archive
Re: Large Primes
daemon@ATHENA.MIT.EDU (William Allen Simpson)
Sun Jan 3 20:30:32 1999
Date: Sun, 3 Jan 99 19:48:11 GMT
From: "William Allen Simpson" <wsimpson@greendragon.com>
To: Andrew Maslar <amaslar@home.com>
Cc: Cryptography List <cryptography@c2.net>
> Does anybody know of a source for large (1000+ bit) STRONG prime
> numbers?
> Generating software or archives would be helpful.
>
There are probably quite a few out there, but here is what Karn and I
worked up over the years to generate safe moduli for Photuris -- a pair
of programs: qsieve.c is memory intensive, qsafe.c is cpu intensive. It
uses an intermediate and final output file format that was discussed on
the photuris mailing list.
This generates "safe" primes, rather than "strong" primes.
I meant to put this in a RFC someday, but it would be nice to know
whether I'd done something wrong first.... Any problems/suggestions?
qsieve.c:
/* Sieve candidates for "safe" primes,
* suitable for use as Diffie-Hellman moduli;
* that is, where q = (p-1)/2 is also prime.
*
* This is the first of two steps. This step is memory intensive.
*
* 1996 May William Allen Simpson
* extracted from earlier code by Phil Karn, April 1994.
* save large primes list for later processing.
* 1998 May William Allen Simpson
* parameterized.
*/
#include <stdio.h>
#include <time.h>
#include <gmp.h> /* after stdio */
//#define SMALL_CHECK 1
//#define LARGE_CHECK 1
/* Chosen for commonly available machines with at least 12 MB of
* available main memory. Using virtual memory can cause thrashing.
* Therefore, this should be the largest number that is supported
* without a large amount of disk activity -- that would increase the
* run time from minutes to days or weeks!
*** Do not increase this number beyond the unsigned integer bit size.
*** Due to a multiple of 4, it must be LESS than 128 (yielding 2**30 bits).
*/
#define LARGE_MEMORY (12UL) /* megabytes */
/* Constant: When used with 32-bit integers, the largest sieve prime
* has to be less than 2**32.
*/
#define SMALL_MAXIMUM (0xffffffffUL)
/* Constant: can sieve all primes less than 2**32, as 65537**2 > 2**32-1.
*/
#define TINY_NUMBER (1UL<<16)
/* Ensure enough bit space for testing 2*q.
*/
#define TEST_DOUBLER (3) /* 2**n, n < 5 */
#define TEST_LIMIT (LARGE_MEMORY << ((18 - 6) - TEST_DOUBLER))
#define TEST_MAXIMUM (TEST_LIMIT < (1L<<16) ? TEST_LIMIT : (1L<<16))
#define TEST_MINIMUM (1L << 6)
/* Bit operations on 32-bit words
*/
#define BIT_CLEAR(a,n) ((a)[(n)>>5] &= ~(1L << ((n) & 31)))
#define BIT_SET(a,n) ((a)[(n)>>5] |= (1L << ((n) & 31)))
#define BIT_TEST(a,n) ((a)[(n)>>5] & (1L << ((n) & 31)))
typedef unsigned long int uint32; /* 32-bit unsigned integer */
uint32 *LargeSieve;
uint32 largewords;
uint32 largetries;
uint32 largenumbers;
uint32 largebits;
MP_INT largebase;
uint32 *SmallSieve;
uint32 smallbits;
uint32 smallbase;
uint32 *TinySieve;
uint32 tinybits;
/* Sieve out p's and q's with small factors
*/
void
sieve_large( uint32 s )
{
uint32 r;
uint32 u;
uint32 v;
#ifdef SMALL_CHECK
printf("%lu\n",s);
#endif
largetries++;
/* r = largebase mod s */
r = mpz_fdiv_ui(&largebase,s);
if(r == 0)
{
/* s divides into largebase exactly */
u = 0;
}
else
{
/* largebase+u is first entry divisible by s */
u = s - r;
}
if ( u < largebits * 2 )
{
/* The sieve omits p's and q's divisible by 2,
* so ensure that largebase+u is odd.
* Then, step through the sieve in increments of 2*s
*/
if (u & 0x1)
{
/* Make largebase+u odd, and u even */
u += s;
}
/* Mark all multiples of 2*s */
for ( u /= 2; u < largebits; u += s )
{
BIT_SET(LargeSieve,u);
};
}
/* r = p mod s */
r = (2*r+1) % s;
if(r == 0)
{
/* s divides p exactly */
u = 0;
}
else
{
/* p+u is first entry divisible by s */
u = s - r;
}
if ( u < largebits * 4 )
{
/* The sieve omits p's divisible by 4,
* so ensure that largebase+u is not.
* Then, step through the sieve in increments of 4*s
*/
while (u & 0x3)
{
if ( SMALL_MAXIMUM - u < s )
{
return;
}
u += s;
};
/* Mark all multiples of 4*s */
for ( u /= 4; u < largebits; u += s )
{
BIT_SET(LargeSieve,u);
};
}
}
main( int argc, char *argv[] )
{
MP_INT one;
MP_INT q;
struct tm *gtm;
FILE *largefile = NULL;
uint32 j;
uint32 power;
uint32 r;
uint32 s;
uint32 smallwords = TINY_NUMBER >> 6;
uint32 t;
time_t time_start;
time_t time_stop;
uint32 tinywords = TINY_NUMBER >> 6;
int i;
if ( argc < 3 )
{
printf("Missing argument: "
"<outfile> <bits> [<start>]\n");
exit(1);
}
/* Set power to the length in bits of the prime to be generated.
* This is generally 1 less than the desired safe prime moduli p.
*/
power = strtoul(argv[2],NULL,10);
if ( power < TEST_MINIMUM )
{
printf("Too few bits: %lu < %lu\n",
power,
TEST_MINIMUM );
exit(1);
}
else if ( power > TEST_MAXIMUM )
{
printf("Too many bits: %lu > %lu\n",
power,
TEST_MAXIMUM );
exit(1);
}
largefile = fopen(argv[1],"a+");
if ( largefile == NULL )
{
printf("Unable to write %s\n", argv[1] );
exit(1);
}
if ( (TinySieve = (uint32*)calloc(tinywords,sizeof(uint32))) == NULL )
{
printf("Insufficient memory for tiny sieve: need %lu bytes\n",
tinywords << 2 );
exit(1);
}
tinybits = tinywords << 5; /* 32-bit words */
if ( (SmallSieve = (uint32*)calloc(smallwords,sizeof(uint32))) == NULL )
{
printf("Insufficient memory for small sieve: need %lu bytes\n",
smallwords << 2 );
free(TinySieve);
exit(1);
}
smallbits = smallwords << 5; /* 32-bit words */
/* The density of ordinary primes is on the order of 1/bits,
* so the density of safe primes should be about (1/bits)**2.
* Set test range to something well above bits**2 to be reasonably
* sure (but not guaranteed) of catching at least one safe prime.
*/
largewords = ((power*power) >> (5-TEST_DOUBLER));
if ( largewords > (LARGE_MEMORY << 18) )
{
largewords = (LARGE_MEMORY << 18);
}
/* dynamically determine available memory */
while ( (LargeSieve = (uint32*)calloc(largewords,sizeof(uint32))) == NULL )
{
largewords -= (1L << 16); /* 1/4 MB chunks */
}
largebits = largewords << 5; /* 32-bit words */
largenumbers = largebits << 1; /* even numbers excluded */
/* validation check: count the number of primes tried */
largetries = 0;
mpz_init_set_ui(&one,1);
mpz_init(&q);
mpz_init(&largebase);
/* Generate random starting point for subprime search,
* or use specified parameter.
*/
if ( argc < 4 )
{
srandom(time_start);
mpz_random(&largebase,(power+32)/32);
/* Trim the random bits to the desired range
*/
mpz_mul_2exp(&q,&one,power);
mpz_sub_ui(&q,&q,1);
mpz_and(&largebase,&largebase,&q);
/* Set the high order bit to ensure the length
*/
mpz_mul_2exp(&q,&one,power-1);
mpz_ior(&largebase,&largebase,&q);
}
else
{
mpz_set_str(&largebase,argv[3],0);
}
/* ensure odd */
if((mpz_get_ui(&largebase) & 1) == 0)
{
mpz_sub_ui(&largebase,&largebase,1);
}
time(&time_start);
printf("%.24s Sieve large primes up to %lu from %lu-bit start point:\n",
ctime(&time_start),
largenumbers,
power );
mpz_out_str(NULL,16,&largebase);
printf("\n");
/* TinySieve */
for ( i = 0; i < tinybits; i++ )
{
if(BIT_TEST(TinySieve,i))
{
/* 2*i+3 is composite */
continue;
}
/* The next tiny prime */
t = 2*i + 3;
/* Mark all multiples of t */
for( j = i+t; j < tinybits; j += t )
{
BIT_SET(TinySieve,j);
};
sieve_large(t);
};
/* start the small block search at the next possible prime.
* to avoid fencepost errors, the last pass is skipped.
*/
for ( smallbase = TINY_NUMBER + 3;
smallbase < (SMALL_MAXIMUM - TINY_NUMBER);
smallbase += TINY_NUMBER )
{
for ( i = 0; i < tinybits; i++ )
{
if(BIT_TEST(TinySieve,i))
{
/* 2*i+3 is composite */
continue;
}
/* The next tiny prime */
t = 2*i + 3;
r = smallbase % t;
if (r == 0)
{
/* t divides into smallbase exactly */
s = 0;
}
else
{
/* smallbase+s is first entry divisible by t */
s = t - r;
}
/* The sieve omits even numbers,
* so ensure that smallbase+s is odd.
* Then, step through the sieve in increments of 2*t
*/
if (s & 1)
{
/* Make smallbase+s odd, and s even */
s += t;
}
/* Mark all multiples of 2*t */
for ( s /= 2; s < smallbits; s += t )
{
BIT_SET(SmallSieve,s);
};
};
/* SmallSieve */
for ( i = 0; i < smallbits; i++ )
{
if(BIT_TEST(SmallSieve,i))
{
/* 2*i+smallbase is composite */
continue;
}
/* The next small prime */
sieve_large( (2 * i) + smallbase );
};
memset( SmallSieve, 0, smallwords << 2 );
};
time(&time_stop);
printf("%.24s Sieved with %lu small primes in %lu seconds\n",
ctime(&time_stop),
largetries,
(long)(time_stop - time_start) );
for(j=r=0;j<largebits;j++)
{
if(BIT_TEST(LargeSieve,j))
{
/* Definitely composite, skip */
continue;
}
#ifdef LARGE_CHECK
printf("test q = largebase+%lu\n",2*j);
#endif
mpz_add_ui(&q,&largebase,2*j);
time(&time_stop);
gtm = gmtime(&time_stop);
fprintf(largefile,"%04d%02d%02d%02d%02d%02d ",
gtm->tm_year + 1900,
gtm->tm_mon + 1,
gtm->tm_mday,
gtm->tm_hour,
gtm->tm_min,
gtm->tm_sec );
fprintf(largefile,"4 "); /* Sophie-Germaine */
fprintf(largefile,"2 "); /* Sieve */
fprintf(largefile,"%lu ",largetries); /* iterations */
fprintf(largefile,"%lu ",power); /* bits */
fprintf(largefile,"0 "); /* generator unknown */
mpz_out_str(largefile,16,&q);
fprintf(largefile,"\n");
r++; /* count q */
}
time(&time_stop);
free(LargeSieve);
free(SmallSieve);
free(TinySieve);
fclose(largefile);
printf("%.24s Found %lu candidates\n",
ctime(&time_stop),
r );
}
qsafe.c:
/* Test probable "safe" primes,
* suitable for use as Diffie-Hellman moduli;
* that is, where q = (p-1)/2 is also prime.
*
* This is the second of two steps. This step is processor intensive.
*
* 1996 May William Allen Simpson
* extracted from earlier code by Phil Karn, April 1994.
* read large prime candidates list (q),
* and check prime probability of (p).
* 1998 May William Allen Simpson
* parameterized.
* optionally limit to a single generator.
*/
#include <stdio.h>
#include <time.h>
#include <gmp.h> /* after stdio */
//#define FINALCHECK 1
/* need line long enough for largest prime and headers */
#define LINESIZE 32768
#define TRIAL_MINIMUM (4)
#define TEST_MILLER_RABIN (0x04)
typedef unsigned long int uint32; /* 32-bit unsigned integer */
main( int argc, char *argv[] )
{
MP_INT q;
MP_INT p;
struct tm *gtm;
FILE *finalfile = NULL;
FILE *largefile = NULL;
char *cp;
char *lp;
uint32 count_candidates = 0;
uint32 count_safe = 0;
uint32 generator_known;
uint32 generator_wanted = 0;
uint32 in_tests;
uint32 in_tries;
uint32 in_type;
uint32 in_size;
uint32 trials;
time_t time_start;
time_t time_stop;
if ( argc < 4 )
{
printf("Missing argument: "
"<infile> <trials> <outfile> [<generator>]\n");
exit(1);
}
largefile = fopen(argv[1],"r");
if ( largefile == NULL )
{
printf("Unable to read %s\n", argv[1] );
exit(1);
}
if ( (trials = strtoul(argv[2],NULL,10)) < TRIAL_MINIMUM )
{
trials = TRIAL_MINIMUM;
}
finalfile = fopen(argv[3],"a+");
if ( finalfile == NULL )
{
printf("Unable to write %s\n", argv[3] );
exit(1);
}
if ( argc > 4 )
{
generator_wanted = strtoul(argv[4],NULL,16);
}
time(&time_start);
mpz_init(&p);
mpz_init(&q);
printf("%.24s Final %lu Miller-Rabin trials (%lx generator)\n",
ctime(&time_start),
trials,
generator_wanted );
lp = (char*)malloc(LINESIZE+1);
while ( fgets(lp,LINESIZE,largefile) != NULL )
{
int ll = strlen(lp);
if ( ll < 14 || *lp == '!' || *lp == '#' )
{
//print("comment or blank\n");
continue;
}
// time
cp = &lp[14]; /* (skip) */
// type
in_type = strtoul( cp, &cp, 10 );
// tests
in_tests = strtoul( cp, &cp, 10 );
if ( in_tests & 1 )
{
//print("known composite\n");
continue;
}
// tries
in_tries = strtoul( cp, &cp, 10 );
// size (bits)
in_size = strtoul( cp, &cp, 10 );
if ( in_size < 63 )
{
//print("short bit size\n");
continue;
}
// generator (hex)
generator_known = strtoul( cp, &cp, 16 );
// modulus (hex)
switch (in_type)
{
case 4:
//print("Sophie-Germaine\n");
mpz_set_str(&q,cp,16);
/* p = 2*q + 1 */
mpz_mul_2exp(&p,&q,1);
mpz_add_ui(&p,&p,1);
in_size += 1;
generator_known = 0;
break;
default:
mpz_set_str(&p,cp,16);
/* q = (p-1) / 2 */
mpz_div_2exp(&q,&p,1);
break;
};
if ( in_tests < TEST_MILLER_RABIN )
{
in_tries = trials;
}
else
{
in_tries += trials;
}
count_candidates++;
/* guess unknown generator */
if ( generator_known == 0 )
{
if ( mpz_fdiv_ui(&p,24) == 11 )
{
generator_known = 2;
}
else if ( mpz_fdiv_ui(&p,12) == 5 )
{
generator_known = 3;
}
else
{
uint32 r = mpz_fdiv_ui(&p,10);
if ( r == 3 || r == 7 )
{
generator_known = 5;
}
}
}
/* skip tests when desired generator doesn't match */
if ( generator_wanted > 0
&& generator_wanted != generator_known )
{
#ifdef FINALCHECK
printf("generator %ld != %ld\n",
generator_known,
generator_wanted );
#endif
continue;
}
/* The (1/4)^N performance bound on Miller-Rabin is
* extremely pessimistic, so don't spend a lot of time
* really verifying that q is prime until after we know
* that p is also prime. A single pass will weed out the
* vast majority of composite q's.
*/
if(!mpz_probab_prime_p(&q,1))
{
#ifdef FINALCHECK
printf("q failed first possible prime test\n");
#endif
continue;
}
/* q is possibly prime, so go ahead and really make
* sure that p is prime. If it is, then we can go back and
* do the same for q. If p is composite, chances are that
* will show up on the first Rabin-Miller iteration so
* it doesn't hurt to specify a high iteration count.
*/
if(!mpz_probab_prime_p(&p,trials))
{
#ifdef FINALCHECK
printf("p is not prime\n");
#endif
continue;
}
#ifdef FINALCHECK
printf("p is almost certainly prime\n");
#endif
/* recheck q more rigorously */
if(!mpz_probab_prime_p(&q,trials/2))
{
#ifdef FINALCHECK
printf("q is not prime\n");
#endif
continue;
}
#ifdef FINALCHECK
printf("q is almost certainly prime\n");
#endif
time(&time_stop);
gtm = gmtime(&time_stop);
fprintf(finalfile,"%04d%02d%02d%02d%02d%02d ",
gtm->tm_year + 1900,
gtm->tm_mon + 1,
gtm->tm_mday,
gtm->tm_hour,
gtm->tm_min,
gtm->tm_sec );
fprintf(finalfile,"2 "); /* Safe */
fprintf(finalfile,"%lu ",in_tests | TEST_MILLER_RABIN);
fprintf(finalfile,"%lu ",in_tries);
fprintf(finalfile,"%lu ",in_size);
fprintf(finalfile,"%lx ",generator_known);
mpz_out_str(finalfile,16,&p);
fprintf(finalfile,"\n");
count_safe++;
}
time(&time_stop);
free(lp);
fclose(finalfile);
fclose(largefile);
printf("%.24s Found %lu safe primes of %lu candidates in %lu seconds\n",
ctime(&time_stop),
count_safe,
count_candidates,
(long)(time_stop - time_start) );
}
WSimpson@UMich.edu
Key fingerprint = 17 40 5E 67 15 6F 31 26 DD 0D B9 9B 6A 15 2C 32
