Factoring Carmichael Numbers

Dear Pari Fans,

I am pleased to announce a new method to factorise Carmichael numbers which is based on Lucas sequences. Whether it is utile or futile I shall leave up to you to judge.

Carmichael numbers, n, have the property that a^n==a mod n for all bases a. For odd prime n over a quadratic polynomial x^2-P*x+Q has the property that Q^((n-1)/2)==kronecker(Q,n) mod n. Thus the polynomial x^2-x+a^2 implies a^n==a mod n. Rather we use the polynomial x^2-(a^2-2)*x+1 over which a quick binary chain can be computed. The idea to factor Carmichael numbers is to take the g.c.d with n of x^(2*(n-1))-1 or of x^(2*n)-x^2. The algorithm takes the factors of non trivial g.c.ds and recurses the work on the two factors and so on over the respective g.c.ds. The function returns a vector of the complete factorisation of n.

{global(nm1);} \\obsolete
                                        
{Cf(n)=nm1=n-1;vecsort(CF(n));} // main function

{CF(N)=my(k,a,lu,lv,i,g,V1,V2);while(1,

k=Mod(random(1000000)+3,N); \\ pool for random numbers

a=k^2-2;lu=2;lv=a; \\ initial Lucas sequence values

forstep(i=exponent(nm1),1,-1,
if(bittest(nm1,i),
lu=lu*lv-a;lv=lv*lv-2,
lv=lu*lv-a;lu=lu*lu-2));

g=lift(gcd(lu*lv-2*a,N)); \\ test with penultimate bit = 0

if(g<2,g=lift(gcd(lv*lv-2-a,N))); \\ test with penultimate bit = 1

if(g>1, \\ if g.c.ds on non trival factors recurse

if(ispseudoprime(g),V1=[g],V1=CF(g));

g=N/g;if(ispseudoprime(g),V2=[g],V2=CF(g));

return(concat(V1,V2))));} \\ send back the factorisations \\ end of program


All timings were done on a single core of a 1.7GHz Celeron laptop with 1GB of memory allocated. For the 3808 digit number [1] with 662 factors:
23.0 seconds using Pari/GP's factor().
21.6 seconds using the above code.

For the 43779 digit number [1] with 5616 factors 5,443.5 seconds. The function factor() ran out of memory after 21 hours and did not complete the task it was given.

With a smaller amount of allocated memory the 16432 digit number [2] with 8 factors takes about 230 seconds.

This new algorithm is easily made parallel. It would be wise to switch to "factor()" below a certain length threshold on recursion.

[1] https://github.com/drazioti/Carmichael/tree/master/Tables/800_to_10000_factors

[2] https://listserv.nodak.edu/cgi-bin/wa.exe?A2=NMBRTHRY;a892cf10.2106&S=

Best,

Paul Underwood