CE(k1,k2,b,n1,n2)=
{
forstep(k=k1,k2,[1,2],
for(n=n1,n2,
N=k*b^n-1;
my(s=Mod(2*polchebyshev(k*b/2,1,polchebyshev(b/2,1,2)),N));
for(i=1,n-2,s=2*polchebyshev(b,1,s/2));
if((s==0) && !ispseudoprime(N),print("k=",k," b=",b," n=",n))))
}