[catacomb] / math / pgen.c

/* -*-c-*-
 *
 * Prime generation glue
 *
 * (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.
 */

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

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "fibrand.h"
#include "grand.h"
#include "mp.h"
#include "mprand.h"
#include "pgen.h"
#include "pfilt.h"
#include "rabin.h"

/*----- Standard prime filter ---------------------------------------------*/

/* --- @pgen_filter@ --- */

int pgen_filter(int rq, pgen_event *ev, void *p)
{
  pgen_filterctx *f = p;
  int rc = PGEN_FAIL;

  switch (rq) {
    case PGEN_BEGIN:
      rc = pfilt_create(&f->f, ev->m);
      mp_drop(ev->m);
      break;
    case PGEN_TRY:
      mp_drop(ev->m);
      break;
    case PGEN_DONE:
      pfilt_destroy(&f->f);
      return (PGEN_DONE);
    default:
      rc = PGEN_ABORT;
      break;
  }

  if (rc == PGEN_FAIL && !((f->step | f->f.m->v[0]) & 1))
    rc = pfilt_step(&f->f, 1);
  while (rc == PGEN_FAIL)
    rc = pfilt_step(&f->f, f->step);
  ev->m = MP_COPY(f->f.m);
  return (rc);
}

/* --- @pgen_jump@ --- *
 *
 * Similar to the standard @pgen_filter@, but jumps in large steps rather
 * than small ones.
 */

int pgen_jump(int rq, pgen_event *ev, void *p)
{
  pgen_jumpctx *f = p;
  int rc = PGEN_ABORT;

  switch (rq) {
    case PGEN_BEGIN: {
      mp *g = MP_NEW;
      mp_gcd(&g, 0, 0, ev->m, f->j->m);
      if (MP_CMP(g, >, MP_ONE)) {
	mp_drop(g);
	return (PGEN_ABORT);
      }
      mp_drop(g);
      rc = pfilt_create(&f->f, ev->m);
      mp_drop(ev->m);
    } break;
    case PGEN_TRY:
      mp_drop(ev->m);
      rc = pfilt_jump(&f->f, f->j);
      break;
    case PGEN_DONE:
      pfilt_destroy(&f->f);
      return (PGEN_DONE);
  }

  while (rc == PGEN_FAIL)
    rc = pfilt_jump(&f->f, f->j);
  ev->m = MP_COPY(f->f.m);
  return (rc);
}

/*----- Standard prime test -----------------------------------------------*/

/* --- @pgen_test@ --- */

int pgen_test(int rq, pgen_event *ev, void *p)
{
  rabin *r = p;
  mp *a = MP_NEW;
  int rc = PGEN_ABORT;

  switch (rq) {
    case PGEN_BEGIN:
      rabin_create(r, ev->m);
      rc = PGEN_TRY;
      break;
    case PGEN_TRY:
      a = mprand_range(a, ev->m, ev->r, 0);
      rc = rabin_rtest(r, a);
      break;
    case PGEN_DONE:
      rabin_destroy(r);
      rc = PGEN_DONE;
      break;
  }

  mp_drop(a);
  return (rc);
}

/* --- @pgen_bailliepswtest@ --- */

int pgen_bailliepswtest(int rq, pgen_event *ev, void *p)
{
  rabin r;
  int rc;

  switch (rq) {
    case PGEN_BEGIN:
      if (ev->tests != 2) rc = PGEN_ABORT;
      else rc = PGEN_TRY;
      break;

    case PGEN_DONE:
      rc = PGEN_DONE;
      break;

    case PGEN_TRY:
      switch (ev->tests) {
	case 2:
	  rabin_create(&r, ev->m);
	  rc = rabin_test(&r, MP_TWO);
	  rabin_destroy(&r);
	  break;
	case 1:
	  rc = pgen_granfrob(ev->m, 0, 0);
	  break;
	default:
	  rc = PGEN_ABORT;
	  break;
      }
      break;

    default:
      rc = PGEN_ABORT;
      break;
  }

  return (rc);
}

/*----- The main driver ---------------------------------------------------*/

/* --- @pgen@ --- *
 *
 * Arguments:	@const char *name@ = name of the value being searched for
 *		@mp *d@ = destination for the result integer
 *		@mp *m@ = start value to pass to stepper
 *		@pgen_proc *event@ = event handler function
 *		@void *ectx@ = context argument for event andler
 *		@unsigned steps@ = number of steps to take in search
 *		@pgen_proc *step@ = stepper function to use
 *		@void *sctx@ = context argument for stepper
 *		@unsigned tests@ = number of tests to make
 *		@pgen_proc *test@ = tester function to use
 *		@void *tctx@ = context argument for tester
 *
 * Returns:	Pointer to final result, or null.
 *
 * Use:		A generalized prime-number search skeleton.  Yes, that's a
 *		scary number of arguments.
 */

mp *pgen(const char *name, mp *d, mp *m, pgen_proc *event, void *ectx,
	 unsigned steps, pgen_proc *step, void *sctx,
	 unsigned tests, pgen_proc *test, void *tctx)
{
  pgen_event ev;
  int rq, rc;
  pgen_proc *proc;
  void *ctx;
  int p;

  enum { P_STEP, P_TEST };

  /* --- Set up the initial event block --- */

  ev.name = name;
  if (m)
    ev.m = MP_COPY(m);
  else
    ev.m = 0;
  ev.steps = steps;
  ev.tests = tests;
  ev.r = fibrand_create(0);

  /* --- Tell the event handler we're under way --- */

  if (event && event(PGEN_BEGIN, &ev, ectx) == PGEN_ABORT) {
    ev.r->ops->destroy(ev.r);
    return (0);
  }

  /* --- Set up for the initial call --- */

  proc = step; ctx = sctx; p = P_STEP; rq = PGEN_BEGIN;

  /* --- Enter the great maelstrom of state transitions --- */

  for (;;) {
    unsigned act = 0;

#define A_STEP 1u
#define A_TEST 2u
#define A_EVENT 4u
#define A_ENDTEST 8u
#define A_ENDSTEP 16u
#define A_DONE 32u

    /* --- Call the procedure and decide what to do next --- */

    rc = proc(rq, &ev, ctx);
    switch (rc) {
      case PGEN_TRY:
	if (p == P_TEST)
	  rq = PGEN_TRY;
	else {
	  act |= A_EVENT;
	  proc = test; ctx = tctx; p = P_TEST;
	  rq = PGEN_BEGIN;
	}
	break;
      case PGEN_PASS:
	act |= A_TEST | A_EVENT;
	if (p == P_TEST)
	  rq = PGEN_TRY;
	else {
	  proc = test; ctx = tctx; p = P_TEST;
	  rq = PGEN_BEGIN;
	}
	break;
      case PGEN_FAIL:
	act |= A_STEP;
	if (p == P_TEST) {
	  act |= A_ENDTEST | A_EVENT;
	  proc = step; ctx = sctx; p = P_STEP;
	}
	rq = PGEN_TRY;
	break;
      case PGEN_DONE:
	act |= A_EVENT | A_DONE | A_ENDSTEP;
	if (p == P_TEST)
	  act |= A_ENDTEST;
	break;
      case PGEN_ABORT:
	act |= A_EVENT | A_DONE;
	if (p == P_TEST || rq == PGEN_TRY)
	  act |= A_ENDSTEP;
	if (p == P_TEST && rq != PGEN_BEGIN)
	  act |= A_ENDTEST;
	break;
      default:
	assert(((void)"Invalid response from function", 0));
	break;
    }

    /* --- If decrementing counters is requested, do that --- */

    if ((act & A_STEP) && steps) {
      ev.steps--;
      if (!ev.steps) {
	act |= A_EVENT | A_ENDSTEP | A_DONE;
	rc = PGEN_ABORT;
      }
      ev.tests = tests;
    }

    if ((act & A_TEST) && tests) {
      ev.tests--;
      if (!ev.tests) {
	act |= A_ENDTEST | A_ENDSTEP | A_DONE;
	rc = PGEN_DONE;
      }
    }

    /* --- Report an event if so directed --- */

    if ((act & A_EVENT) && event && event(rc, &ev, ectx) == PGEN_ABORT) {
      rc = PGEN_ABORT;
      if (!(act & A_DONE)) {
	act |= A_ENDSTEP | A_DONE;
	if (p == P_TEST && rq != PGEN_BEGIN)
	  act |= A_ENDTEST;
      }
    }

    /* --- Close down tester and stepper functions --- */

    if (act & A_ENDTEST)
      test(PGEN_DONE, &ev, tctx);
    if (act & A_ENDSTEP)
      step(PGEN_DONE, &ev, sctx);

    /* --- Stop the entire test if necessary --- */

    if (act & A_DONE)
      break;
  }

  /* --- Tidy up and return --- */

  if (rc == PGEN_ABORT) {
    mp_drop(ev.m);
    ev.m = 0;
  }
  ev.r->ops->destroy(ev.r);
  mp_drop(d);

  return (ev.m);
}

/* --- @pgen_primep@ --- *
 *
 * Arguments:	@mp *p@ = a number to check
 *		@grand *gr@ = a random number source
 *
 * Returns:	Nonzero if @p@ is really prime.
 *
 * Use:		Checks the primality of @p@.  If @p@ is prime, then this
 *		function returns nonzero; if @p@ is really composite then it
 *		%%\emph{probably}%% returns zero, but might not.
 *
 *		Currently, this function uses the Baillie--PSW test, which
 *		combines a single Miller--Rabin test with witness 2 with a
 *		single Frobenius test with parameters chosen using
 *		Selfridge's `Method A'.  No composites are known which pass
 *		this test, though it's conjectured that infinitely many
 *		exist.
 */

int pgen_primep(mp *p, grand *gr)
{
  rabin r;
  int rc;

  if (MP_NEGP(p)) return (0);
  switch (pfilt_smallfactor(p)) {
    case PGEN_DONE: return (1);
    case PGEN_FAIL: return (0);
  }
  rabin_create(&r, p); rc = rabin_test(&r, MP_TWO); rabin_destroy(&r);
  if (rc == PGEN_FAIL) return (0);
  rc = pgen_granfrob(p, 0, 0); if (rc == PGEN_FAIL) return (0);
  return (1);
}

/*----- Test rig ----------------------------------------------------------*/

#ifdef TEST_RIG

#include <mLib/testrig.h>

static int t_primep(dstr *v)
{
  mp *m = *(mp **)v[0].buf;
  int e = *(int *)v[1].buf;
  int r;
  grand *rng;
  int ok = 1;

  rng = fibrand_create(0);
  r = pgen_primep(m, rng);
  GR_DESTROY(rng);
  if (e != r) {
    fputs("\n*** primep failed", stderr);
    fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
    fprintf(stderr, "\nexpected %d", e);
    fprintf(stderr, "\nreported %d", r);
    fputc('\n', stderr);
    ok = 0;
  }

  mp_drop(m);
  assert(mparena_count(MPARENA_GLOBAL) == 0);
  return (ok);
}

static int verify(dstr *v)
{
  mp *m = *(mp **)v[0].buf;
  mp *q = *(mp **)v[1].buf;
  mp *p;
  int ok = 1;

  pgen_filterctx pf;
  rabin r;

  pf.step = 2;
  p = pgen("p", MP_NEW, m, pgen_evspin, 0, 0, pgen_filter, &pf,
	   rabin_iters(mp_bits(m)), pgen_test, &r);
  if (!p || !MP_EQ(p, q)) {
    fputs("\n*** pgen failed", stderr);
    fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
    fputs("\np = ", stderr); mp_writefile(p, stderr, 10);
    fputs("\nq = ", stderr); mp_writefile(q, stderr, 10);
    fputc('\n', stderr);
    ok = 0;
  }

  mp_drop(m);
  mp_drop(q);
  mp_drop(p);
  assert(mparena_count(MPARENA_GLOBAL) == 0);
  return (ok);
}

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

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

/*----- That's all, folks -------------------------------------------------*/
Commit	Line	Data
0f5ec153	1	/* --c--
0f5ec153	2	*
581c854e	3	* Prime generation glue
0f5ec153	4	*
	5	* (c) 1999 Straylight/Edgeware
	6	*/
	7
45c0fd36	8	/----- Licensing notice --------------------------------------------------
0f5ec153	9	*
	10	* This file is part of Catacomb.
	11	*
	12	* Catacomb is free software; you can redistribute it and/or modify
	13	* it under the terms of the GNU Library General Public License as
	14	* published by the Free Software Foundation; either version 2 of the
	15	* License, or (at your option) any later version.
45c0fd36	16	*
0f5ec153	17	* Catacomb is distributed in the hope that it will be useful,
	18	* but WITHOUT ANY WARRANTY; without even the implied warranty of
	19	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	20	* GNU Library General Public License for more details.
45c0fd36	21	*
0f5ec153	22	* You should have received a copy of the GNU Library General Public
	23	* License along with Catacomb; if not, write to the Free
	24	* Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
	25	* MA 02111-1307, USA.
	26	*/
	27
0f5ec153	28	/----- Header files ------------------------------------------------------/
0f5ec153	29
581c854e	30	#include <assert.h>
	31	#include <stdio.h>
	32	#include <stdlib.h>
	33	#include <string.h>
	34
	35	#include "fibrand.h"
	36	#include "grand.h"
0f5ec153	37	#include "mp.h"
581c854e	38	#include "mprand.h"
0f5ec153	39	#include "pgen.h"
581c854e	40	#include "pfilt.h"
0f5ec153	41	#include "rabin.h"
0f5ec153	42
581c854e	43	/----- Standard prime filter ---------------------------------------------/
0f5ec153	44
581c854e	45	/* --- @pgen_filter@ --- */
0f5ec153	46
581c854e	47	int pgen_filter(int rq, pgen_event ev, void p)
0f5ec153	48	{
581c854e	49	pgen_filterctx *f = p;
57541e97	50	int rc = PGEN_FAIL;
581c854e	51
	52	switch (rq) {
	53	case PGEN_BEGIN:
	54	rc = pfilt_create(&f->f, ev->m);
	55	mp_drop(ev->m);
	56	break;
	57	case PGEN_TRY:
	58	mp_drop(ev->m);
581c854e	59	break;
	60	case PGEN_DONE:
	61	pfilt_destroy(&f->f);
	62	return (PGEN_DONE);
57541e97	63	default:
	64	rc = PGEN_ABORT;
	65	break;
0f5ec153	66	}
45c0fd36	67
57541e97	68	if (rc == PGEN_FAIL && !((f->step \| f->f.m->v[0]) & 1))
57541e97	69	rc = pfilt_step(&f->f, 1);
581c854e	70	while (rc == PGEN_FAIL)
	71	rc = pfilt_step(&f->f, f->step);
	72	ev->m = MP_COPY(f->f.m);
0f5ec153	73	return (rc);
	74	}
	75
581c854e	76	/* --- @pgen_jump@ --- *
0f5ec153	77	*
581c854e	78	* Similar to the standard @pgen_filter@, but jumps in large steps rather
581c854e	79	* than small ones.
0f5ec153	80	*/
0f5ec153	81
581c854e	82	int pgen_jump(int rq, pgen_event ev, void p)
0f5ec153	83	{
581c854e	84	pgen_jumpctx *f = p;
581c854e	85	int rc = PGEN_ABORT;
0f5ec153	86
581c854e	87	switch (rq) {
8b021c3f	88	case PGEN_BEGIN: {
	89	mp *g = MP_NEW;
	90	mp_gcd(&g, 0, 0, ev->m, f->j->m);
	91	if (MP_CMP(g, >, MP_ONE)) {
	92	mp_drop(g);
	93	return (PGEN_ABORT);
	94	}
	95	mp_drop(g);
581c854e	96	rc = pfilt_create(&f->f, ev->m);
581c854e	97	mp_drop(ev->m);
8b021c3f	98	} break;
581c854e	99	case PGEN_TRY:
	100	mp_drop(ev->m);
	101	rc = pfilt_jump(&f->f, f->j);
	102	break;
	103	case PGEN_DONE:
	104	pfilt_destroy(&f->f);
	105	return (PGEN_DONE);
	106	}
45c0fd36	107
581c854e	108	while (rc == PGEN_FAIL)
	109	rc = pfilt_jump(&f->f, f->j);
	110	ev->m = MP_COPY(f->f.m);
	111	return (rc);
	112	}
0f5ec153	113
581c854e	114	/----- Standard prime test -----------------------------------------------/
0f5ec153	115
581c854e	116	/* --- @pgen_test@ --- */
0f5ec153	117
581c854e	118	int pgen_test(int rq, pgen_event ev, void p)
	119	{
	120	rabin *r = p;
a938be51	121	mp *a = MP_NEW;
581c854e	122	int rc = PGEN_ABORT;
0f5ec153	123
581c854e	124	switch (rq) {
	125	case PGEN_BEGIN:
	126	rabin_create(r, ev->m);
	127	rc = PGEN_TRY;
	128	break;
283b9af0	129	case PGEN_TRY:
a938be51 MW	130	a = mprand_range(a, ev->m, ev->r, 0);
a938be51 MW	131	rc = rabin_rtest(r, a);
283b9af0	132	break;
581c854e	133	case PGEN_DONE:
	134	rabin_destroy(r);
	135	rc = PGEN_DONE;
	136	break;
0f5ec153	137	}
0f5ec153	138
a938be51	139	mp_drop(a);
0f5ec153	140	return (rc);
	141	}
	142
a889a2b3 MW	143	/* --- @pgen_bailliepswtest@ --- */
	144
	145	int pgen_bailliepswtest(int rq, pgen_event ev, void p)
	146	{
	147	rabin r;
	148	int rc;
	149
	150	switch (rq) {
	151	case PGEN_BEGIN:
	152	if (ev->tests != 2) rc = PGEN_ABORT;
	153	else rc = PGEN_TRY;
	154	break;
	155
	156	case PGEN_DONE:
	157	rc = PGEN_DONE;
	158	break;
	159
	160	case PGEN_TRY:
	161	switch (ev->tests) {
	162	case 2:
	163	rabin_create(&r, ev->m);
	164	rc = rabin_test(&r, MP_TWO);
	165	rabin_destroy(&r);
	166	break;
	167	case 1:
	168	rc = pgen_granfrob(ev->m, 0, 0);
	169	break;
	170	default:
	171	rc = PGEN_ABORT;
	172	break;
	173	}
	174	break;
	175
	176	default:
	177	rc = PGEN_ABORT;
	178	break;
	179	}
	180
	181	return (rc);
	182	}
	183
581c854e	184	/----- The main driver ---------------------------------------------------/
	185
	186	/* --- @pgen@ --- *
bd98b2df	187	*
581c854e	188	* Arguments: @const char *name@ = name of the value being searched for
	189	* @mp *d@ = destination for the result integer
	190	* @mp *m@ = start value to pass to stepper
	191	* @pgen_proc *event@ = event handler function
	192	* @void *ectx@ = context argument for event andler
	193	* @unsigned steps@ = number of steps to take in search
	194	* @pgen_proc *step@ = stepper function to use
	195	* @void *sctx@ = context argument for stepper
	196	* @unsigned tests@ = number of tests to make
	197	* @pgen_proc *test@ = tester function to use
	198	* @void *tctx@ = context argument for tester
bd98b2df	199	*
581c854e	200	* Returns: Pointer to final result, or null.
bd98b2df	201	*
581c854e	202	* Use: A generalized prime-number search skeleton. Yes, that's a
581c854e	203	* scary number of arguments.
bd98b2df	204	*/
bd98b2df	205
581c854e	206	mp pgen(const char name, mp d, mp m, pgen_proc event, void ectx,
	207	unsigned steps, pgen_proc step, void sctx,
	208	unsigned tests, pgen_proc test, void tctx)
bd98b2df	209	{
581c854e	210	pgen_event ev;
	211	int rq, rc;
	212	pgen_proc *proc;
	213	void *ctx;
4108c8d2	214	int p;
	215
	216	enum { P_STEP, P_TEST };
bd98b2df	217
581c854e	218	/* --- Set up the initial event block --- */
bd98b2df	219
581c854e	220	ev.name = name;
	221	if (m)
	222	ev.m = MP_COPY(m);
	223	else
	224	ev.m = 0;
1222a44a MW	225	ev.steps = steps;
1222a44a MW	226	ev.tests = tests;
581c854e	227	ev.r = fibrand_create(0);
bd98b2df	228
581c854e	229	/* --- Tell the event handler we're under way --- */
bd98b2df	230
4108c8d2	231	if (event && event(PGEN_BEGIN, &ev, ectx) == PGEN_ABORT) {
4108c8d2	232	ev.r->ops->destroy(ev.r);
581c854e	233	return (0);
4108c8d2	234	}
bd98b2df	235
581c854e	236	/* --- Set up for the initial call --- */
bd98b2df	237
4108c8d2	238	proc = step; ctx = sctx; p = P_STEP; rq = PGEN_BEGIN;
0f5ec153	239
581c854e	240	/* --- Enter the great maelstrom of state transitions --- */
0f5ec153	241
581c854e	242	for (;;) {
	243	unsigned act = 0;
	244
4108c8d2	245	#define A_STEP 1u
	246	#define A_TEST 2u
	247	#define A_EVENT 4u
	248	#define A_ENDTEST 8u
	249	#define A_ENDSTEP 16u
	250	#define A_DONE 32u
581c854e	251
	252	/* --- Call the procedure and decide what to do next --- */
	253
	254	rc = proc(rq, &ev, ctx);
	255	switch (rc) {
	256	case PGEN_TRY:
4108c8d2	257	if (p == P_TEST)
581c854e	258	rq = PGEN_TRY;
	259	else {
	260	act \|= A_EVENT;
4108c8d2	261	proc = test; ctx = tctx; p = P_TEST;
581c854e	262	rq = PGEN_BEGIN;
	263	}
	264	break;
	265	case PGEN_PASS:
	266	act \|= A_TEST \| A_EVENT;
4108c8d2	267	if (p == P_TEST)
581c854e	268	rq = PGEN_TRY;
581c854e	269	else {
4108c8d2	270	proc = test; ctx = tctx; p = P_TEST;
581c854e	271	rq = PGEN_BEGIN;
	272	}
	273	break;
	274	case PGEN_FAIL:
	275	act \|= A_STEP;
4108c8d2	276	if (p == P_TEST) {
581c854e	277	act \|= A_ENDTEST \| A_EVENT;
4108c8d2	278	proc = step; ctx = sctx; p = P_STEP;
581c854e	279	}
	280	rq = PGEN_TRY;
	281	break;
	282	case PGEN_DONE:
	283	act \|= A_EVENT \| A_DONE \| A_ENDSTEP;
4108c8d2	284	if (p == P_TEST)
581c854e	285	act \|= A_ENDTEST;
	286	break;
	287	case PGEN_ABORT:
	288	act \|= A_EVENT \| A_DONE;
4108c8d2	289	if (p == P_TEST \|\| rq == PGEN_TRY)
581c854e	290	act \|= A_ENDSTEP;
4108c8d2	291	if (p == P_TEST && rq != PGEN_BEGIN)
581c854e	292	act \|= A_ENDTEST;
	293	break;
	294	default:
	295	assert(((void)"Invalid response from function", 0));
	296	break;
	297	}
0f5ec153	298
581c854e	299	/* --- If decrementing counters is requested, do that --- */
0f5ec153	300
581c854e	301	if ((act & A_STEP) && steps) {
1222a44a MW	302	ev.steps--;
1222a44a MW	303	if (!ev.steps) {
581c854e	304	act \|= A_EVENT \| A_ENDSTEP \| A_DONE;
	305	rc = PGEN_ABORT;
	306	}
1222a44a	307	ev.tests = tests;
581c854e	308	}
0f5ec153	309
581c854e	310	if ((act & A_TEST) && tests) {
1222a44a MW	311	ev.tests--;
1222a44a MW	312	if (!ev.tests) {
581c854e	313	act \|= A_ENDTEST \| A_ENDSTEP \| A_DONE;
	314	rc = PGEN_DONE;
	315	}
0f5ec153	316	}
0f5ec153	317
581c854e	318	/* --- Report an event if so directed --- */
bd98b2df	319
581c854e	320	if ((act & A_EVENT) && event && event(rc, &ev, ectx) == PGEN_ABORT) {
	321	rc = PGEN_ABORT;
	322	if (!(act & A_DONE)) {
	323	act \|= A_ENDSTEP \| A_DONE;
8501f5f0	324	if (p == P_TEST && rq != PGEN_BEGIN)
581c854e	325	act \|= A_ENDTEST;
581c854e	326	}
45c0fd36	327	}
bd98b2df	328
581c854e	329	/* --- Close down tester and stepper functions --- */
0f5ec153	330
581c854e	331	if (act & A_ENDTEST)
	332	test(PGEN_DONE, &ev, tctx);
	333	if (act & A_ENDSTEP)
	334	step(PGEN_DONE, &ev, sctx);
	335
	336	/* --- Stop the entire test if necessary --- */
	337
	338	if (act & A_DONE)
	339	break;
	340	}
	341
	342	/* --- Tidy up and return --- */
	343
	344	if (rc == PGEN_ABORT) {
	345	mp_drop(ev.m);
	346	ev.m = 0;
	347	}
	348	ev.r->ops->destroy(ev.r);
eab06f16	349	mp_drop(d);
581c854e	350
581c854e	351	return (ev.m);
0f5ec153	352	}
0f5ec153	353
34e4f738	354	/* --- @pgen_primep@ --- *
	355	*
	356	* Arguments: @mp *p@ = a number to check
	357	* @grand *gr@ = a random number source
	358	*
	359	* Returns: Nonzero if @p@ is really prime.
6001a9ff MW	360	*
	361	* Use: Checks the primality of @p@. If @p@ is prime, then this
	362	* function returns nonzero; if @p@ is really composite then it
	363	* %%\emph{probably}%% returns zero, but might not.
	364	*
	365	* Currently, this function uses the Baillie--PSW test, which
	366	* combines a single Miller--Rabin test with witness 2 with a
	367	* single Frobenius test with parameters chosen using
	368	* Selfridge's `Method A'. No composites are known which pass
	369	* this test, though it's conjectured that infinitely many
	370	* exist.
34e4f738	371	*/
	372
	373	int pgen_primep(mp p, grand gr)
	374	{
34e4f738	375	rabin r;
6001a9ff	376	int rc;
34e4f738	377
a69a3efd	378	if (MP_NEGP(p)) return (0);
34e4f738	379	switch (pfilt_smallfactor(p)) {
	380	case PGEN_DONE: return (1);
	381	case PGEN_FAIL: return (0);
	382	}
6001a9ff MW	383	rabin_create(&r, p); rc = rabin_test(&r, MP_TWO); rabin_destroy(&r);
	384	if (rc == PGEN_FAIL) return (0);
	385	rc = pgen_granfrob(p, 0, 0); if (rc == PGEN_FAIL) return (0);
	386	return (1);
34e4f738	387	}
34e4f738	388
581c854e	389	/----- Test rig ----------------------------------------------------------/
0f5ec153	390
	391	#ifdef TEST_RIG
	392
	393	#include <mLib/testrig.h>
	394
d31f4a79 MW	395	static int t_primep(dstr *v)
	396	{
	397	mp m = (mp **)v[0].buf;
	398	int e = (int )v[1].buf;
	399	int r;
	400	grand *rng;
	401	int ok = 1;
	402
	403	rng = fibrand_create(0);
	404	r = pgen_primep(m, rng);
	405	GR_DESTROY(rng);
	406	if (e != r) {
	407	fputs("\n*** primep failed", stderr);
	408	fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
	409	fprintf(stderr, "\nexpected %d", e);
	410	fprintf(stderr, "\nreported %d", r);
	411	fputc('\n', stderr);
	412	ok = 0;
	413	}
	414
	415	mp_drop(m);
	416	assert(mparena_count(MPARENA_GLOBAL) == 0);
45c0fd36	417	return (ok);
d31f4a79 MW	418	}
d31f4a79 MW	419
0f5ec153	420	static int verify(dstr *v)
	421	{
	422	mp m = (mp **)v[0].buf;
581c854e	423	mp q = (mp **)v[1].buf;
581c854e	424	mp *p;
0f5ec153	425	int ok = 1;
0f5ec153	426
581c854e	427	pgen_filterctx pf;
581c854e	428	rabin r;
0f5ec153	429
581c854e	430	pf.step = 2;
	431	p = pgen("p", MP_NEW, m, pgen_evspin, 0, 0, pgen_filter, &pf,
	432	rabin_iters(mp_bits(m)), pgen_test, &r);
22bab86c	433	if (!p \|\| !MP_EQ(p, q)) {
581c854e	434	fputs("\n*** pgen failed", stderr);
0f5ec153	435	fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
581c854e	436	fputs("\np = ", stderr); mp_writefile(p, stderr, 10);
581c854e	437	fputs("\nq = ", stderr); mp_writefile(q, stderr, 10);
0f5ec153	438	fputc('\n', stderr);
	439	ok = 0;
	440	}
	441
	442	mp_drop(m);
581c854e	443	mp_drop(q);
eab06f16	444	mp_drop(p);
d94f85ac	445	assert(mparena_count(MPARENA_GLOBAL) == 0);
0f5ec153	446	return (ok);
	447	}
	448
	449	static test_chunk tests[] = {
	450	{ "pgen", verify, { &type_mp, &type_mp, 0 } },
d31f4a79	451	{ "primep", t_primep, { &type_mp, &type_int, 0 } },
0f5ec153	452	{ 0, 0, { 0 } }
	453	};
	454
	455	int main(int argc, char *argv[])
	456	{
	457	sub_init();
0f00dc4c	458	test_run(argc, argv, tests, SRCDIR "/t/pgen");
0f5ec153	459	return (0);
0f5ec153	460	}
0f5ec153	461	#endif
	462
	463	/----- That's all, folks -------------------------------------------------/