Comparison to ECPP?

[carpetpool]

For PARI/GP:

rset(n)={
prime_fact_r=[];
cr=1;
cr_max=round(log(n))^2;
q=2*round(log(log(n)));
	while(cr<cr_max,
	q=nextprime(q+1);
	if(n^2%q!=1 & n%q!=0,
		o=znorder(Mod(n,q));
		if((cr%(q-1))!=0,
			cr=lcm(o,cr);
			prime_fact_r=concat(prime_fact_r,q);
		); 
	); 
);
return(prime_fact_r)
}

verify_rset(n)={
new_rset=[];
old_rset=rset(n);
j=length(old_rset);
cr=1;
cr_max=round(log(n))^2;
	while(cr<cr_max,
	q=old_rset[j];
	o=znorder(Mod(n,q));
	if((cr%(q-1))!=0,
		cr=lcm(o,cr);
		new_rset=concat(new_rset,q);
	);
	j = j-1
);
return(vecsort(new_rset))
}

my_isprime(n)={
qlist=verify_rset(n);
l=length(qlist);
if(n==4,
	return(0)
);
for(i=1,l,
	q=qlist[i];
	U=floor(sqrt(eulerphi(q))*log(n));
	if( Mod(Mod(1,n)*x+U,x^q-1)^n!=x^(n%q)+U,
		return(0)
	);
);
return(1)
}

Is the test faster than ECPP and or practical for large (multi-thousand digit) numbers?

Correctness for algorithm relies on the truth of some of the findings here.

Also note I replaced steps [2] and [5] in the algorithm demonstrated in the paper as follows:

Replacement for Step [2]: r is not necessarily prime, and ord(n,r) > log(n)^2
Replacement for Step [5]: floor(sqrt(eulerphi(r))*log(n)) is floor(sqrt(eulerphi(q))*log(n)) for each prime q | r.

[CRGreathouse]

It seems to be rather fast. I generated a random thousand-digit number and it took only 47 seconds, compared to other ECPP implementations I had laying around: 355 seconds with Primo, 234 seconds with PARI/GP using primecert(,0), and 225 seconds with PARI/GP using isprime(,3).

[carpetpool]

I've modified the code with some print statements:

rset(n)={
prime_fact_r=[];
cr=1;
cr_max=round(log(n))^2;
q=2*round(log(log(n)));
	while(cr<cr_max,
	q=nextprime(q+1);
	if(n^2%q!=1 & n%q!=0,
		o=znorder(Mod(n,q));
		if((cr%(q-1))!=0,
			cr=lcm(o,cr);
			prime_fact_r=concat(prime_fact_r,q);
		); 
	); 
);
return(prime_fact_r)
}

verify_rset(n)={
new_rset=[];
old_rset=rset(n);
j=length(old_rset);
cr=1;
cr_max=round(log(n))^2;
	while(cr<cr_max,
	q=old_rset[j];
	o=znorder(Mod(n,q));
	if((cr%(q-1))!=0,
		cr=lcm(o,cr);
		new_rset=concat(new_rset,q);
	);
	j = j-1
);
return(vecsort(new_rset))
}

my_isprime(n)={
qlist=verify_rset(n);
printf("Testing using q=\n");
print(qlist);
l=length(qlist);
if(n==4,
	return(0)
);
for(i=1,l,
	q=qlist[i];
	printf("Computing (x+U)^n mod(x^%d-1,n)...\n",q);
	U=floor(sqrt(eulerphi(q))*log(n));
	if( Mod(Mod(1,n)*x+U,x^q-1)^n!=x^(n%q)+U,
		printf("%d is composite. \n",n);
		return(0)
	);
);
printf("%d is prime! \n",n);
return(1)
}

I've run it on a random 5000-digit prime number and it took under an hour:

> my_isprime(n)
Testing using q=
[47, 53, 59, 67, 73, 83]
Computing (x+U)^n mod(x^47-1,n)...
Computing (x+U)^n mod(x^53-1,n)...
Computing (x+U)^n mod(x^59-1,n)...
Computing (x+U)^n mod(x^67-1,n)...
Computing (x+U)^n mod(x^73-1,n)...
Computing (x+U)^n mod(x^83-1,n)...
... is prime!
= 1
time = 3162.419 seconds

The computational aspect of it is quite simple: If the set of primes q that are used is relatively small (like what the code produced above), the test will run faster, without as much memory. The set of primes used depends on n modulo many small primes. Fortunately for most primes, the set consists of very small primes (so odds are, most inputs will work!). However, the worst case scenario appears to be when polynomials such as n-1, n+1, n^2-1, n^4-1, n^6-1, n^12-1 have "many" small prime factors. I've generated a 5000-digit number n (here) which would be one of the "worst case" scenarios for this test: n-1 is divisible by every prime from 2 to 3500, (3499#), yet the factorization is not substantial for a 33.33% factor n-1 proof. As expected:

> my_isprime(n)
Testing using q=
[3511, 3517, 3527]
Computing (x+U)^n mod(x^3511-1,n)...
  ***   at top-level: my_isprime(p)
  ***                 ^-------------
  ***   in function my_isprime: ...f(Mod(Mod(1,n)*x+U,x^q-1)^n!=x^(n%q)+U,printf(
  ***                                                       ^---------------------
  *** _^_: the PARI stack overflows !
  current stack size: 128000000 (122.070 Mbytes)
  [hint] set 'parisizemax' to a non-zero value in your GPRC

As I mentioned before, this algorithm is heuristic (relies on unproven assumptions), but it should be an indicator that this test would be practical for (most) large numbers once it or a similar test is theoretically proven to be correct. To fix the memory problem, one should use polynomials such that the degree d is dependent on only the size of the input, NOT congruences modulo small prime numbers. A simple solution would be a variant of one of the tests here (which I have yet to work out). I'd be willing to take a shot at proving some of these ideas if I had more knowledge of number theory and group theory.

If anyone's is interested I will hopefully post another (heuristic of course) test, where the cyclotomic polynomials used in this test are replaced with polynomials defining cyclotomic subfields (and thus are cyclic polynomials).

I've also been debating about writing a C++ version of this algorithm, but I'm not sure if it's worth the effort if what the achievement is speed and efficiency.

(All suggestions and modifications to the current GP code are welcome!)