[u/mdw/catacomb] / pgen.c

/* -*-c-*-
 *
 * $Id: pgen.c,v 1.3 1999/12/10 23:28:35 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.3  1999/12/10 23:28:35  mdw
 * Track suggested destination changes.
 *
 * Revision 1.2  1999/11/20 22:23:05  mdw
 * Add multiply-and-add function for Diffie-Hellman safe prime generation.
 *
 * 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)
{
  int rc = PGEN_MAYBE;
  int i;

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

  p->m = mp_split(p->m);
  mp_ensure(p->m, MP_LEN(p->m) + 1);
  mpx_uaddn(p->m->v, p->m->vl, step);
  mp_shrink(p->m);

  /* --- 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;
    }
  }

  /* --- Small numbers must be prime --- */

  if (rc == PGEN_MAYBE && MP_LEN(p->m) == 1 &&
      p->m->v[0] < MAXPRIME * MAXPRIME)
    rc = PGEN_PRIME;

  /* --- Done --- */

  return (rc);
}

/* --- @pgen_muladd@ --- *
 *
 * Arguments:	@pgen *p@ = destination prime generation context
 *		@const pgen *q@ = source prime generation context
 *		@mpw m@ = number to multiply by
 *		@mpw a@ = number to add
 *
 * Returns:	One of the @PGEN@ constants above.
 *
 * Use:		Multiplies the number in a prime generation context by a
 *		small value and then adds a small value.  The destination
 *		should either be uninitialized or the same as the source.
 *
 *		Common things to do include multiplying by 2 and adding 0 to
 *		turn a prime into a jump for finding other primes with @q@ as
 *		a factor of @p - 1@, or multiplying by 2 and adding 1.
 */

int pgen_muladd(pgen *p, const pgen *q, mpw m, mpw a)
{
  int rc = PGEN_MAYBE;
  int i;

  /* --- Multiply the big number --- */

  {
    mp *d = mp_create(MP_LEN(q->m) + 2);
    mpx_umuln(d->v, d->vl, q->m->v, q->m->vl, m);
    mpx_uaddn(d->v, d->vl, a);
    d->f = q->m->f;
    if (p == q)
      mp_drop(p->m);
    mp_shrink(d);
    p->m = d;
  }

  /* --- Gallivant through the residue table --- */
      
  for (i = 0; i < NPRIME; i++) {
    p->r[i] = (q->r[i] * m + a) % 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;
    }
  }

  /* --- Small numbers must be prime --- */

  if (rc == PGEN_MAYBE && MP_LEN(p->m) == 1 &&
      p->m->v[0] < MAXPRIME * MAXPRIME)
    rc = PGEN_PRIME;

  /* --- Finished --- */

  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;
    }
  }

  /* --- Small numbers must be prime --- */

  if (rc == PGEN_MAYBE && MP_LEN(p->m) == 1 &&
      p->m->v[0] < MAXPRIME * MAXPRIME)
    rc = PGEN_PRIME;

  /* --- 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);
  assert(mparena_count(MPARENA_GLOBAL) == 0);
  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 -------------------------------------------------*/
Commit	Line	Data
0f5ec153	1	/* --c--
0f5ec153	2	*
d94f85ac	3	* $Id: pgen.c,v 1.3 1999/12/10 23:28:35 mdw Exp $
0f5ec153	4	*
	5	* Finding and testing prime numbers
	6	*
	7	* (c) 1999 Straylight/Edgeware
	8	*/
	9
	10	/----- Licensing notice --------------------------------------------------
	11	*
	12	* This file is part of Catacomb.
	13	*
	14	* Catacomb is free software; you can redistribute it and/or modify
	15	* it under the terms of the GNU Library General Public License as
	16	* published by the Free Software Foundation; either version 2 of the
	17	* License, or (at your option) any later version.
	18	*
	19	* Catacomb is distributed in the hope that it will be useful,
	20	* but WITHOUT ANY WARRANTY; without even the implied warranty of
	21	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	22	* GNU Library General Public License for more details.
	23	*
	24	* You should have received a copy of the GNU Library General Public
	25	* License along with Catacomb; if not, write to the Free
	26	* Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
	27	* MA 02111-1307, USA.
	28	*/
	29
	30	/----- Revision history --------------------------------------------------
	31	*
	32	* $Log: pgen.c,v $
d94f85ac	33	* Revision 1.3 1999/12/10 23:28:35 mdw
	34	* Track suggested destination changes.
	35	*
bd98b2df	36	* Revision 1.2 1999/11/20 22:23:05 mdw
	37	* Add multiply-and-add function for Diffie-Hellman safe prime generation.
	38	*
0f5ec153	39	* Revision 1.1 1999/11/19 13:17:57 mdw
	40	* Prime number generator and tester.
	41	*
	42	*/
	43
	44	/----- Header files ------------------------------------------------------/
	45
	46	#include "mp.h"
	47	#include "mpmont.h"
	48	#include "pgen.h"
	49	#include "ptab.h"
	50	#include "rabin.h"
	51
	52	/----- Main code ---------------------------------------------------------/
	53
	54	/* --- @pgen_create@ --- *
	55	*
	56	* Arguments: @pgen *p@ = pointer to prime generation context
	57	* @mp *m@ = pointer to initial number to test
	58	*
	59	* Returns: One of the @PGEN@ constants above.
	60	*
	61	* Use: Tests an initial number for primality by computing its
	62	* residue modulo various small prime numbers. This is fairly
	63	* quick, but not particularly certain. If a @PGEN_MAYBE@
	64	* result is returned, perform Rabin-Miller tests to confirm.
	65	*/
	66
	67	int pgen_create(pgen p, mp m)
	68	{
	69	int rc = PGEN_MAYBE;
	70	int i;
	71	mp *r = MP_NEW;
	72	mpw qw;
	73	mp q;
	74
	75	/* --- Take a copy of the number --- */
	76
	77	mp_shrink(m);
	78	p->m = MP_COPY(m);
	79
	80	/* --- Fill in the residues --- */
	81
	82	mp_build(&q, &qw, &qw + 1);
	83	for (i = 0; i < NPRIME; i++) {
	84	qw = ptab[i];
	85	mp_div(0, &r, m, &q);
	86	p->r[i] = r->v[0];
	87	if (!p->r[i] && rc == PGEN_MAYBE) {
	88	if (MP_LEN(m) == 1 && m->v[0] == ptab[i])
	89	rc = PGEN_PRIME;
	90	else
	91	rc = PGEN_COMPOSITE;
	92	}
	93	}
	94
	95	/* --- Done --- */
	96
	97	mp_drop(r);
	98	return (rc);
	99	}
	100
	101	/* --- @pgen_destroy@ --- *
	102	*
103	* Arguments: @pgen *p@ = pointer to prime generation context
104	*
105	* Returns: ---
106	*
107	* Use: Discards a context and all the resources it holds.
108	*/
109
110	void pgen_destroy(pgen *p)
111	{
112	mp_drop(p->m);
113	}
114
115	/* --- @pgen_step@ --- *
116	*
117	* Arguments: @pgen *p@ = pointer to prime generation context
118	* @mpw step@ = how much to step the number
119	*
120	* Returns: One of the @PGEN@ constants above.
121	*
122	* Use: Steps a number by a small amount. Stepping is much faster
123	* than initializing with a new number. The test performed is
124	* the same simple one used by @ptab_create@, so @PGEN_MAYBE@
125	* results should be followed up by a Rabin-Miller test.
126	*/
127
128	int pgen_step(pgen *p, mpw step)
129	{
0f5ec153	130	int rc = PGEN_MAYBE;
	131	int i;
	132
	133	/* --- Add the step on to the number --- */
	134
d94f85ac	135	p->m = mp_split(p->m);
d94f85ac	136	mp_ensure(p->m, MP_LEN(p->m) + 1);
bd98b2df	137	mpx_uaddn(p->m->v, p->m->vl, step);
bd98b2df	138	mp_shrink(p->m);
0f5ec153	139
	140	/* --- Update the residue table --- */
	141
	142	for (i = 0; i < NPRIME; i++) {
	143	p->r[i] = (p->r[i] + step) % ptab[i];
	144	if (!p->r[i] && rc == PGEN_MAYBE) {
	145	if (MP_LEN(p->m) == 1 && p->m->v[0] == ptab[i])
	146	rc = PGEN_PRIME;
	147	else
	148	rc = PGEN_COMPOSITE;
	149	}
	150	}
	151
bd98b2df	152	/* --- Small numbers must be prime --- */
	153
	154	if (rc == PGEN_MAYBE && MP_LEN(p->m) == 1 &&
	155	p->m->v[0] < MAXPRIME * MAXPRIME)
	156	rc = PGEN_PRIME;
	157
0f5ec153	158	/* --- Done --- */
	159
	160	return (rc);
	161	}
	162
bd98b2df	163	/* --- @pgen_muladd@ --- *
	164	*
	165	* Arguments: @pgen *p@ = destination prime generation context
	166	* @const pgen *q@ = source prime generation context
	167	* @mpw m@ = number to multiply by
	168	* @mpw a@ = number to add
	169	*
	170	* Returns: One of the @PGEN@ constants above.
	171	*
	172	* Use: Multiplies the number in a prime generation context by a
	173	* small value and then adds a small value. The destination
	174	* should either be uninitialized or the same as the source.
	175	*
	176	* Common things to do include multiplying by 2 and adding 0 to
	177	* turn a prime into a jump for finding other primes with @q@ as
	178	* a factor of @p - 1@, or multiplying by 2 and adding 1.
	179	*/
	180
	181	int pgen_muladd(pgen p, const pgen q, mpw m, mpw a)
	182	{
	183	int rc = PGEN_MAYBE;
	184	int i;
	185
	186	/* --- Multiply the big number --- */
	187
	188	{
	189	mp *d = mp_create(MP_LEN(q->m) + 2);
	190	mpx_umuln(d->v, d->vl, q->m->v, q->m->vl, m);
	191	mpx_uaddn(d->v, d->vl, a);
	192	d->f = q->m->f;
	193	if (p == q)
	194	mp_drop(p->m);
	195	mp_shrink(d);
	196	p->m = d;
	197	}
	198
	199	/* --- Gallivant through the residue table --- */
	200
	201	for (i = 0; i < NPRIME; i++) {
	202	p->r[i] = (q->r[i] * m + a) % ptab[i];
	203	if (!p->r[i] && rc == PGEN_MAYBE) {
	204	if (MP_LEN(p->m) == 1 && p->m->v[0] == ptab[i])
	205	rc = PGEN_PRIME;
	206	else
	207	rc = PGEN_COMPOSITE;
	208	}
	209	}
	210
	211	/* --- Small numbers must be prime --- */
	212
	213	if (rc == PGEN_MAYBE && MP_LEN(p->m) == 1 &&
	214	p->m->v[0] < MAXPRIME * MAXPRIME)
	215	rc = PGEN_PRIME;
	216
	217	/* --- Finished --- */
	218
	219	return (rc);
	220	}
	221
	222
0f5ec153	223	/* --- @pgen_jump@ --- *
	224	*
	225	* Arguments: @pgen *p@ = pointer to prime generation context
	226	* @pgen *j@ = pointer to another generation context
	227	*
	228	* Returns: One of the @PGEN@ constants above.
	229	*
	230	* Use: Steps a number by a large amount. Even so, jumping is much
	231	* faster than initializing a new number. The test peformed is
	232	* the same simple one used by @ptab_create@, so @PGEN_MAYBE@
	233	* results should be followed up by a Rabin-Miller test.
	234	*
	235	* Note that the number stored in the @j@ context is probably
	236	* better off being even than prime. The important thing is
	237	* that all of the residues for the number have already been
	238	* computed.
	239	*/
	240
	241	int pgen_jump(pgen p, pgen j)
	242	{
	243	int rc = PGEN_MAYBE;
	244	int i;
	245
	246	/* --- Add the step on --- */
	247
	248	p->m = mp_add(p->m, p->m, j->m);
	249
	250	/* --- Update the residue table --- */
	251
	252	for (i = 0; i < NPRIME; i++) {
	253	p->r[i] = p->r[i] + j->r[i];
	254	if (p->r[i] > ptab[i])
	255	p->r[i] -= ptab[i];
	256	if (!p->r[i] && rc == PGEN_MAYBE) {
	257	if (MP_LEN(p->m) == 1 && p->m->v[0] == ptab[i])
	258	rc = PGEN_PRIME;
	259	else
	260	rc = PGEN_COMPOSITE;
	261	}
	262	}
	263
bd98b2df	264	/* --- Small numbers must be prime --- */
	265
	266	if (rc == PGEN_MAYBE && MP_LEN(p->m) == 1 &&
	267	p->m->v[0] < MAXPRIME * MAXPRIME)
	268	rc = PGEN_PRIME;
	269
0f5ec153	270	/* --- Done --- */
	271
	272	return (rc);
	273	}
	274
	275	/----- Test code ---------------------------------------------------------/
	276
	277	#ifdef TEST_RIG
	278
	279	#include <mLib/testrig.h>
	280
	281	static int verify(dstr *v)
	282	{
	283	mp m = (mp **)v[0].buf;
	284	mp r = (mp **)v[1].buf;
	285	pgen p;
	286	mpw gw;
	287	mp g;
	288	int rc;
	289	static char baton[] = "/-\\\|";
	290	char *b = baton;
	291	int ok = 1;
	292
	293	mp_build(&g, &gw, &gw + 1);
	294	rc = pgen_create(&p, m);
	295	for (;;) {
	296	if (rc == PGEN_PRIME)
	297	break;
	298	if (rc == PGEN_MAYBE) {
	299	int i;
	300	rabin r;
	301	rabin_create(&r, p.m);
	302	for (i = 0; i < 5; i++) {
	303	gw = rand();
	304	if (rabin_test(&r, &g) == PGEN_COMPOSITE)
	305	break;
	306	}
	307	rabin_destroy(&r);
	308	if (i == 5)
	309	break;
	310	putc(*b++, stderr);
	311	putc('\b', stderr);
	312	if (!*b)
	313	b = baton;
	314	}
	315	rc = pgen_step(&p, 2);
	316	}
	317
	318	if (MP_CMP(p.m, !=, r)) {
	319	fputs("\n*** failed.", stderr);
	320	fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
	321	fputs("\np = ", stderr); mp_writefile(p.m, stderr, 10);
	322	fputs("\nr = ", stderr); mp_writefile(r, stderr, 10);
	323	fputc('\n', stderr);
	324	ok = 0;
	325	}
	326
	327	mp_drop(m);
	328	mp_drop(r);
	329	pgen_destroy(&p);
d94f85ac	330	assert(mparena_count(MPARENA_GLOBAL) == 0);
0f5ec153	331	return (ok);
	332	}
	333
	334	static test_chunk tests[] = {
	335	{ "pgen", verify, { &type_mp, &type_mp, 0 } },
	336	{ 0, 0, { 0 } }
	337	};
	338
	339	int main(int argc, char *argv[])
	340	{
	341	sub_init();
	342	test_run(argc, argv, tests, SRCDIR "/tests/pgen");
	343	return (0);
	344	}
	345
	346	#endif
	347
	348	/----- That's all, folks -------------------------------------------------/