Prime number generator and tester.

[u/mdw/catacomb] / pgen.c

/* -*-c-*-
 *
 * $Id: pgen.c,v 1.1 1999/11/19 13:17:57 mdw Exp $
 *
 * Finding and testing prime numbers
 *
 * (c) 1999 Straylight/Edgeware
 */

/*----- Licensing notice --------------------------------------------------* 
 *
 * This file is part of Catacomb.
 *
 * Catacomb 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.
 * 
 * Catacomb 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 Catacomb; if not, write to the Free
 * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 * MA 02111-1307, USA.
 */

/*----- Revision history --------------------------------------------------* 
 *
 * $Log: pgen.c,v $
 * Revision 1.1  1999/11/19 13:17:57  mdw
 * Prime number generator and tester.
 *
 */

/*----- Header files ------------------------------------------------------*/

#include "mp.h"
#include "mpmont.h"
#include "pgen.h"
#include "ptab.h"
#include "rabin.h"

/*----- Main code ---------------------------------------------------------*/

/* --- @pgen_create@ --- *
 *
 * Arguments:	@pgen *p@ = pointer to prime generation context
 *		@mp *m@ = pointer to initial number to test
 *
 * Returns:	One of the @PGEN@ constants above.
 *
 * Use:		Tests an initial number for primality by computing its
 *		residue modulo various small prime numbers.  This is fairly
 *		quick, but not particularly certain.  If a @PGEN_MAYBE@
 *		result is returned, perform Rabin-Miller tests to confirm.
 */

int pgen_create(pgen *p, mp *m)
{
  int rc = PGEN_MAYBE;
  int i;
  mp *r = MP_NEW;
  mpw qw;
  mp q;

  /* --- Take a copy of the number --- */

  mp_shrink(m);
  p->m = MP_COPY(m);

  /* --- Fill in the residues --- */

  mp_build(&q, &qw, &qw + 1);
  for (i = 0; i < NPRIME; i++) {
    qw = ptab[i];
    mp_div(0, &r, m, &q);
    p->r[i] = r->v[0];
    if (!p->r[i] && rc == PGEN_MAYBE) {
      if (MP_LEN(m) == 1 && m->v[0] == ptab[i])
	rc = PGEN_PRIME;
      else
	rc = PGEN_COMPOSITE;
    }
  }

  /* --- Done --- */

  mp_drop(r);
  return (rc);
}

/* --- @pgen_destroy@ --- *
 *
 * Arguments:	@pgen *p@ = pointer to prime generation context
 *
 * Returns:	---
 *
 * Use:		Discards a context and all the resources it holds.
 */

void pgen_destroy(pgen *p)
{
  mp_drop(p->m);
}

/* --- @pgen_step@ --- *
 *
 * Arguments:	@pgen *p@ = pointer to prime generation context
 *		@mpw step@ = how much to step the number
 *
 * Returns:	One of the @PGEN@ constants above.
 *
 * Use:		Steps a number by a small amount.  Stepping is much faster
 *		than initializing with a new number.  The test performed is
 *		the same simple one used by @ptab_create@, so @PGEN_MAYBE@
 *		results should be followed up by a Rabin-Miller test.
 */

int pgen_step(pgen *p, mpw step)
{
  mp s;
  int rc = PGEN_MAYBE;
  int i;

  /* --- Add the step on to the number --- */

  mp_build(&s, &step, &step + 1);
  p->m = mp_add(p->m, p->m, &s);

  /* --- Update the residue table --- */

  for (i = 0; i < NPRIME; i++) {
    p->r[i] = (p->r[i] + step) % ptab[i];
    if (!p->r[i] && rc == PGEN_MAYBE) {
      if (MP_LEN(p->m) == 1 && p->m->v[0] == ptab[i])
	rc = PGEN_PRIME;
      else
	rc = PGEN_COMPOSITE;
    }
  }

  /* --- Done --- */

  return (rc);
}

/* --- @pgen_jump@ --- *
 *
 * Arguments:	@pgen *p@ = pointer to prime generation context
 *		@pgen *j@ = pointer to another generation context
 *
 * Returns:	One of the @PGEN@ constants above.
 *
 * Use:		Steps a number by a large amount.  Even so, jumping is much
 *		faster than initializing a new number.  The test peformed is
 *		the same simple one used by @ptab_create@, so @PGEN_MAYBE@
 *		results should be followed up by a Rabin-Miller test.
 *
 *		Note that the number stored in the @j@ context is probably
 *		better off being even than prime.  The important thing is
 *		that all of the residues for the number have already been
 *		computed.
 */

int pgen_jump(pgen *p, pgen *j)
{
  int rc = PGEN_MAYBE;
  int i;

  /* --- Add the step on --- */

  p->m = mp_add(p->m, p->m, j->m);

  /* --- Update the residue table --- */

  for (i = 0; i < NPRIME; i++) {
    p->r[i] = p->r[i] + j->r[i];
    if (p->r[i] > ptab[i])
      p->r[i] -= ptab[i];
    if (!p->r[i] && rc == PGEN_MAYBE) {
      if (MP_LEN(p->m) == 1 && p->m->v[0] == ptab[i])
	rc = PGEN_PRIME;
      else
	rc = PGEN_COMPOSITE;
    }
  }

  /* --- Done --- */

  return (rc);
}

/*----- Test code ---------------------------------------------------------*/

#ifdef TEST_RIG

#include <mLib/testrig.h>

static int verify(dstr *v)
{
  mp *m = *(mp **)v[0].buf;
  mp *r = *(mp **)v[1].buf;
  pgen p;
  mpw gw;
  mp g;
  int rc;
  static char baton[] = "/-\\|";
  char *b = baton;
  int ok = 1;

  mp_build(&g, &gw, &gw + 1);
  rc = pgen_create(&p, m);
  for (;;) {
    if (rc == PGEN_PRIME)
      break;
    if (rc == PGEN_MAYBE) {
      int i;
      rabin r;
      rabin_create(&r, p.m);
      for (i = 0; i < 5; i++) {
	gw = rand();
	if (rabin_test(&r, &g) == PGEN_COMPOSITE)
	  break;
      }
      rabin_destroy(&r);
      if (i == 5)
	break;
      putc(*b++, stderr);
      putc('\b', stderr);
      if (!*b)
	b = baton;
    }
    rc = pgen_step(&p, 2);
  }

  if (MP_CMP(p.m, !=, r)) {
    fputs("\n*** failed.", stderr);
    fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
    fputs("\np = ", stderr); mp_writefile(p.m, stderr, 10);
    fputs("\nr = ", stderr); mp_writefile(r, stderr, 10);
    fputc('\n', stderr);
    ok = 0;
  }

  mp_drop(m);
  mp_drop(r);
  pgen_destroy(&p);
  return (ok);
}

static test_chunk tests[] = {
  { "pgen", verify, { &type_mp, &type_mp, 0 } },
  { 0, 0, { 0 } }
};

int main(int argc, char *argv[])
{
  sub_init();
  test_run(argc, argv, tests, SRCDIR "/tests/pgen");
  return (0);
}

#endif

/*----- That's all, folks -------------------------------------------------*/