[u/mdw/catacomb] / pfilt.c

/* -*-c-*-
 *
 * $Id: pfilt.c,v 1.2 2000/06/17 11:54:27 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: pfilt.c,v $
 * Revision 1.2  2000/06/17 11:54:27  mdw
 * Use new MP memory management functions.
 *
 * Revision 1.1  1999/12/22 15:49:39  mdw
 * Renamed from `pgen'.  Reworking for new prime-search system.
 *
 * 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 "pfilt.h"
#include "pgen.h"
#include "primetab.h"

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

/* --- @pfilt_create@ --- *
 *
 * Arguments:	@pfilt *p@ = pointer to prime filtering context
 *		@mp *m@ = pointer to initial number to test
 *
 * Returns:	One of the @PGEN@ result codes.
 *
 * 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_TRY@
 *		result is returned, perform Rabin-Miller tests to confirm.
 */

int pfilt_create(pfilt *p, mp *m)
{
  int rc = PGEN_TRY;
  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 = primetab[i];
    mp_div(0, &r, m, &q);
    p->r[i] = r->v[0];
    if (!p->r[i] && rc == PGEN_TRY) {
      if (MP_LEN(m) == 1 && m->v[0] == primetab[i])
	rc = PGEN_DONE;
      else
	rc = PGEN_FAIL;
    }
  }

  /* --- Done --- */

  mp_drop(r);
  return (rc);
}

/* --- @pfilt_destroy@ --- *
 *
 * Arguments:	@pfilt *p@ = pointer to prime filtering context
 *
 * Returns:	---
 *
 * Use:		Discards a context and all the resources it holds.
 */

void pfilt_destroy(pfilt *p)
{
  mp_drop(p->m);
}

/* --- @pfilt_step@ --- *
 *
 * Arguments:	@pfilt *p@ = pointer to prime filtering context
 *		@mpw step@ = how much to step the number
 *
 * Returns:	One of the @PGEN@ result codes.
 *
 * 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 @primetab_create@, so @PGEN_TRY@
 *		results should be followed up by a Rabin-Miller test.
 */

int pfilt_step(pfilt *p, mpw step)
{
  int rc = PGEN_TRY;
  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) % primetab[i];
    if (!p->r[i] && rc == PGEN_TRY) {
      if (MP_LEN(p->m) == 1 && p->m->v[0] == primetab[i])
	rc = PGEN_DONE;
      else
	rc = PGEN_FAIL;
    }
  }

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

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

  /* --- Done --- */

  return (rc);
}

/* --- @pfilt_muladd@ --- *
 *
 * Arguments:	@pfilt *p@ = destination prime filtering context
 *		@const pfilt *q@ = source prime filtering context
 *		@mpw m@ = number to multiply by
 *		@mpw a@ = number to add
 *
 * Returns:	One of the @PGEN@ result codes.
 *
 * Use:		Multiplies the number in a prime filtering 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 pfilt_muladd(pfilt *p, const pfilt *q, mpw m, mpw a)
{
  int rc = PGEN_TRY;
  int i;

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

  {
    mp *d = mp_new(MP_LEN(q->m) + 2, q->m->f);
    mpx_umuln(d->v, d->vl, q->m->v, q->m->vl, m);
    mpx_uaddn(d->v, d->vl, a);
    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) % primetab[i];
    if (!p->r[i] && rc == PGEN_TRY) {
      if (MP_LEN(p->m) == 1 && p->m->v[0] == primetab[i])
	rc = PGEN_DONE;
      else
	rc = PGEN_FAIL;
    }
  }

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

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

  /* --- Finished --- */

  return (rc);
}

/* --- @pfilt_jump@ --- *
 *
 * Arguments:	@pfilt *p@ = pointer to prime filtering context
 *		@const pfilt *j@ = pointer to another filtering context
 *
 * Returns:	One of the @PGEN@ result codes.
 *
 * 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 @primetab_create@, so @PGEN_TRY@
 *		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 pfilt_jump(pfilt *p, const pfilt *j)
{
  int rc = PGEN_TRY;
  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] > primetab[i])
      p->r[i] -= primetab[i];
    if (!p->r[i] && rc == PGEN_TRY) {
      if (MP_LEN(p->m) == 1 && p->m->v[0] == primetab[i])
	rc = PGEN_DONE;
      else
	rc = PGEN_FAIL;
    }
  }

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

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

  /* --- Done --- */

  return (rc);
}

/*----- That's all, folks -------------------------------------------------*/
Commit	Line	Data
9a8b0c8d	1	/* --c--
9a8b0c8d	2	*
87b63d41	3	* $Id: pfilt.c,v 1.2 2000/06/17 11:54:27 mdw Exp $
9a8b0c8d	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: pfilt.c,v $
87b63d41	33	* Revision 1.2 2000/06/17 11:54:27 mdw
	34	* Use new MP memory management functions.
	35	*
9a8b0c8d	36	* Revision 1.1 1999/12/22 15:49:39 mdw
	37	* Renamed from `pgen'. Reworking for new prime-search system.
	38	*
	39	* Revision 1.3 1999/12/10 23:28:35 mdw
	40	* Track suggested destination changes.
	41	*
	42	* Revision 1.2 1999/11/20 22:23:05 mdw
	43	* Add multiply-and-add function for Diffie-Hellman safe prime generation.
	44	*
	45	* Revision 1.1 1999/11/19 13:17:57 mdw
	46	* Prime number generator and tester.
	47	*
	48	*/
	49
	50	/----- Header files ------------------------------------------------------/
	51
	52	#include "mp.h"
	53	#include "mpmont.h"
	54	#include "pfilt.h"
	55	#include "pgen.h"
	56	#include "primetab.h"
	57
	58	/----- Main code ---------------------------------------------------------/
	59
	60	/* --- @pfilt_create@ --- *
	61	*
	62	* Arguments: @pfilt *p@ = pointer to prime filtering context
	63	* @mp *m@ = pointer to initial number to test
	64	*
	65	* Returns: One of the @PGEN@ result codes.
	66	*
	67	* Use: Tests an initial number for primality by computing its
	68	* residue modulo various small prime numbers. This is fairly
	69	* quick, but not particularly certain. If a @PGEN_TRY@
	70	* result is returned, perform Rabin-Miller tests to confirm.
	71	*/
	72
	73	int pfilt_create(pfilt p, mp m)
	74	{
	75	int rc = PGEN_TRY;
	76	int i;
	77	mp *r = MP_NEW;
	78	mpw qw;
	79	mp q;
	80
	81	/* --- Take a copy of the number --- */
	82
	83	mp_shrink(m);
	84	p->m = MP_COPY(m);
	85
	86	/* --- Fill in the residues --- */
	87
	88	mp_build(&q, &qw, &qw + 1);
	89	for (i = 0; i < NPRIME; i++) {
	90	qw = primetab[i];
	91	mp_div(0, &r, m, &q);
	92	p->r[i] = r->v[0];
	93	if (!p->r[i] && rc == PGEN_TRY) {
	94	if (MP_LEN(m) == 1 && m->v[0] == primetab[i])
	95	rc = PGEN_DONE;
	96	else
	97	rc = PGEN_FAIL;
	98	}
	99	}
100
101	/* --- Done --- */
102
103	mp_drop(r);
104	return (rc);
105	}
106
107	/* --- @pfilt_destroy@ --- *
108	*
109	* Arguments: @pfilt *p@ = pointer to prime filtering context
110	*
111	* Returns: ---
112	*
113	* Use: Discards a context and all the resources it holds.
114	*/
115
116	void pfilt_destroy(pfilt *p)
117	{
118	mp_drop(p->m);
119	}
120
121	/* --- @pfilt_step@ --- *
122	*
123	* Arguments: @pfilt *p@ = pointer to prime filtering context
124	* @mpw step@ = how much to step the number
125	*
126	* Returns: One of the @PGEN@ result codes.
127	*
128	* Use: Steps a number by a small amount. Stepping is much faster
129	* than initializing with a new number. The test performed is
130	* the same simple one used by @primetab_create@, so @PGEN_TRY@
131	* results should be followed up by a Rabin-Miller test.
132	*/
133
134	int pfilt_step(pfilt *p, mpw step)
135	{
136	int rc = PGEN_TRY;
137	int i;
138
139	/* --- Add the step on to the number --- */
140
141	p->m = mp_split(p->m);
142	mp_ensure(p->m, MP_LEN(p->m) + 1);
143	mpx_uaddn(p->m->v, p->m->vl, step);
144	mp_shrink(p->m);
145
146	/* --- Update the residue table --- */
147
148	for (i = 0; i < NPRIME; i++) {
149	p->r[i] = (p->r[i] + step) % primetab[i];
150	if (!p->r[i] && rc == PGEN_TRY) {
151	if (MP_LEN(p->m) == 1 && p->m->v[0] == primetab[i])
152	rc = PGEN_DONE;
153	else
154	rc = PGEN_FAIL;
155	}
156	}
157
158	/* --- Small numbers must be prime --- */
159
160	if (rc == PGEN_TRY && MP_LEN(p->m) == 1 &&
161	p->m->v[0] < MAXPRIME * MAXPRIME)
162	rc = PGEN_DONE;
163
164	/* --- Done --- */
165
166	return (rc);
167	}
168
169	/* --- @pfilt_muladd@ --- *
170	*
171	* Arguments: @pfilt *p@ = destination prime filtering context
172	* @const pfilt *q@ = source prime filtering context
173	* @mpw m@ = number to multiply by
174	* @mpw a@ = number to add
175	*
176	* Returns: One of the @PGEN@ result codes.
177	*
178	* Use: Multiplies the number in a prime filtering context by a
179	* small value and then adds a small value. The destination
180	* should either be uninitialized or the same as the source.
181	*
182	* Common things to do include multiplying by 2 and adding 0 to
183	* turn a prime into a jump for finding other primes with @q@ as
184	* a factor of @p - 1@, or multiplying by 2 and adding 1.
185	*/
186
187	int pfilt_muladd(pfilt p, const pfilt q, mpw m, mpw a)
188	{
189	int rc = PGEN_TRY;
190	int i;
191
192	/* --- Multiply the big number --- */
193
194	{
87b63d41	195	mp *d = mp_new(MP_LEN(q->m) + 2, q->m->f);
9a8b0c8d	196	mpx_umuln(d->v, d->vl, q->m->v, q->m->vl, m);
9a8b0c8d	197	mpx_uaddn(d->v, d->vl, a);
9a8b0c8d	198	if (p == q)
	199	mp_drop(p->m);
	200	mp_shrink(d);
	201	p->m = d;
	202	}
	203
	204	/* --- Gallivant through the residue table --- */
	205
	206	for (i = 0; i < NPRIME; i++) {
	207	p->r[i] = (q->r[i] * m + a) % primetab[i];
	208	if (!p->r[i] && rc == PGEN_TRY) {
	209	if (MP_LEN(p->m) == 1 && p->m->v[0] == primetab[i])
	210	rc = PGEN_DONE;
	211	else
	212	rc = PGEN_FAIL;
	213	}
	214	}
	215
	216	/* --- Small numbers must be prime --- */
	217
	218	if (rc == PGEN_TRY && MP_LEN(p->m) == 1 &&
	219	p->m->v[0] < MAXPRIME * MAXPRIME)
	220	rc = PGEN_DONE;
	221
	222	/* --- Finished --- */
	223
	224	return (rc);
	225	}
	226
	227	/* --- @pfilt_jump@ --- *
	228	*
	229	* Arguments: @pfilt *p@ = pointer to prime filtering context
	230	* @const pfilt *j@ = pointer to another filtering context
	231	*
	232	* Returns: One of the @PGEN@ result codes.
	233	*
	234	* Use: Steps a number by a large amount. Even so, jumping is much
	235	* faster than initializing a new number. The test peformed is
	236	* the same simple one used by @primetab_create@, so @PGEN_TRY@
	237	* results should be followed up by a Rabin-Miller test.
	238	*
	239	* Note that the number stored in the @j@ context is probably
	240	* better off being even than prime. The important thing is
	241	* that all of the residues for the number have already been
	242	* computed.
	243	*/
	244
	245	int pfilt_jump(pfilt p, const pfilt j)
	246	{
	247	int rc = PGEN_TRY;
	248	int i;
	249
	250	/* --- Add the step on --- */
	251
	252	p->m = mp_add(p->m, p->m, j->m);
	253
	254	/* --- Update the residue table --- */
	255
	256	for (i = 0; i < NPRIME; i++) {
	257	p->r[i] = p->r[i] + j->r[i];
	258	if (p->r[i] > primetab[i])
	259	p->r[i] -= primetab[i];
	260	if (!p->r[i] && rc == PGEN_TRY) {
	261	if (MP_LEN(p->m) == 1 && p->m->v[0] == primetab[i])
262	rc = PGEN_DONE;
263	else
264	rc = PGEN_FAIL;
265	}
266	}
267
268	/* --- Small numbers must be prime --- */
269
270	if (rc == PGEN_TRY && MP_LEN(p->m) == 1 &&
271	p->m->v[0] < MAXPRIME * MAXPRIME)
272	rc = PGEN_DONE;
273
274	/* --- Done --- */
275
276	return (rc);
277	}
278
279	/----- That's all, folks -------------------------------------------------/