IETF-SSH archive
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index][Old Index]
Re: [TLS] MODP group modulus derivation [was: Re: I can has SHA-1 hashes for RFC 2409/3526 MODP groups?]
Henrick Hellström <henrick%streamsec.se@localhost> writes:
>> Secondly look at (p-1)/2=q instead. If this is known prime determining
>> if 2*q+1 is prime is quick using Pocklington's Theorem: you don't need
>> to do the extensive primality testing.
>
> Exactly, have you tried this?
I couldn't resist, so I have now tried that. Below is a lite program to
search for strong primes of the given form. It generates a 4096-bit
prime in about 30s on my office machine (an intel core i5). I haven't had
the patience to let it complete search for the 8192 prime.
Some comments:
1. This is a quick hack, not thoroughly tested or profiled, etc.
2. It may be advantageous to use an even larger prime table for the
sieving. Current code overflows somewhere if it is made substantially
larger, though. Generating the prime table is more or less
instantaneous.
3. The 8192-bit prime definitely seems a bit expensive to generate.
4. I let the pari/gp compute the pi approximation for me. The program
could be made more self-contained if one adds some suitable
function for computing pi to arbitrary precision, but I'm not
really familiar with that problem.
Regards,
/Niels
/*
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 <http://www.gnu.org/licenses/>.
*/
#include <assert.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gmp.h>
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])
{
/* fprintf (stderr, "%u\n", p); */
if (nprimes <= table_size)
{
table_size = 50 + 3*table_size/2;
primes = realloc (primes, table_size);
}
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);
}
/* Find the first strong prime in the sequence start + k step */
void
next_strong_prime(mpz_t r, const mpz_t start, const mpz_t step)
{
mpz_t t, q, p, pm1, h;
uint32_t *qmod = xalloc (nprimes * sizeof(*qmod));
uint32_t *hmod = xalloc (nprimes * sizeof(*hmod));
unsigned i;
mpz_inits (t, 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 (q, start, 1);
mpz_fdiv_q_2exp (q, q, 1);
mpz_fdiv_q_2exp (h, step, 1);
for (i = 0; i < nprimes; i++)
hmod[i] = mpz_fdiv_ui (h, primes[i].p);
for (;;)
{
/* New q, recompute qmod */
for (i = 0; i < nprimes; i++)
qmod[i] = mpz_fdiv_ui (q, primes[i].p);
/* Limit to avoid overflow in i * hmod[j] */
for (i = 0; i < 100000; i++)
{
unsigned j;
for (j = 0; j < nprimes; j++)
{
if (DIVISIBLE_P(qmod[j] + i * hmod[j], j))
/* q not prime */
goto not_prime;
if (DIVISIBLE_P(2*qmod[j] + 1 + 2*i * hmod[j], j))
/* p not prime */
goto not_prime;
}
fprintf (stderr, ".");
mpz_addmul_ui (q, h, i);
mpz_mul_2exp (pm1, q, 1);
mpz_add_ui (p, pm1, 1);
/* Since q is updated, we need to reset i and update qmod */
for (j = 0; j < nprimes; j++)
qmod[j] = (qmod[j] + i * hmod[j]) % primes[j].p;
i = 0;
/* 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)
goto not_prime;
fprintf (stderr, "p");
if (!mpz_millerrabin (q, 25))
goto not_prime;
assert (mpz_probab_prime_p (p, 25));
fprintf (stderr, "pq\n");
mpz_swap (r, p);
mpz_clears (t, q, p, pm1, h, NULL);
free (hmod);
free (qmod);
return;
not_prime:
/* Go on with the same i and q */
;
}
/* Can this ever happen? */
mpz_addmul_ui (q, h, i);
fprintf (stderr, "100000 candidates passed trial division!\n");
}
}
int
main (int argc, char **argv)
{
mpz_t pi8192, start, step, t, p;
int bits;
int limit = 1000000;
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;
}
}
init_table(limit);
fprintf (stderr, "Using %u primes for sieving, largest p = %u\n",
nprimes, primes[nprimes-1].p);
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);
next_strong_prime (p, start, step);
mpz_sub (t, p, start);
mpz_fdiv_q_2exp (t, t, 64);
gmp_printf ("p = 0x%Zx\nt = %Zd\n", p, t);
mpz_clears (pi8192, start, step, t, p, NULL);
free (primes);
return EXIT_SUCCESS;
}
--
Niels Möller. PGP-encrypted email is preferred. Keyid C0B98E26.
Internet email is subject to wholesale government surveillance.
Home |
Main Index |
Thread Index |
Old Index