Frobenius and Lucas Test Pari/gp Code

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

Ldouble =
  (X)->u=vector(2);u[1]=X[1]*X[2];u[2]=(X[2]*X[2]+D*X[1]*X[1])/Mod(2,N);return(u)

Linc =
  (X)->u=vector(2);u[1]=(P*X[1]+X[2])/Mod(2,N);u[2]=(D*X[1]+P*X[2])/Mod(2,N);return(u)

Linit =
  (X,Y)->P=X;Q=Y;D=P^2-4*Q;result1=vector(2);result1[1]=Mod(1,N);result1[2]=Mod(P,N);bin=digits(N+1,2)

Lucastest =
  ()->l=length(bin);index=1;result=result1;for(Y=2,l,result=Ldouble(result);index=2*index;if(bin[Y],index=index+1;result=Linc(result)));return(result)

Phi =
  (P,X)->(X^P-1)/(X-1)

frank =
  (X)->u=1;while((fibonacci(u)%X)>0,u=u+1);return(u)

g =
  (X)->d=divisors(X);if(numdiv(X)<4,return(X));u=d[3];return(u)

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)

is =
  (n)->((Mod([1,1;1,0],n))^(n+1))[1,2]==0&&!isprime(n)&&n>1

isLucasCarmichaael =
  (X)->if(!issquarefree(X),return(0));a=factor(X);b=matsize(a);l=b[1];t=1;for(Y=1,l,t=t&&(0==((X-kronecker(5,X))%(b[Y,1]-kronecker(5,b[Y,1])))));return(t)

isLucasCarmichael =
  (X)->if(!issquarefree(X),return(0));a=factor(X);b=matsize(a);l=b[1];t=1;for(Y=1,l,t=t&&(0==((X-kronecker(5,X))%(a[Y,1]-kronecker(5,a[Y,1])))));return(t)

isPSW =
  (X)->t=(((X%5)==2)||((X%5)==3));t=t&&(lift(Mod(2,X)^(X-1))==1);t=t&&((fibonacci(X+1)%X)==0);return(t)

lcf =
  (X)->X/lpf(X)

leftover =
  (X)->u=fibonacci(X);for(Y=1,X-1,u=u/gcd(u,fibonacci(Y)));return(u)

lpf =
  (X)->u=factor(X);return(u[1,1])

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

pell =
  (X)->if(X==0,return(0));if(X==1,return(1));u0=0;u1=1;k=1;while(k<X,k=k+1;unew=P*u1+u0;u0=u1;u1=unew);return(u1)

rank =
  (X)->u=1;while((fibonacci(u)%X)>0,u=u+1);return(u)

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

test =
  (X)->u=lpf(X);u2=(u-1)/gcd(u-1,X-1);return(u2)

test2 =
  (X)->u=factor(X);u2=u[2,1];u3=lpf(X);return((u2-1)%(u3-1))

?