MultiSource/Benchmarks/MallocBench/cfrac/pcfrac.c - third_party/llvm-test-suite - Git at Google

 /*
  * 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;
 }
	/*
	* 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(kn) /
	pset(&twog, padd(r, r)); /* constant: 2sqrt(kn) */
	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;
	}