(pfilt_smallfactor): New function for doing trial division the hard

way.

(pfilt_create): Use @mpx_udivn@ for computing residues, for improved
performance.

Pull the `small prime' test into a separate function, and do it
properly.

diff --git a/pfilt.c b/pfilt.c
index eda5116..4b58304 100644 (file)
--- a/pfilt.c
+++ b/pfilt.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: pfilt.c,v 1.2 2000/06/17 11:54:27 mdw Exp $
+ * $Id: pfilt.c,v 1.3 2000/08/15 21:44:27 mdw Exp $
  *
  * Finding and testing prime numbers
  *
@@ -30,6 +30,16 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: pfilt.c,v $
+ * Revision 1.3  2000/08/15 21:44:27  mdw
+ * (pfilt_smallfactor): New function for doing trial division the hard
+ * way.
+ *
+ * (pfilt_create): Use @mpx_udivn@ for computing residues, for improved
+ * performance.
+ *
+ * Pull the `small prime' test into a separate function, and do it
+ * properly.
+ *
  * Revision 1.2  2000/06/17 11:54:27  mdw
  * Use new MP memory management functions.
  *
@@ -50,13 +60,81 @@
 /*----- Header files ------------------------------------------------------*/
 
 #include "mp.h"
-#include "mpmont.h"
+#include "mpint.h"
 #include "pfilt.h"
 #include "pgen.h"
 #include "primetab.h"
+#include "primorial.h"
 
 /*----- Main code ---------------------------------------------------------*/
 
+/* --- @smallenough@ --- *
+ *
+ * Arguments:  @mp *m@ = integer to test
+ *
+ * Returns:    One of the @PGEN@ result codes.
+ *
+ * Use:                Assuming that @m@ has been tested by trial division on every
+ *             prime in the small-primes array, this function will return
+ *             @PGEN_DONE@ if the number is less than the square of the
+ *             largest small prime.
+ */
+
+static int smallenough(mp *m)
+{
+  static mp *max = 0;
+  int rc = PGEN_TRY;
+
+  if (!max) {
+    max = mp_fromuint(MP_NEW, MAXPRIME);
+    max = mp_sqr(max, max);
+    max->a->n--; /* Permanent allocation */
+  }
+  if (MP_CMP(m, <, max))
+    rc = PGEN_DONE;
+  return (rc);
+}
+
+/* --- @pfilt_smallfactor@ --- *
+ *
+ * Arguments:  @mp *m@ = integer to test
+ *
+ * Returns:    One of the @PGEN@ result codes.
+ *
+ * Use:                Tests a number by dividing by a number of small primes.  This
+ *             is a useful first step if you're testing random primes; for
+ *             sequential searches, @pfilt_create@ works better.
+ */
+
+int pfilt_smallfactor(mp *m)
+{
+  int rc = PGEN_TRY;
+  int i;
+  size_t sz = MP_LEN(m);
+  mpw *v = mpalloc(m->a, sz);
+
+  /* --- Fill in the residues --- */
+
+  for (i = 0; i < NPRIME; i++) {
+    if (!mpx_udivn(v, v + sz, m->v, m->vl, primetab[i])) {
+      if (MP_LEN(m) == 1 && m->v[0] == primetab[i])
+       rc = PGEN_DONE;
+      else
+       rc = PGEN_FAIL;
+    }
+  }
+
+  /* --- Check for small primes --- */
+
+  if (rc == PGEN_TRY)
+    rc = smallenough(m);
+
+  /* --- Done --- */
+
+  mpfree(m->a, v);
+  return (rc);
+}
+
 /* --- @pfilt_create@ --- *
  *
  * Arguments:  @pfilt *p@ = pointer to prime filtering context
@@ -74,9 +152,8 @@ int pfilt_create(pfilt *p, mp *m)
 {
   int rc = PGEN_TRY;
   int i;
-  mp *r = MP_NEW;
-  mpw qw;
-  mp q;
+  size_t sz = MP_LEN(m);
+  mpw *v = mpalloc(m->a, sz);
 
   /* --- Take a copy of the number --- */
 
@@ -85,11 +162,8 @@ int pfilt_create(pfilt *p, mp *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];
+    p->r[i] = mpx_udivn(v, v + sz, m->v, m->vl, primetab[i]);
     if (!p->r[i] && rc == PGEN_TRY) {
       if (MP_LEN(m) == 1 && m->v[0] == primetab[i])
        rc = PGEN_DONE;
@@ -98,9 +172,14 @@ int pfilt_create(pfilt *p, mp *m)
     }
   }
 
+  /* --- Check for small primes --- */
+
+  if (rc == PGEN_TRY)
+    rc = smallenough(m);
+
   /* --- Done --- */
 
-  mp_drop(r);
+  mpfree(m->a, v);
   return (rc);
 }
 
@@ -155,11 +234,10 @@ int pfilt_step(pfilt *p, mpw step)
     }
   }
 
-  /* --- Small numbers must be prime --- */
+  /* --- Check for small primes --- */
 
-  if (rc == PGEN_TRY && MP_LEN(p->m) == 1 &&
-      p->m->v[0] < MAXPRIME * MAXPRIME)
-    rc = PGEN_DONE;
+  if (rc == PGEN_TRY)
+    rc = smallenough(p->m);
 
   /* --- Done --- */
 
@@ -213,11 +291,10 @@ int pfilt_muladd(pfilt *p, const pfilt *q, mpw m, mpw a)
     }
   }
 
-  /* --- Small numbers must be prime --- */
+  /* --- Check for small primes --- */
 
-  if (rc == PGEN_TRY && MP_LEN(p->m) == 1 &&
-      p->m->v[0] < MAXPRIME * MAXPRIME)
-    rc = PGEN_DONE;
+  if (rc == PGEN_TRY)
+    rc = smallenough(p->m);
 
   /* --- Finished --- */
 
@@ -265,11 +342,10 @@ int pfilt_jump(pfilt *p, const pfilt *j)
     }
   }
 
-  /* --- Small numbers must be prime --- */
+  /* --- Check for small primes --- */
 
-  if (rc == PGEN_TRY && MP_LEN(p->m) == 1 &&
-      p->m->v[0] < MAXPRIME * MAXPRIME)
-    rc = PGEN_DONE;
+  if (rc == PGEN_TRY)
+    rc = smallenough(p->m);
 
   /* --- Done --- */