/*
Copyright (C) 2014 Niels Möller
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see .
*/
#include
#include
#include
#include
#include
#include
#include
static void *
xalloc (size_t size)
{
void *p = malloc (size);
if (!p)
{
fprintf (stderr, "Memory exhausted.\n");
abort();
}
return p;
}
static void *
xrealloc (void *p, size_t size)
{
p = realloc (p, size);
if (!p)
{
fprintf (stderr, "Memory exhausted.\n");
abort();
}
return p;
}
struct prime_info
{
uint32_t p; /* An odd prime. */
uint32_t inv; /* p^{-1} mod 2^32 */
uint32_t limit; /* floor ((2^32-1) / p) */
};
static struct prime_info *primes;
static unsigned nprimes;
#define DIVISIBLE_P(x, i) ((x) * primes[i].inv <= primes[i].limit)
static uint32_t binvert(uint32_t a)
{
unsigned char tab[8] = { 1, 11, 13, 7, 9, 3, 5, 15 };
uint32_t x = tab[(a >> 1) & 7]; /* 4-bit inverse */
x = 2*x - a*x*x; /* 8 bits */
x = 2*x - a*x*x; /* 16 bits */
x = 2*x - a*x*x; /* 32 bits */
assert (a*x == 1);
return x;
}
static void
init_table (unsigned limit)
{
/* Number of odd numbers in the range 3 <= j < limit */
unsigned sieve_size = limit/2-1;
char *sieve = xalloc (sieve_size);
unsigned table_size = 0;
unsigned i;
primes = NULL;
nprimes = 0;
memset (sieve, 1, sieve_size);
for (i = 0; i < sieve_size; i++)
{
unsigned p = 3 + 2*i;
unsigned j;
if (sieve[i])
{
if (nprimes >= table_size)
{
table_size = 50 + 3*table_size/2;
primes = xrealloc (primes, table_size * sizeof(*primes));
}
primes[nprimes].p = p;
primes[nprimes].inv = binvert(p);
primes[nprimes].limit = (~(uint32_t) 0) / p;
nprimes++;
}
/* Clear elements corrsponding to p^2, (p+2)^p, ... */
for (j = (p*p-3)/2; j < sieve_size; j += p)
sieve[j] = 0;
}
free (sieve);
}
/* Returns 1 if q + offset and 2(q + offset)+1 appear to be prime */
static int
trial_division (const uint32_t *qmod, const uint32_t *hmod,
uint32_t offset, unsigned n)
{
unsigned i;
for (i = 0; i < nprimes; i++)
{
if (DIVISIBLE_P(qmod[i] + offset * hmod[i], i))
/* q not prime */
return 0;
if (DIVISIBLE_P(2*qmod[i] + 1 + 2*offset * hmod[i], i))
/* p not prime */
return 0;
}
return 1;
}
/* Find the first strong pseudoprime in the sequence start + k step */
static void
next_strong_prime(mpz_t r, const mpz_t start, const mpz_t step)
{
mpz_t t, q0, q, p, pm1, h;
uint32_t *qmod = xalloc (nprimes * sizeof(*qmod));
uint32_t *hmod = xalloc (nprimes * sizeof(*hmod));
unsigned i, j;
unsigned max_i;
mpz_inits (t, q0, q, p, pm1, h, NULL);
/* Some other variants are possible, but exclude them, for
simplicity */
assert (mpz_odd_p (start));
assert (mpz_even_p (step));
mpz_gcd (t, start, step);
assert (mpz_cmp_ui (t, 1) == 0);
mpz_sub_ui (q0, start, 1);
mpz_fdiv_q_2exp (q0, q0, 1);
mpz_fdiv_q_2exp (h, step, 1);
for (i = 0; i < nprimes; i++)
{
hmod[i] = mpz_fdiv_ui (h, primes[i].p);
qmod[i] = mpz_fdiv_ui (q0, primes[i].p);
}
/* Limit to avoid overflow in 2*offest * hmod[j] */
max_i = primes[nprimes-1].limit / 2;
for (;;)
{
for (i = 0; i < max_i; i++)
{
if (!trial_division (qmod, hmod, i, nprimes))
continue;
fprintf (stderr, ".");
mpz_mul_ui (q, h, i);
mpz_add (q, q, q0);
mpz_mul_2exp (pm1, q, 1);
mpz_add_ui (p, pm1, 1);
/* First check if q prime implies p prime. By pocklington,
check that gcd (a^2-1, p) = 1 and a^{2q} = 1 (mod p).
Need use only a == 2, and then we already know that a^2-1
= 3 doesn't have any common factor with p. */
mpz_set_ui (t, 2);
mpz_powm (t, t, q, p);
/* If t = +/- 1 (mod p), then t^2 = 1, and the conditions of
Pocklingtons's theorem (except q prime, which remain to
be checked) are satisfied.
On the other hand, if p and q are indeed prime, then 1 or
-1 are the only possibilities, since t^2 = a^{p-1} = 1
(mod p), and those values are the only square roots of 1.
*/
if (mpz_cmp_ui (t, 1) != 0 && mpz_cmp (t, pm1) != 0)
continue;
fprintf (stderr, "p");
if (!mpz_millerrabin (q, 25))
continue;
fprintf (stderr, "q\n");
assert (mpz_probab_prime_p (p, 25));
mpz_swap (r, p);
mpz_clears (t, q0, q, p, pm1, h, NULL);
free (hmod);
free (qmod);
return;
}
/* Update q0 and qmod */
mpz_addmul_ui (q0, h, i);
for (j = 0; j < nprimes; j++)
{
/* Limit is an approximate reciprocal, floor (2^32 / p) */
uint32_t t = qmod[j] + i * hmod[j];
uint32_t a = ((uint64_t) t * primes[j].limit) >> 32;
qmod[j] = t - a*primes[j].p;
}
}
}
#define BENCH_LIMIT 1000000
#define BENCH_PRIME 0xfffffffb
#define TIME(t, code) do { \
clock_t time_start, time_end; \
unsigned time_reps, time_i; \
for (time_reps = 1;; time_reps *= 2) \
{ \
time_start = clock (); \
for (time_i = 0; time_i < time_reps; time_i++) \
{code;} \
time_end = clock (); \
(t) = (time_end - time_start) / (double) CLOCKS_PER_SEC; \
if ((t) > 0.01) \
break; \
} \
(t) /= time_reps; \
} while (0)
int
main (int argc, char **argv)
{
mpz_t pi8192, start, step, t, p;
int bits;
int limit = 0;
clock_t start_time, end_time;
if (argc < 2)
{
fprintf (stderr,"%s BITS [SIEVE-LIMIT]\n", argv[0]);
return EXIT_FAILURE;
}
bits = atoi(argv[1]);
if (bits < 256 || bits > 8192)
{
fprintf (stderr,"Supported range is 256 <= bits < 8192.\n");
return EXIT_FAILURE;
}
if (argc > 2)
{
limit = atoi(argv[2]);
if (limit < 4)
{
fprintf (stderr,"Invalid prime limit.\n");
return EXIT_FAILURE;
}
}
start_time = clock();
mpz_inits (pi8192, start, step, t, p, NULL);
/* floor (pi * 2^8192), computed with pari/gp. This is an 8194 bit number */
mpz_set_str (pi8192,
"342668632977872056340614340318593589633845496407221537975716"
"613737152546230132865696034849115540419392405119782364448996"
"505597009061206839355038690511597649633340755295662641609317"
"535642924345851483170828622195425552862989984357380546470618"
"654558998494810401062770045070958319118968772694396124916431"
"894580767686315880694138861428875902873600542294826287748558"
"777680794805812784463450507747144885813404423833320387798321"
"609280015466841325521574345612804300925224268961852838014727"
"344949392988773343665477596160680355250188118765396205271825"
"594335340262310201378513530164474118131690724342410799910757"
"546982847641291056024097727491901066776286401345290418249625"
"628084349896086834732113313412019501787434146151912994099903"
"795562822901556466065108423731653397934665628418225503906500"
"252356756723136478156580557983552276992132926102607206902926"
"217226186459755888047793134565096032547470753739485581485056"
"095310312774412642091527392486496373093667048791992749203017"
"039714338536895459352056225848086680604116795009774640833217"
"634173864865403767984673414411303165787683355096108212088449"
"911140581907950288312790601048394648396324285541310288154248"
"299687810226823038493790254875059185044994880417764356982840"
"836688145748032677215052979431605921845872542983666091618301"
"641191123769382721257794537529591295620218822973870974507625"
"710916011007688849347790863162497296922421716445647812915537"
"174656134053117238096786067011739714153581530618764175406523"
"316433003785370745957635127014688173775323182910793894295533"
"869774970638119508426058911918833814428165195037324835873478"
"264560759508465057903722354035637978838385867006178927628223"
"655893247618439468688098957702197415152483332250099523049205"
"947113206854728848346101557957368387613667773778399604898871"
"571620512445391665708436224753527427390770571921354058146861"
"534703545990059875592521055483476922407124171344793745585445"
"612078729092110740281827121644569208704341570355808361384144"
"320068653876547846047724572021904389025572441127188782597462"
"975206221252835136986819772656964303920606775754918430878535"
"831653872047705625068521848678547523189492202180161565293546"
"004893947202674258029340791722950503743099899657901284673761"
"082932857148614112547296837932346608390996266386438042203034"
"975635130970578410146576920597437945752613183601468824949930"
"014430184646230304464605951116975119907446615541881480270560"
"054980165635722725085125496591712178571547448110006296573138"
"030681065709805770745255837670620283490202933751846093959476"
"8886090", 10);
mpz_setbit (step, 64);
mpz_sub_ui (t, step, 1);
mpz_mul_2exp (start, t, bits - 64);
mpz_add (start, start, t);
mpz_fdiv_q_2exp (t, pi8192, (8194 + 128 - bits));
assert (mpz_sizeinbase (t, 2) == bits - 128);
mpz_mul_2exp (t, t, 64);
mpz_add (start, start, t);
if (limit)
init_table(limit);
else
{
/* Try to choose table size automagically */
uint32_t *qmod;
unsigned i;
int res = 0;
double powm_time, div_time;
double ratio;
TIME (powm_time, (mpz_set_ui (t, 2), mpz_powm (t, t, start, start)));
init_table (BENCH_LIMIT);
qmod = xalloc (BENCH_LIMIT * sizeof(*qmod));
for (i = 0; i < nprimes; i++)
qmod[i] = BENCH_PRIME % primes[i].p;
TIME (div_time, res = trial_division (qmod, qmod, 0, nprimes));
assert (res == 1);
free (qmod);
div_time /= nprimes;
ratio = powm_time / div_time;
if (ratio < 10.0)
limit = 10;
else if (ratio > 1e8)
limit = 100000000;
else
limit = (int) ratio;
if (primes[nprimes-1].p < limit)
{
free (primes);
init_table (limit);
}
else
while (primes[nprimes-1].p > limit)
nprimes--;
}
fprintf (stderr, "Using %u primes for sieving, largest p = %u\n",
nprimes, primes[nprimes-1].p);
next_strong_prime (p, start, step);
end_time = clock();
mpz_sub (t, p, start);
mpz_fdiv_q_2exp (t, t, 64);
gmp_printf ("p = 0x%Zx\nk = %Zd\nelapsed time %.2fs\n",
p, t, (end_time - start_time) / (double) CLOCKS_PER_SEC);
mpz_clears (pi8192, start, step, t, p, NULL);
free (primes);
return EXIT_SUCCESS;
}