\\ All references are to
\\ Alexander Abatzoglou, Alice Silverberg, Andrew V. Sutherland, and Angela Wong, Deterministic elliptic curve primality proving for a special sequence of numbers, Tenth Algorithmic Number Theory Symposium (ANTS X, 2012), pp. 1-20.

\\ Computes J_n by modular matrix exponentiation.
A267437(n)=([0,1,0,0;0,0,1,0;0,0,0,1;-4,8,-7,4]^n*[9;11;11;23])[1,1];

\\ Given k mod 72, find the x-coordinate of the point P_a from Table 1.
x0(n)=my(n214);if(n%3!=1,1,n24=n%24;n24==10,21,n24==4||n24==7||n24==13||n24==22,15,n==25||n==43,-633,81);

\\ Given k mod 72 and k mod 24, find the twisting parameter a from Table 1.
twist(n,n24=n%24)=-if(n24%3!=1,1,n24==10,6,n24==4||n24==7||n24==13||n24==22,5,n==25||n==43,111,17)

\\ Returns 1 if J_n is prime and 0 otherwise. Used to compute sequence A267439 in the OEIS.
is(n)=
{
	my(n72=n%72,n24=n72%24);
	if(n24%8==0 || n24==6, return(0));
	my(J=A267437(n),d=Mod(7,J)^((J+1)/4),a);
	if(d^2!=-7, return(0));
	a=twist(n72,n24);
	my(C=lift((1-3*d)/32),x=lift((3*d+7)/(56*a)*(x0(n72)-(d-7)*a/2)),z=1);
	d=0;
	for(i=1,n,
		my(s2=(x+z)^2%J,d2=(x-z)^2%J,t=s2-d2);
		if(t<0,t+=J);
		[x, z] = [s2*d2, t*(d2+C*t)]%J
	);
	if (z==0, return(0));	\\ Remark 4.2
	my(d2=(x-z)^2,t=((x+z)^2-d2)%J);
	t*(d2+C*t)%J==0;
}

\\ Example of use
select(is, [2..1000])