Pari/gp code for Khashin’s Frobenius primality test

I’ve collected in one spot the functions for using Sergey Khashin’s Frobenius primality test. The idea of a Frobenius test is due to Jon Grantham.

?  for(Y=1,100000,N=1+2*Y; galoisinit2(N,1,1);if(!issquare(n), c=findc();C=Mod(c,n);aaa=Frobeniustest(); t= ((lift(aaa[1])-A)==(0))&&((lift(aaa[2])+B)==n); if(!t&&(isprime(n)), print(n))    ))
?
? galoisinit2
%1092 = (N,a,b)->n=N;identity=vector(2);identity[1]=Mod(1,n);identity[2]=Mod(0,n);base=vector(2);base[1]=Mod(a,n);base[2]=Mod(b,n);bin=digits(n,2);A=a;B=b
? findc
%1093 = ()->s=3;while(kronecker(s,n)>-1,s=s+2);return(s)
? Frobeniustest
%1094 = ()->l=length(bin);product=identity;for(X=1,l,product=square(product,C);if(bin[X],product=mulbase(product,C)));return(product)
? square
%1095 = (X,C)->galoisproduct(X,X,C)
? galoisproduct
%1096 = (X,Y,C)->v=vector(2);v[1]=X[1]*Y[1]+C*X[2]*Y[2];v[2]=X[1]*Y[2]+X[2]*Y[1];return(v)
? mulbase
%1097 = (X,C)->galoisproduct(X,base,C)
?

I've updated my code so that it chooses a quadratic non-residue at random in [1, N/2], where N is the number to be tested:

Frobeniustest =
  ()->l=length(bin);product=identity;for(X=1,l,product=square(product,C);if(bin[X],product=mulbase(product,C)));return(product)

findc =
  ()->s=1+random((N\2)-1);while(kronecker(s,n)>-1,s=1+random((N\2)-1));return(s)

galoisinit2 =
  (N,a,b)->n=N;identity=vector(2);identity[1]=Mod(1,n);identity[2]=Mod(0,n);base=vector(2);base[1]=Mod(a,n);base[2]=Mod(b,n);bin=digits(n,2);A=a;B=b

galoisproduct =
  (X,Y,C)->v=vector(2);v[1]=X[1]*Y[1]+C*X[2]*Y[2];v[2]=X[1]*Y[2]+X[2]*Y[1];return(v)

mulbase =
  (X,C)->galoisproduct(X,base,C)

square =
  (X,C)->galoisproduct(X,X,C)