add cfrac benchmark

1 files changed, 731 insertions, 0 deletions

diff --git a/src/benchmarks/cfrac/pcfrac.c b/src/benchmarks/cfrac/pcfrac.c
new file mode 100644
index 0000000..2a0d032
--- /dev/null
+++ b/src/benchmarks/cfrac/pcfrac.c
@@ -0,0 +1,731 @@
+/*
+ * pcfrac: Implementation of the continued fraction factoring algoritm
+ *
+ * Every two digits additional appears to double the factoring time
+ * 
+ * Written by Dave Barrett (barrett%asgard@boulder.Colorado.EDU)
+ */
+#include <string.h>
+#include <stdio.h>
+#include <math.h>
+
+#ifdef __STDC__
+#include <stdlib.h>
+#endif
+#include "precision.h"
+#include "pfactor.h"
+
+extern int verbose;
+
+unsigned cfracNabort  = 0;
+unsigned cfracTsolns  = 0;
+unsigned cfracPsolns  = 0;
+unsigned cfracT2solns = 0;
+unsigned cfracFsolns  = 0;
+
+extern unsigned short primes[];
+extern unsigned primesize;
+
+typedef unsigned *uptr;
+typedef uptr	 uvec;
+typedef unsigned char *solnvec;
+typedef unsigned char *BitVector;
+
+typedef struct SolnStruc {
+   struct SolnStruc *next;
+   precision x;		/* lhs of solution */
+   precision t;		/* last large prime remaining after factoring */
+   precision r;		/* accumulated root of pm for powers >= 2 */
+   BitVector e;		/* bit vector of factorbase powers mod 2 */
+} Soln;
+
+typedef Soln *SolnPtr;
+
+#define BPI(x)	((sizeof x[0]) << 3)
+
+void setBit(bv, bno, value)
+   register BitVector bv;
+   register unsigned bno, value;
+{
+   bv  += bno / BPI(bv);
+   bno %= BPI(bv);
+   *bv |= ((value != 0) << bno);
+}
+
+unsigned getBit(bv, bno)
+   register BitVector bv;
+   register unsigned bno;
+{
+   register unsigned res;
+
+   bv  += bno / BPI(bv);
+   bno %= BPI(bv);
+   res  = (*bv >> bno) & 1;
+
+   return res;
+}
+
+BitVector newBitVector(value, size)
+   register solnvec value;
+   unsigned size;
+{
+   register BitVector res;
+   register solnvec  vp = value + size;
+   unsigned msize = ((size + BPI(res)-1) / BPI(res)) * sizeof res[0];
+
+#ifdef BWGC 		 
+   res = (BitVector) gc_malloc(msize);
+#else		 
+   res = (BitVector) malloc(msize);
+#endif
+   if (res == (BitVector) 0) return res;
+
+   memset(res, '\0', msize);
+   do {
+      if (*--vp) {
+	 setBit(res, vp - value, (unsigned) *vp);
+      }
+   } while (vp != value);
+   return res;
+}
+
+void printSoln(stream, prefix, suffix, pm, m, p, t, e)
+   FILE *stream;
+   char *prefix, *suffix;
+   register unsigned *pm, m;
+   precision p, t;
+   register solnvec e;
+{
+   register unsigned i, j = 0;
+
+   for (i = 1; i <= m; i++) j += (e[i] != 0);
+
+   fputs(prefix, stream); 
+   fputp(stream, pparm(p)); fputs(" = ", stream);
+   if (*e & 1) putc('-', stream);  else putc('+', stream);
+   fputp(stream, pparm(t)); 
+
+   if (j >= 1) fputs(" *", stream);
+   do {
+      e++;
+      switch (*e) {
+      case 0: break;
+      case 1: fprintf(stream, " %u", *pm); break;
+      default:
+	 fprintf(stream, " %u^%u", *pm, (unsigned) *e);
+      }
+      pm++;
+   } while (--m);
+
+   fputs(suffix, stream); 
+   fflush(stream);
+   pdestroy(p); pdestroy(t);
+}
+
+/*
+ * Combine two solutions 
+ */
+void combineSoln(x, t, e, pm, m, n, bp)
+   precision *x, *t, n;
+   uvec	     pm;
+   register  solnvec e;
+   unsigned  m;
+   SolnPtr   bp;
+{
+   register unsigned j;
+
+   (void) pparm(n);
+   if (bp != (SolnPtr) 0) {
+      pset(x, pmod(pmul(bp->x, *x), n));
+      pset(t, pmod(pmul(bp->t, *t), n));
+      pset(t, pmod(pmul(bp->r, *t), n));
+      e[0] += getBit(bp->e, 0);
+   }
+   e[0] &= 1;
+   for (j = 1; j <= m; j++) {
+      if (bp != (SolnPtr) 0) e[j] += getBit(bp->e, j);
+      if (e[j] > 2) {
+	 pset(t, pmod(pmul(*t, 
+	    ppowmod(utop(pm[j-1]), utop((unsigned) e[j]>>1), n)), n));
+	 e[j] &= 1;
+      } else if (e[j] == 2) {
+	 pset(t, pmod(pmul(*t, utop(pm[j-1])), n));
+	 e[j] = 0;
+      }
+   }
+   pdestroy(n);
+}
+
+/*
+ * Create a normalized solution structure from the given inputs
+ */
+SolnPtr newSoln(n, pm, m, next, x, t, e)
+   precision n;
+   unsigned  m;
+   uvec     pm;
+   SolnPtr next;
+   precision x, t;
+   solnvec e;
+{
+#ifdef BWGC 		 
+   SolnPtr bp = (SolnPtr) gc_malloc(sizeof (Soln));
+#else		 
+   SolnPtr bp = (SolnPtr) malloc(sizeof (Soln));
+#endif
+
+   if (bp != (SolnPtr) 0) {
+      bp->next = next;
+      bp->x = pnew(x);
+      bp->t = pnew(t);
+      bp->r = pnew(pone);
+      /*
+       * normalize e, put the result in bp->r and e
+       */
+      combineSoln(&bp->x, &bp->r, e, pm, m, pparm(n), (SolnPtr) 0);
+      bp->e = newBitVector(e, m+1);			/* BitVector */
+   }
+
+   pdestroy(n);
+   return bp;
+}
+
+void freeSoln(p)
+   register SolnPtr p;
+{
+   if (p != (SolnPtr) 0) {
+      pdestroy(p->x);
+      pdestroy(p->t);
+      pdestroy(p->r);
+#ifndef IGNOREFREE
+      free(p->e);			/* BitVector */
+      free(p);
+#endif
+   }
+}
+
+void freeSolns(p)
+   register SolnPtr p;
+{
+   register SolnPtr l;
+   
+   while (p != (SolnPtr) 0) {
+      l = p;
+      p = p->next;
+      freeSoln(l);
+   }
+}
+
+SolnPtr findSoln(sp, t)
+   register SolnPtr sp;
+   precision t;
+{
+   (void) pparm(t);
+   while (sp != (SolnPtr) 0) {
+      if peq(sp->t, t) break;
+      sp = sp->next;
+   }
+   pdestroy(t);
+   return sp;
+}
+
+static unsigned pcfrac_k      = 1;
+static unsigned pcfrac_m      = 0;
+static unsigned pcfrac_aborts = 3;
+
+/*
+ * Structure for early-abort.  Last entry must be <(unsigned *) 0, uUndef>
+ */
+typedef struct {
+   unsigned  *pm;	/* bound check occurs before using this pm entry */
+   precision bound;	/* max allowable residual to prevent abort */
+} EasEntry;
+
+typedef EasEntry *EasPtr;
+
+void freeEas(eas)
+   EasPtr eas;
+{
+   register EasPtr ep = eas;
+
+   if (ep != (EasPtr) 0) {
+      while (ep->pm != 0) {
+	 pdestroy(ep->bound);
+	 ep++;
+      }
+#ifndef IGNOREFREE
+      free(eas);
+#endif
+   }
+}
+
+/*
+ * Return Pomerance's L^alpha (L = exp(sqrt(log(n)*log(log(n)))))
+ */
+double pomeranceLpow(n, y) 
+   double n;
+   double y;
+{
+   double lnN = log(n);
+   double res = exp(y * sqrt(lnN * log(lnN)));
+   return res;
+}
+
+/*
+ * Pomerance's value 'a' from page 122 "of Computational methods in Number
+ * Theory", part 1, 1982.
+ */
+double cfracA(n, aborts)
+   double n;
+   unsigned aborts;
+{
+   return 1.0 / sqrt(6.0 + 2.0 / ((double) aborts + 1.0));
+}
+
+/*
+ * Returns 1 if a is a quadratic residue of odd prime p, 
+ * p-1 if non-quadratic residue, 0 otherwise (gcd(a,p)<>1)
+ */
+#define plegendre(a,p) ppowmod(a, phalf(psub(p, pone)), p)
+
+/*
+ * Create a table of small primes of quadratic residues of n
+ *
+ * Input:   
+ *    n      - the number to be factored
+ *    k      - the multiple of n to be factored
+ *    *m     - the number of primes to generate (0 to select best)
+ *    aborts - the number of early aborts
+ *
+ * Assumes that plegendre # 0, for if it is, that pm is a factor of n.
+ * This algorithm already assumes you've used trial division to eliminate
+ * all of these!
+ *
+ * Returns: the list of primes actually generated (or (unsigned *) 0 if nomem)
+ *          *m changed to reflect the number of elements in the list
+ */
+uvec pfactorbase(n, k, m, aborts)
+   precision n;
+   unsigned k;
+   unsigned *m, aborts;
+{
+   double   dn, a;
+   register unsigned short *primePtr = primes;
+   register unsigned count = *m;
+   unsigned maxpm = primes[primesize-1];
+   unsigned *res  = (uvec) 0, *pm;
+   precision nk   = pnew(pmul(pparm(n), utop(k)));
+
+   if (*m == 0) {			/* compute a suitable m */
+      dn    = ptod(nk);
+      a     = cfracA(dn, aborts);
+      maxpm = (unsigned) (pomeranceLpow(dn, a) + 0.5);
+      do {
+	 if ((unsigned) *primePtr++ >= maxpm) break;
+      } while ((unsigned) *primePtr != 1);
+      count = primePtr - primes;
+      primePtr = primes;
+   }
+   /*
+    * This m tends to be too small for small n, and becomes closer to 
+    * optimal as n goes to infinity.  For 30 digits, best m is ~1.5 this m.
+    * For 38 digits, best m appears to be ~1.15 this m. It's appears to be
+    * better to guess too big than too small.
+    */
+#ifdef BWGC 		 
+   res = (uvec) gc_malloc(count * sizeof (unsigned));
+#else		 
+   res = (uvec) malloc(count * sizeof (unsigned));
+#endif
+   if (res == (uvec) 0) goto doneMk;
+
+   pm    = res;
+   *pm++ = (unsigned) *primePtr++;		/* two is first element */
+   count = 1;
+   if (count != *m) do {
+      if (picmp(plegendre(nk, utop((unsigned) *primePtr)), 1) <= 0) { /* 0,1 */
+	 *pm++ = *primePtr;
+	 count++;
+	 if (count == *m) break;
+	 if ((unsigned) *primePtr >= maxpm) break;
+      }
+      ++primePtr;
+   } while (*primePtr != 1);
+   *m = count;
+
+doneMk:
+   pdestroy(nk);
+   pdestroy(n);
+   return res;
+}
+
+/*
+ * Compute Pomerance's early-abort-stragegy
+ */
+EasPtr getEas(n, k, pm, m, aborts)
+   precision n;
+   unsigned k, *pm, m, aborts;
+{
+   double x    = 1.0 / ((double) aborts + 1.0);
+   double a    = 1.0 / sqrt(6.0 + 2.0 * x);
+   double ax   = a * x, csum = 1.0, tia = 0.0;
+   double dn, dpval, dbound, ci;
+   unsigned i, j, pval;
+
+   precision bound = pUndef;
+   EasPtr   eas;
+   
+   if (aborts == 0) return (EasPtr) 0;
+
+#ifdef BWGC 		 
+   eas = (EasPtr) gc_malloc((aborts+1) * sizeof (EasEntry));
+#else		 
+   eas = (EasPtr) malloc((aborts+1) * sizeof (EasEntry));
+#endif
+   if (eas == (EasPtr) 0) return eas;
+
+   dn   = ptod(pmul(utop(k), pparm(n)));	/* should this be n ? */
+   for (i = 1; i <= aborts; i++) {
+      eas[i-1].pm = (unsigned *) 0;
+      eas[i-1].bound = pUndef;
+      tia        += ax;
+      ci          = 4.0 * tia * tia / (double) i;
+      csum       -= ci;
+      dpval 	  = pomeranceLpow(dn, tia);
+      dbound      = pow(dn, 0.5 * csum);
+
+      pval	  = (unsigned) (dpval + 0.5);
+      pset(&bound, dtop(dbound));
+      for (j = 0; j < m; j++) {
+	 if (pm[j] >= pval) goto foundpm;
+      }
+      break;
+foundpm:
+      if (verbose > 1) {
+	 printf(" Abort %u on p = %u (>=%u) and q > ", i, pm[j], pval);
+	 fputp(stdout, bound);  putc('\n', stdout);
+	 fflush(stdout);
+      }
+      eas[i-1].pm = &pm[j];
+      pset(&eas[i-1].bound, bound);
+   }
+   eas[i-1].pm = (unsigned *) 0;
+   eas[i-1].bound = pUndef;
+
+   pdestroy(bound);
+   pdestroy(n);
+
+   return eas;
+}
+
+/*
+ * Factor the argument Qn using the primes in pm.  Result stored in exponent
+ * vector e, and residual factor, f.  If non-null, eas points to a list of
+ * early-abort boundaries.
+ *
+ * e is set to the number of times each prime in pm divides v.
+ *
+ * Returns: 
+ *    -2 - if factoring aborted because of early abort
+ *    -1 - factoring failed
+ *     0 - if result is a "partial" factoring
+ *     1 - normal return (a "full" factoring)
+ */
+int pfactorQ(f, t, pm, e, m, eas)
+   precision *f;
+   precision t;
+   register unsigned *pm;
+   register solnvec  e;
+   register unsigned m;
+   EasEntry *eas;
+{
+   precision maxp  = pUndef; 
+   unsigned  maxpm = pm[m-1], res = 0;
+   register unsigned *pp = (unsigned *) 0;
+   
+   (void) pparm(t);
+
+   if (eas != (EasEntry *) 0) {
+      pp = eas->pm;
+      pset(&maxp, eas->bound);
+   }
+
+   memset((char *) e, '\0', m * sizeof e[0]);  /* looks slow here, but isn't */
+
+   while (peven(t)) {		/* assume 2 1st in pm; save time */
+      pset(&t, phalf(t));
+      (*e)++;
+   }
+   --m;
+
+   do {
+      e++; pm++;
+      if (pm == pp) {		/* check for early abort */
+	 if (pgt(t, maxp)) {
+	    res = -2;
+	    goto gotSoln;
+	 }
+	 eas++;
+	 pp = eas->pm;
+	 pset(&maxp, eas->bound);
+      }
+      while (pimod(t, (int) *pm) == 0) {
+	 pset(&t, pidiv(t, (int) *pm));
+	 (*e)++;
+      }
+   } while (--m != 0);
+   res = -1;
+   if (picmp(t, 1) == 0) {
+      res = 1;
+   } else if (picmp(pidiv(t, (int) *pm), maxpm) <= 0) {
+#if 0 		/* it'll never happen;  Honest! If so, pm is incorrect. */
+      if (picmp(t, maxpm) <= 0) {
+	 fprintf(stderr, "BUG: partial with t < maxpm! t = "); 
+	 fputp(stderr, t); putc('\n', stderr);
+      }
+#endif
+      res = 0;
+   }
+gotSoln:
+   pset(f, t);
+   pdestroy(t); pdestroy(maxp);
+   return res;
+}
+
+/*
+ * Attempt to factor n using continued fractions (n must NOT be prime)
+ *
+ * n        - The number to attempt to factor
+ * maxCount - if non-null, points to the maximum number of iterations to try.
+ *
+ * This algorithm may fail if it get's into a cycle or maxCount expires
+ * If failed, n is returned.
+ *
+ * This algorithm will loop indefinitiely in n is prime.
+ *
+ * This an implementation of Morrison and Brillhart's algorithm, with 
+ * Pomerance's early abort strategy, and Knuth's method to find best k.
+ */
+precision pcfrac(n, maxCount)
+   precision n;
+   unsigned *maxCount;
+{
+   unsigned k      = pcfrac_k;
+   unsigned m      = pcfrac_m;
+   unsigned aborts = pcfrac_aborts;
+   SolnPtr  oddt   = (SolnPtr) 0, sp, bp, *b;
+   EasPtr   eas    = (EasPtr) 0;
+   uvec     pm     = (uvec) 0; 
+   solnvec  e	   = (solnvec) 0;
+   unsigned bsize, s = 0, count = 0;
+   register unsigned h, j;
+   int i;
+
+   precision t = pUndef,
+      r  = pUndef, twog   = pUndef,  u = pUndef, lastU  = pUndef,
+      Qn = pUndef, lastQn = pUndef, An = pUndef, lastAn = pUndef,
+      x  = pUndef, y      = pUndef, qn = pUndef,     rn = pUndef;
+
+   precision res = pnew(pparm(n));		/* default res is argument */
+   
+   pm = pfactorbase(n, k, &m, aborts); 		/* m may have been reduced */
+
+   bsize  = (m+2) * sizeof (SolnPtr);
+#ifdef BWGC 		 
+   b      = (SolnPtr *) gc_malloc(bsize);
+#else		 
+   b      = (SolnPtr *) malloc(bsize);
+#endif
+   if (b == (SolnPtr *) 0) goto nomem;
+
+#ifdef BWGC 		 
+   e = (solnvec) gc_malloc((m+1) * sizeof e[0]); 
+#else		 
+   e = (solnvec) malloc((m+1) * sizeof e[0]); 
+#endif
+   if (e == (solnvec) 0) {
+nomem:
+      errorp(PNOMEM, "pcfrac", "out of memory");
+      goto bail;
+   }
+
+   memset(b, '\0', bsize);				 /* F1: Initialize */
+   if (maxCount != (unsigned *) 0) count = *maxCount;
+   cfracTsolns = cfracPsolns = cfracT2solns = cfracFsolns = cfracNabort = 0;
+
+   eas = getEas(n, k, pm, m, aborts);		/* early abort strategy */
+
+   if (verbose > 1) {
+      fprintf(stdout, "factorBase[%u]: ", m);
+      for (j = 0; j < m; j++) {
+	 fprintf(stdout, "%u ", pm[j]);
+      }
+      putc('\n', stdout);
+      fflush(stdout);
+   }
+
+   pset(&t,      pmul(utop(k), n)); 		/* E1: Initialize */
+   pset(&r,      psqrt(t));			/* constant: sqrt(k*n) */
+   pset(&twog,   padd(r, r));			/* constant: 2*sqrt(k*n) */
+   pset(&u,      twog);				/* g + Pn */
+   pset(&lastU,  twog);
+   pset(&Qn,     pone);
+   pset(&lastQn, psub(t, pmul(r, r)));
+   pset(&An,     pone);
+   pset(&lastAn, r);
+   pset(&qn,     pzero);
+
+   do {
+F2:
+      do {
+	 if (--count == 0) goto bail;
+	 pset(&t, An);
+	 pdivmod(padd(pmul(qn, An), lastAn), n, pNull, &An);  	/* (5) */
+	 pset(&lastAn, t);
+
+	 pset(&t, Qn);
+	 pset(&Qn, padd(pmul(qn, psub(lastU, u)), lastQn));	/* (7) */
+	 pset(&lastQn, t);
+
+	 pset(&lastU, u);
+
+	 pset(&qn, pone);		/* eliminate 40% of next divmod */
+	 pset(&rn, psub(u, Qn));
+	 if (pge(rn, Qn)) {
+	    pdivmod(u, Qn, &qn, &rn);				/* (4) */
+	 }
+	 pset(&u, psub(twog, rn));				/* (6) */
+	 s = 1-s;
+
+	 e[0] = s;
+	 i = pfactorQ(&t, Qn, pm, &e[1], m, eas);  	/* E3: Factor Qn */
+	 if (i < -1) cfracNabort++;
+	 /* 
+	  * We should (but don't, yet) check to see if we can get a 
+	  * factor by a special property of Qn = 1
+	  */
+	 if (picmp(Qn, 1) == 0) {
+	    errorp(PDOMAIN, "pcfrac", "cycle encountered; pick bigger k");
+	    goto bail;			/* we ran into a cycle; give up */
+	 }
+      } while (i < 0); 			/* while not a solution */
+
+      pset(&x, An); 	/* End of Algorithm E; we now have solution: <x,t,e> */
+
+      if (i == 0) {	/* if partial */
+	 if ((sp = findSoln(oddt, t)) == (SolnPtr) 0) {
+	    cfracTsolns++;
+	    if (verbose >= 2) putc('.', stderr);
+	    if (verbose > 3) printSoln(stdout, "Partial: ","\n", pm,m,x,t,e);
+	    oddt = newSoln(n, pm, m, oddt, x, t, e);
+	    goto F2;			/* wait for same t to occur again */
+	 }
+	 if (verbose > 3) printSoln(stdout, "Partial: ", " -->\n", pm,m,x,t,e);
+	 pset(&t, pone);		/* take square root */
+         combineSoln(&x, &t, e, pm, m, n, sp);
+	 cfracT2solns++;
+	 if (verbose) putc('#', stderr);
+	 if (verbose > 2) printSoln(stdout, "PartSum: ", "", pm, m, x, t, e);
+      } else {
+	 combineSoln(&x, &t, e, pm, m, n, (SolnPtr) 0);		/* normalize */
+	 cfracPsolns++;
+	 if (verbose) putc('*', stderr);
+	 if (verbose > 2) printSoln(stdout, "Full:    ", "", pm, m, x, t, e);
+      }
+
+      /* 
+       * Crude gaussian elimination.  We should be more effecient about the
+       * binary vectors here, but this works as it is.
+       *
+       * At this point, t must be pone, or t occurred twice
+       *
+       * Loop Invariants:  e[0:h] even
+       *                   t^2 is a product of squares of primes
+       *                   b[h]->e[0:h-1] even and b[h]->e[h] odd
+       */
+      h = m+1;
+      do {
+	 --h;
+	 if (e[h]) {				/* F3: Search for odd */
+	    bp=b[h];
+	    if (bp == (SolnPtr) 0) { 		/* F4: Linear dependence? */
+	       if (verbose > 3)  {
+		  printSoln(stdout, " -->\nFullSum: ", "", pm, m, x, t, e); 
+	       }
+	       if (verbose > 2) putc('\n', stdout);
+	       b[h] = newSoln(n, pm, m, bp, x, t, e);
+	       goto F2;
+	    }
+	    combineSoln(&x, &t, e, pm, m, n, bp);
+	 }
+      } while (h != 0);
+      /*
+       * F5: Try to Factor: We have a perfect square (has about 50% chance)
+       */
+      cfracFsolns++;
+      pset(&y, t);				/* t is already sqrt'd */
+
+      switch (verbose) {
+      case 0: break;
+      case 1: putc('/', stderr); break;
+      case 2: putc('\n', stderr); break;
+      default: ;
+	 putc('\n', stderr);
+	 printSoln(stdout, " -->\nSquare:  ", "\n", pm, m, x, t, e); 
+	 fputs("x,y:     ", stdout); 
+	 fputp(stdout, x); fputs("  ", stdout);
+	 fputp(stdout, y); putc('\n', stdout);
+	 fflush(stdout);
+      }
+   } while (peq(x, y) || peq(padd(x, y), n));	/* while x = +/- y */
+
+   pset(&res, pgcd(padd(x, y), n));		/* factor found at last */
+
+   /*
+    * Check for degenerate solution.  This shouldn't happen. Detects bugs.
+    */
+   if (peq(res, pone) || peq(res, n)) {
+      fputs("Error!  Degenerate solution:\n", stdout);
+      fputs("x,y:   ", stdout); 
+      fputp(stdout, x); fputs(" ", stdout);
+      fputp(stdout, y); putc('\n', stdout);
+      fflush(stdout);
+      abort();
+   }
+
+bail:
+   if (maxCount != (unsigned *) 0) *maxCount = count;
+
+   if (b != (SolnPtr *) 0) for (j = 0; j <= m; j++) freeSoln(b[j]);
+   freeEas(eas);
+   freeSolns(oddt);
+#ifndef IGNOREFREE
+   free(e);
+   free(pm);
+#endif   
+
+   pdestroy(r);  pdestroy(twog);   pdestroy(u);  pdestroy(lastU);
+   pdestroy(Qn); pdestroy(lastQn); pdestroy(An); pdestroy(lastAn);
+   pdestroy(x);  pdestroy(y);      pdestroy(qn); pdestroy(rn);
+   pdestroy(t);  pdestroy(n);
+
+   return presult(res);
+}
+
+/*
+ * Initialization for pcfrac factoring method
+ *
+ * k      - An integer multiplier to use for n (k must be < n)
+ *            you can use findk to get a good value.  k should be squarefree
+ * m      - The number of primes to use in the factor base
+ * aborts - the number of early aborts to use
+ */
+int pcfracInit(m, k, aborts)
+   unsigned m;
+   unsigned k;
+   unsigned aborts;
+{
+   pcfrac_m      = m;
+   pcfrac_k      = k;
+   pcfrac_aborts = aborts;
+   return 1;
+}