Re: MODP group modulus derivation [was: Re: I can has SHA-1 hashes for RFC 2409/3526 MODP groups?]

Henrick Hellström <[email protected]> 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.