/*
 *          Copyright (c) mjh-EDV Beratung, 1996-1999
 *     mjh-EDV Beratung - 63263 Neu-Isenburg - Rosenstrasse 12
 *          Tel +49 6102 328279 - Fax +49 6102 328278
 *                Email info@mjh.teddy-net.com
 *
 *   Author: Jordan Hrycaj <jordan@mjh.teddy-net.com>
 *
 *   $Id: make-primes.c,v 1.1 2000/08/14 20:51:36 jordan Exp $
 *
 *   This library is free software; you can redistribute it and/or
 *   modify it under the terms of the GNU Library General Public
 *   License as published by the Free Software Foundation; either
 *   version 2 of the License, or (at your option) any later version.
 *
 *   This library 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
 *   Library General Public License for more details.
 *
 *   You should have received a copy of the GNU Library General Public
 *   License along with this library; if not, write to the Free
 *   Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include "common-stuff.h"

#if TIME_WITH_SYS_TIME
# include <sys/time.h>
# include <time.h>
#else
# if HAVE_SYS_TIME_H
#  include <sys/time.h>
# else
#  include <time.h>
# endif
#endif

#ifdef USE_PTHREADS
#ifdef HAVE_PTHREAD_H
#include <pthread.h>
#endif
#endif

#ifdef HAVE_ASSERT_H
#include <assert.h>
#else
#define assert(x)
#endif

#ifndef HAVE_RAND
#define rand() random ()
#endif

#include "baseXX.h"
#include "peks-baseXX.h"
#include "rand/crandom.h"
#include "cipher/cipher.h"
#include "rand/rnd-pool.h"

#define __MAKE_PRIME_INTERNAL__
#include "make-primes.h"

/* ------------------------------------------------------------------------ *
 *                private variables                                         *
 * ------------------------------------------------------------------------ */

static int spDIM = 0 ;

/* ------------------------------------------------------------------------ *
 *                private functions                                         *
 * ------------------------------------------------------------------------ */

static void 
xprint 
  (FILE     *stream, 
   const char *text)
{
  fputs (text, stream);
  fflush (stream) ;
}

static unsigned short
get_a_random_smallprime_or_1 
  (void)
{
  extern unsigned short small_prime_numbers [] ;
  short int inx ;

  if (spDIM == 0)
    while (small_prime_numbers [++ spDIM])
      ;

  fast_random_bytes ((char *)&inx, sizeof (inx)) ;

  /* I need the number 1 or a prime from list of small primes 
     starting at 3. */
  inx %= (spDIM + 2) ;

  /* 1 and 2 are not in small_prime_numbers [] */
  inx -= 2 ;

  return inx < 0 ? -inx : small_prime_numbers [inx] ;
}

/* ------------------------------------------------------------------------ *
 *              global functions: i/o call back                             *
 * ------------------------------------------------------------------------ */

void xprint2 (const char *text) {xprint (stderr, text);}
void xprint1 (const char *text) {xprint (stdout, text);}

/* ------------------------------------------------------------------------ *
 *                private functions: pthreads support                       *
 * ------------------------------------------------------------------------ */

#undef  SYSNAME /* solaris wants it */
#define SYSNAME "make-primes"
#include "peks-sema.h"

/* ------------------------------------------------------------------------ *
 *               private functions: random number generator                 *
 * ------------------------------------------------------------------------ */

/* Some algorithm below tests a list of random numbers for some properties.
   Assuming, that the random generator has enough entropy, we do not destill
   a new random number, all the time. We just fiddle around with the given
   random bits, instead. */

typedef struct _hashy { 
  unsigned  length;	/* bytes represented by the rnd string */
  char  buffer [1];	/* random bits */
  /*VARIABLE SIZE*/
} hashy ;


static void
hashy_random_num
  (hashy   **RND,
   mpz_t    *ROP,
   unsigned bits)
{
  static frame_desc *md ;
  hashy *rnd = *RND ;
  unsigned len;
  char *p ;

  __enter_lock_semaphore ();

  if (md == 0) {
    md = create_frame (find_frame_class (MAKE_PRIME_HASH_TYPE), 0);
    assert (md != 0);
  }
  
  if (rnd == 0) {			/* create descriptor */
    len         = (bits + 7) / 8 ;
    rnd         = vmalloc (sizeof (hashy) + len - 1) ;
    rnd->length = len ;
    *RND = rnd ;

    /* get a mouth full of random bits */
    fast_random_bytes (rnd->buffer, rnd->length);
  }
  
  if (ROP == 0) {			/* ROP == 0: clean up */
    memset (rnd, 0, sizeof (hashy)+rnd->length-1) ;
    xfree  (rnd) ;
    __release_lock_semaphore ();
    return ;
  }

  len = rnd->length ;
  p   = rnd->buffer ;

  /* rehash the whole buffer while adding some easy randomness */
  while (len --) {
    int chunk = (md->mdlen < len) ? md->mdlen : len ;
    
    /* hash output from the random generator */
    { int r = rand () ; XCRCFIRST (md, (char*)&r, sizeof (r)); }
    
#ifdef HAVE_GETTIMEOFDAY2
    /* hash current time in millisecs, seconds etc. */
    { struct timeval tv ; gettimeofday2 (&tv, 0);
    XCRCNEXT (md, (char*)&tv, sizeof (tv)); }
#else
    { time_t ti = time (0);
    XCRCNEXT (md, (char*)&ti, sizeof (ti)); }
#endif

    /* re-hash the old contents */
    XCRCNEXT (md, p, chunk) ;

#if 1
    /* quickly add some more bits */
    make_random_bytes (p, chunk);
    XCRCNEXT (md, p, chunk) ;
#endif

    /* copy one chunk of data, set new internal state */
    memcpy (p, XCRCRESULT0 (md), chunk) ;

    /* let a quarter of the fragments overlap */
    len -= (chunk >> 2) * 3 ;
    p   += (chunk >> 2) * 3 ;
  }

  __release_lock_semaphore ();

  if ((len = (bits + 7) / 8) > rnd->length)
    len = rnd->length ;

  /* convert buffer to base 32 */
  bin2mpz (ROP, rnd->buffer, len);

  /* don't let the random number become too small */
  if (mpz_sizeinbase (*ROP, 2) < bits)
    mpz_setbit (*ROP, bits-1) ;
}

/* ------------------------------------------------------------------------ *
 *               private function: prime test                               *
 * ------------------------------------------------------------------------ */

static unsigned
is_not_miller_rabin_prime
  (hashy     **c,
   mpz_t      *p,
   int num_tests)
{
  /* References: 
   *
   *  - Bruce Schneier; Applied Cryprography; NY, Wiley & Sons, 2nd ed 1996, 
   *    chapter 11.5, pg 260 ff.
   *  - D.R. Stinson; Cryptography; Boca Raton, CRC Press 1995, chapter 4,
   *    Fig 4.8 pg 137 ff. 
   */
  
  mpz_t p_min_1, m, a, z;
  unsigned b, is_definitely_composite;

  unsigned rlen = 32 ;			/* max size of a random number */
  mpz_init (a);				/* will a (small) random a < p */

  if ((b = mpz_sizeinbase (*p, 2) - 1) < rlen)
    rlen = b ;				/* length for a randum mumber */

  mpz_init (p_min_1) ;
  mpz_sub_ui (p_min_1, *p, 1) ;
  b = mpz_scan1 (p_min_1, 0) ;		/* max b with 2^b | p-1 */
  
  mpz_init (m);
  mpz_tdiv_q_2exp (m, p_min_1, b);	/* so p-1 = m * 2^b */

  mpz_init (z);
  do {					/* do this test at least once */
    unsigned j ;
      
    is_definitely_composite = 1;	/* default assumption: no prime */
    hashy_random_num (c, &a, rlen);	/* get new random number */

    mpz_powm (z, a, m, *p) ;		/* z := a^m mod p */
    if (mpz_cmp_ui (z, 1) == 0)	/* z == 1 mod p */
      goto p_may_be_prime ;
      
    for (j = 0; j < b; j ++) {		/* see p-1 == m * 2^b */
      if (mpz_cmp (z, p_min_1) ==0)	/* z == -1 mod p */
	goto p_may_be_prime ;
	  
      mpz_powm_ui (z, z, 2, *p) ;	/* z := z^2 mod p */

      POINT_OF_RANDOM_VAR_COND (j % 11 == 0, z);  

      if (mpz_cmp_ui (z, 1) == 0)	/* if z == 1 here, then for rest */
	break ;				/* definitely not prime */
    }
    goto clean_up_and_return ;		/* not a prime when this is done */
    /*@NOTREACHED@*/
    
  p_may_be_prime:
    is_definitely_composite = 0;	/* maybe a prime, next test or done */
  } while ( -- num_tests > 0) ;		/* repeat that test */
  
 clean_up_and_return:

  mpz_clear         (p_min_1);
  mpz_clear               (m);
  mpz_clear               (a);
  mpz_clear               (z);

  return is_definitely_composite ;
}


static unsigned
this_number_is_a_prime
  (hashy            **c,
   mpz_t           *num,
   unsigned prob_weight)
   
{
  extern unsigned short small_prime_numbers [] ;
  unsigned p, r = 0, n = 0 ;

  mpz_t OP;
  mpz_init (OP);

  /* check for divisibility by some small primes, first */
  while ((p = small_prime_numbers [n ++]) != 0) {
    /* does p | num <=> num == 0(mod p) ? */
    r = mpz_mod_ui (OP, *num, (unsigned long)p) ;
    if (r == 0) 
      break ;
  }

  if (r) /* Rabin-Miller test */
    r = (is_not_miller_rabin_prime (c, num, prob_weight) == 0) ;
  
#ifdef _seems_not_to_be_better_than_the_miller_rabin_test_
  if (r) /* some probabilistic test due to Knuth */
    r = mpz_probab_prime_p (*num, prob_weight) ;
#endif

  mpz_clear (OP);
  return r;
}


/* ------------------------------------------------------------------------ *
 *           private function: make module and generator                    *
 * ------------------------------------------------------------------------ */

static unsigned
get_gen_prime_module_for_given_prime
  (hashy          **c,
   mpz_t           *P,			/* prime module to be found */
   unsigned       *gp,			/* generator (mod p) to be found */
   mpz_t           *Q,			/* input: a large prime */
   unsigned     pprob,			/* prime test prob weight */
   unsigned max_tries,
   void (*prt)(const char*))
{
  int run = 0;				/* indices in small_prime_numbers [] */

  unsigned p_is_1_mod_4 ;		/* set, if p == 1 (mod 4) */
  unsigned g ;
  mpz_t G, OP;				/* G := (mpz_t) g */

  mpz_init (G) ;
  mpz_init (OP) ;

  do {
    unsigned r = get_a_random_smallprime_or_1 () ;
    unsigned s = get_a_random_smallprime_or_1 () ;

    unsigned t = r * s ;	/* t = rs, for small primes r, s */

    if (t == 1)
      continue ;

    /* print keep alive message */
    if (prt != 0 && (run) % 44 == 0) (*prt) (".") ;

    /* What we got here, are a product t of small small primes, and a large 
       prime Q.  We check, whether P := (Q * t + 1) is prime.  If so, the
       Euler function 
       
           phi (P) = p - 1 = Q * t 
       
       is the the product of one large, and at most two small primes.  In
       order to show that a given value G, say is a generator of the module
       generated by Q, it suffices to show, that G^Q*r != 1 (mod P), 
       G^Q*s != 1 (mod P), and G^t != 1 (mod P). */
	 
    mpz_mul_ui (OP, *Q, t);			/* OP := q * t */
    mpz_add_ui (*P, OP, 1);			/* p := q * t + 1 */

    /* retry if this number is not prime */
    if (!this_number_is_a_prime (c, P, pprob))
      continue ;

    /* print found message */
    if (prt != 0) (*prt) ("p") ;
    
    /* Check, whether p is 1 (mod 4). There reason for that is an attack on
       the El Gamal signature scheme possible if:
       
           p == 1 (mod 4) and g | (p-1)  ... (more restrictions)

       For details, see Note 11.67 in: A. Menezes, P.van Oorschot,
       S. Vanstone (Eds.), "Handbook of Applied Cryptography", CRC 
       Press 1997, ISBN 0-8493-8523-7, pg 458 ff. */

    p_is_1_mod_4 = (mpz_mod_ui (OP, OP, 4) == 0) ;	/* OP was (p - 1) */

    /* find generator g (mod p) */
    for (g = 2; g < UINT_MAX; g ++) {

      mpz_set_ui (G, g);

      /* check whether g | p-1 (see the comment, above) */
      if (p_is_1_mod_4 && mpz_mod_ui (OP, *P, g) == 1) 
	continue ;

      /* print keep alive message */
      if (prt != 0 && g % 11 == 0) (*prt) (".") ;

      POINT_OF_RANDOM_VAR_COND (g % 13 == 0, OP);

      mpz_powm_ui (OP, G, t, *P) ;		/* OP := g^t (mod p) */
      if (mpz_cmp_ui (OP, 1) == 0)		/* g^t == 1 (mod p) ? */
	continue ;

      /* for convenience we need: s != 1 */
      if (s == 1) { s = r ; r = 1 ; }

      if (r == 1) {
	/* this case saves one multiplication operation */
	mpz_powm (OP, G, *Q, *P) ;
      } else {				/* means r != 1 */
	mpz_mul_ui (OP, *Q, r) ;		/* OP := q * r */
	mpz_powm (OP, G, OP, *P) ;		/* OP := g^qr (mod p) */
      }
      if (mpz_cmp_ui (OP, 1) == 0)		/* g^qr == 1 (mod p) ? */
	continue ;
	
      /* note that s != 1 */
      mpz_mul_ui (OP, *Q, s) ;		/* OP := q * s */
      mpz_powm (OP, G, OP, *P) ;		/* OP := g^qs (mod p) */
      if (mpz_cmp_ui (OP, 1) == 0)		/* g^qr == 1 (mod p) ? */
	continue ;

      /* now, g^x != 1 (mod p) for all d | (p-1), so p is generator */
      if (prt != 0) (*prt) ("g") ;
      mpz_clear (G);
      mpz_clear (OP);
      return *gp = g ;
    }
  } while (++ run < (int)max_tries); 
  
  if (prt != 0) (*prt) (";\n") ;

  mpz_clear (G);
  mpz_clear (OP);

  /* not found: r-code 0 */
  return 0;
}


static unsigned
find_a_random_prime
  (hashy                 **c,
   mpz_t              *prime,
   unsigned             bits,
#if 0 /* experimental, e.g. check 3 (mod 4) */
   int (*cond)(mpz_t*,mpz_t*),
#endif
   unsigned       prob_weight,
   unsigned         max_tries,
   void   (*prt)(const char*))
{
  mpz_t OP ;

  /* random number must be 'bits' bits large */
  int run = max_tries;
  
  mpz_init  (OP);

  do						/* do at least once */
    {
      hashy_random_num (c, prime, bits);	/* get fresh randum number */
      mpz_setbit         (*prime, 0);		/* num |= 1 makes it odd */

      /* print keep alive message, sometimes */
      if (prt != 0 && (run % 20) == 0) (* prt) (".") ;
      
#if 0 /* experimental, e.g. check 3 (mod 4) */
      /* check, whether there is a particular condition */
      if (cond != 0 && (*cond) (prime, &OP) != 0)
	continue ;
#endif
      /* retry if this number is not prime */
      if (!this_number_is_a_prime (c, prime, prob_weight))
	continue ;

      /* print found message, return */
      if (prt != 0) (* prt) ("q") ;

      mpz_clear (OP);
      return 1 + max_tries - run ;
    }
  while (-- run > 0) ;

  if (prt != 0) (* prt) (";\n") ;

  mpz_clear (OP);
  return 0;
}

/* ------------------------------------------------------------------------ *
 *               global functions: make primes                              *
 * ------------------------------------------------------------------------ */

#ifdef USE_PTHREADS
void 
prime_maker_sema_create
  (void *attr)
{
  if (spDIM == 0)
    get_a_random_smallprime_or_1 () ;
  __sema_create (attr); 
}

void 
prime_maker_sema_destroy 
  (void)
{
  __sema_destroy (); 
}
#endif
 
unsigned
number_is_a_prime
  (mpz_t           *num,
   unsigned prob_weight)
{
  hashy *c = 0;
  int n = this_number_is_a_prime (&c, num, prob_weight) ;
  hashy_random_num (&c, 0, 0);
  return n;
}
  
unsigned
get_random_prime
  (mpz_t             *P,
   unsigned        bits,
   unsigned prob_weight,
   unsigned   max_tries,
   void (*prt)(const char*))
{
  hashy *c = 0;
  int n = find_a_random_prime (&c, P, bits, prob_weight, max_tries, prt);
  hashy_random_num (&c, 0, 0);
  return n;
}

void
get_random_num
  (mpz_t    *ROP,
   unsigned bits,
   mpz_t    *MOD)
{
  mpz_t gcd;
  unsigned k;
  hashy *c;
  char *s;

  if (bits == 0)
    bits = 8 ;

  /* find any random number with at least "bits" bits */

  if (MOD == 0) {
    k = (bits + 7) / 8 ;
    s = ALLOCA (k) ;
    prime_random_bytes (s, k) ;
    bin2mpz (ROP, s, k) ;
    DEALLOCA (s) ;
    POINT_OF_RANDOM_VAR (s);
    /* don't let the random number become too small */
    if (mpz_sizeinbase (*ROP, 2) < bits)
      mpz_setbit (*ROP, bits-1) ;
    return ;
    /* @NOTREACHED@ */
  }

  mpz_init (gcd) ;
  do {
    k = 100 ;
    c =   0 ;
    do {
      /* find a random number *ROP relatively prime to *MOD */
      hashy_random_num (&c, ROP, bits);
      mpz_gcd (gcd, *ROP, *MOD) ;
    } while (mpz_cmp_ui (gcd, 1) != 0 && -- k) ;
    hashy_random_num (&c, 0, 0); /* reset */
  } while (k == 0) ;
  
  mpz_clear (gcd) ;  
  POINT_OF_RANDOM_VAR (gcd);
}

unsigned
get_generated_prime_module
  (mpz_t       *P,			/* prime module to be found */
   unsigned   *gp,			/* generator (mod p) to be found */
   mpz_t       *Q,			/* a large prime to be found */
   unsigned  bits,
   void (*prt)(const char*))
{
  int   run = TRY_GENERATE_PRIME_MODULE ;
  int try_p = TRY_RANDOM_IS_PRIME ;
  int try_g = TRY_PRIME_HAS_GENERATOR ;
  int pprob = PRIME_PROBABBILTY_WEIGHT;
  hashy  *c = 0;

  if (bits < 20)
    bits = 20 ;

  if (prt != 0) /* cosmetics */
    (*prt) (GENERATING_PRIMES_TXT) ;

  for (;;) {

    int ok = find_a_random_prime 
      (&c, Q, bits, pprob, try_p, prt) ;
    
    if (ok) 
      ok = get_gen_prime_module_for_given_prime 
	(&c, P, gp, Q, pprob, try_g, prt) ;

    if (ok) {
      /* reset random string recycler */
      hashy_random_num (&c, 0, 0);
      return *gp ;
    }

    if (run -- == 0) { /* loop control */
      hashy_random_num (&c, 0, 0);
      if (prt != 0) /* cosmetics */
	(*prt) (NOMORE_NEXT_TXT) ;
      return 0 ;
    }
    
    if (prt != 0) /* cosmetics */
      (*prt) (ADVANCING_NEXT_TXT) ;

    POINT_OF_RANDOM_VAR (Q);
  }
  /*@NOTREACHED@*/

  /* make stupid compilers happy */
  return 0; 
}