[dyrnrrin] Quick partial factorization and isprimroot in Pari/GP

\\ return the last composite in a possibly incomplete factorization, or 2 if it is a complete prime factorization.  We hope/assume that the last composite is the largest composite.
lastcomposite(ff)= my(g) ; /* sigil */ my(l=2) ; my(i) ; /* permit the possibility that the unfactored composite is not the last factor */ for(i=1 , matsize(ff)[1] , g=ff[i , 1] ; if(!ispseudoprime(g) , l=g)) ; l

\\ fast approximate factorization: try to factor, but do not try too hard, and return a possibly incomplete factorization.
\\ comments give example numbers which take the given branch for thresholds my(mpqs=2^220) ; my(ecm=2^256)
factor_partial(x)= my(mpqs=2^220) ; my(ecm=2^256) ; if(x<mpqs , factor(x) , my(ff=factorint(x , 11)) ; my(l=lastcomposite(ff)) ; if(ispseudoprime(l) , /* 2^280+10 */ ff , if(l<mpqs , /* 2^280+1 */ factor(x) , if(l>ecm , /* 2^280+5 , incomplete factorization. It would be nice if there were a way to try a little bit of ECM here. */ ff , ff=factorint(x , 9) ; l=lastcomposite(ff) ; if(ispseudoprime(l) , /* 2^280+9 */ ff , if(l<mpqs , /* 2^280+3 */ factor(x) , /* 2^280+12 incomplete factorization */ ff))))))

\\ return 0 if p-1 cannot be factored quickly, otherwise returns whether q is a primitive root of p.  p must be prime, not anything more complicated.
isprimroot_approx(q , p)= if(!ispseudoprime(p) , return(0)) ; my(ff=factor_partial(p-1)) ; my(l=lastcomposite(ff)) ; if(!ispseudoprime(l) , 0 , for(i=1 , matsize(ff)[1] , if(Mod(1 , p)==Mod(q , p)^((p-1)/ff[i , 1]) , return(0))) ; return(1))

Previously, a basic version of isprimroot.

The threshold of mpqs=2^220 in the code means don't try the quadratic sieve on numbers larger than 2^220.  Similarly the threshold ecm=2^256 means don't try the elliptic curve method on numbers larger than 2^256.  These thresholds can easily be changed.

It's unfortunate that factor_partial has to restart from the beginning when partial results look promising.

The first few numbers that factor_partial gives up on are 2^220 + {13, 133, 447, 471, 583, 625, 871, 975, 981, 1003}.  These numbers require MPQS.

? for(i=0 , +oo , p=2^220+i ; if(!ispseudoprime(lastcomposite(factor_partial(p))) , print(i)))