######################################################################
## Perrin.txt: Save this file as  Perrin.txt  to use it,             #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read Perrin.txt <Enter>                           #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Robert Dougherty-Bliss and                             #
#and Doron Zeilberger, Rutgers University ,                          #
##  DoronZeil at gmail dot com                                       # 
######################################################################

#with(LinearAlgebra):
print(`First Written: Feb. 2023: tested for Maple 2020 `):

print():
print(`This is Perrin.txt, A Maple package`):
print(`accompanying Robert Dougherty-Bliss and Doron Zeilberger's  article: `):

print(` Lots and Lots of Perrin-Type Primality Tests and Their Pseudo-Primes `):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/Perrin.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):

print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):

print(`------------------------------------`):

print(`For a list of the supporting functions type: ezra1();`):
print():
print(`For help with a specific procedure, type:`):
print(`"ezra(procedure_name);"`):

 print(`------------------------------------`):
print(`For a list of the Pre-Computed functions type: ezraPC();`):
print():
print(`For help with a specific procedure, type:`):
print(`"ezra(procedure_name);"`):

 print(`------------------------------------`):



 print(`------------------------------------`):
print(`For a list of the functions trying to prove infinite families of pseudo-primes, type:  ezraProof();`):
print():
print(`For help with a specific procedure, type:`):
print(`"ezra(procedure_name);"`):

 print(`------------------------------------`):





ezraProof:=proc()
if args=NULL then

print(`The Proof procedures for proving infinite families of pseudo-primes, are`):
print(`  Binet, Binet1, Find22, Find23 `):

else
ezra(args):

fi:
end:


ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are`):
print(`  Comp, MatC, FindPP2g,  FindPP3g, FindPP4g,GuessRec, Pow1,  PPgF, RtoC, SeqFromRec, ValC, ValCp,`):

else
ezra(args):

fi:
end:



ezraOld:=proc()
if args=NULL then

print(`The OLD procedures no longer needed, are: `):
print(`  PatternDB`):

else
ezra(args):

fi:
end:

ezraPC:=proc()
if args=NULL then

print(`The procedures with Pre-computed data are`):
print(`  DBK, DBZ, Champions1, HitParade1, PDB `):

else
ezra(args):

fi:
end:


ezra:=proc()
if args=NULL then

print(` Perrin.txt: A Maple package for generating Perrin-like Primality Tests and for detecting pseudo-primes  `):

print(`The MAIN procedures are: FindPP1, FindPP2, FindPP3, GenPerrin, PPg, PtestG, GenPerrinPatternSearch, PT, PTall `):
print(` `):



elif nargs=1 and args[1]=Binet then
print(`Binet(f,x,n): the explicit formula that handles the n-th coefficient of the Taylor expansion of f if it is a primaility test. Try:`):
print(`Binet( (2-2*x)/(1-2*x-x^2),x,n);`):


elif nargs=1 and args[1]=Binet1 then
print(`Binet1(f,x,n): Same as Binet(f,x,n) but in more convnenient format for SBE. Try: `):
print(`Binet1( (2-2*x)/(1-2*x-x^2),x,n);`):


elif nargs=1 and args[1]=Champions1 then
print(`Champions1(x): The pre-computed output of PTall(5,5,1000,1000000,1,x): In other words all the primality tests, given in terms of their generating functions in the variable x, that`):
print(`have at most one pseudo-prime <=1000, followed by those less than a million, inspired by compositions with at most`):
print(`five parts that add0-up to <=5. Type:`):
print(`Champions1(x);`):

elif nargs=1 and args[1]=Comp then
print(`Comp(k,N): The set of vectors [a1,...,a[k]] with non-negative integers that add-up to N. Try:`):
print(`Comp(2,3);`):



elif nargs=1 and args[1]=DBK then
print(`DBK(x): The even newer primality test, due to Robert Dougherty-Bliss and Manuel Kauers  whose smallest pseudo-prime is    2,260,550,373 = 3 * 103 * 107 * 68371 . Try:`):
print(`DBK(x);`):

elif nargs=1 and args[1]=DBZ then
print(`DBZ(x): The new primality test whose smallest pseudo-prime is 1531398, namely (-3*x^4-5*x^2-6*x+7)/(-4*x^7-x^4-x^2-x+1). Try:`):
print(`DBZ(x);`):

elif nargs=1 and args[1]=Find22 then
print(`Find22(f,x): if a(n) has the form a1^n+a2^n then it finds an expression for a(2*n) in terms of other things. Try:`):
print(`Find22(PDB(x)[3][2],x);`):


elif nargs=1 and args[1]=Find23 then
print(`Find23(f,x): if a(n) has the form a1^n+a2^n then it finds an expression for a(3*n) in terms of other things. Try:`):
print(`Find23(PDB(x)[3][2],x);`):

elif nargs=1 and args[1]=FindPP1 then
print(`FindPP1(f,s,k,K): All the pseudo-primes for the test inspired by the rational function f of s that happens to be L[n],of the form p^k, where p is a one of the first K primes.`):
print(`Try:`):
print(`FindPP1(GenPerrin([0,1,1],s),s,2,1000);`):

elif nargs=1 and args[1]=FindPP2 then
print(`FindPP2(f,x,K): All the pseudo-primes for the test inspired by f(x)  p*q, where p,q<=ithprime(K)  are  `):
print(`FindPP2( (x^2-2*x+3)/(-4*x^3+x^2-x+1),x,1000);`):

elif nargs=1 and args[1]=FindPP2g then
print(`FindPP2g(f,s,K): Given a rational function f of s and a list of two integers K=[K1,K2] finds`):
print(`all the pseudo-primes that are of the form ithprime(i1)*ithprime(i2) with i1<=i2 and i1<=K1 and i2<=K2. Try:`):
print(`FindPP2g(GenPerrin([0,2,3,1],s),s,[10,20]);`):


elif nargs=1 and args[1]=FindPP3 then
print(`FindPP3(f,x,K): All the pseudo-primes for the test inspired by f(x) of the form p*q*r, where p,q,r<=ithprime(K) `):
print(`FindPP3( (-x^3-2*x^2+4)/(-x^4-x^3-x^2+1),x,1000);`):


elif nargs=1 and args[1]=FindPP3g then
print(`FindPP3g(f,s,K): Given a rational function f of s and a list of three integers K=[K1,K2,K3] finds`):
print(`all the pseudo-primes that are of the form ithprime(i1)*ithprime(i2)*ithprime(i3) with i1<=i2<=i3 and i<=K1 and i2<=K2, i3<=K3. Try:`):
print(`FindPP3g(GenPerrin([0,2,3,1],s),s,[10,20,30]);`):


elif nargs=1 and args[1]=FindPP4g then
print(`FindPP4g(f,s,K): Given a rational function f of s and a list of four integers K=[K1,K2,K3,K4] finds`):
print(`all the pseudo-primes that are of the form ithprime(i1)*ithprime(i2)*ithprime(i3)*ithprime(i4) with i1<=i2<=i3<=i4 and i<=K1 and i2<=K2, i3<=K3 and i4<=K4. Try:`):
print(`FindPP4g(GenPerrin([0,2,3,1],s),s,[10,10,10,10]);`):

elif nargs=1 and args[1]=GenPerrin then
print(`GenPerrin(L, x): The general Perrin-style rational function with a monic denominator. `):
print(`To get the original Perrin generating function, try:`):
print(`GenPerrin([0, 1, 1],x);`):
print(`To get the generic generating function for a quartic denominator, try:`):
print(`GenPerrin([a, b, c, d], x);`):


elif nargs=1 and args[1]=GuessRec then
print(`GuessRec(L): inputs a sequence L and tries to guess`):
print(`a recurrence operator with constant cofficients `):
print(`satisfying it. It returns the initial  values and the operator`):
print(` as a list. For example try:`):
print(`GuessRec([1,1,1,1,1,1]);`):

elif nargs=1 and args[1]=HitParade1 then
print(`HitParade1(x): The ordered list according to the number of pseudo-primes coming out of Champions1(x) (q.v.) according to the number of pseudo-primes<=100000.`):
print(`As you can see, the original Perrin test, with 2 pseudo-primes  is the champion. Try:`):
print(`HitParade1(x);`):


elif nargs=1 and args[1]=Pow1 then
print(`Pow1(M,k,p): M^k mod p. Try:`):
print(`Pow1([[1,1],[0,1]],100,1002);`):


elif nargs=1 and args[1]=PPg then
print(`PPg(f,s,L): All the pseudo-primes <=L in the primality set induced by the rational function f of s`):
print(`PPg(GenPerrin([0,1,1],s),s,10000);`):


elif nargs=1 and args[1]=PPgF then
print(`PPgF(f,s,L): Memory efficient version of Ppg(f,s,L). Try:`):
print(`PPgF(GenPerrin([0,1,1],s),s,10000);`):


elif nargs=1 and args[1]=PT then
print(`PT(k,N,K1,K2,ma,x): All the compositions C of N into k parts such that  f:=GenPerrin(C,x) has <=ma pseudo-primes <=K1. It returns the set`):
print(`of [C,f,Pseudoprimes<=K1, PseuPrimes<-K2]. Try:`):
print(`PT(3,2,2000,10000,3,x);`):

elif nargs=1 and args[1]=PTall then
print(`PTall(k,N,K1,K2,ma,x): inputs positive intgeres k, N, K1,K2,ma, x and applies PT(k1,N1,K1,K2,ma,x) for all 2<=k1<k, 1<=N1<=N. Try:`):
print(`PTall(3,3,2000,10000,1,x);`):



elif nargs=1 and args[1]=PtestG then
print(`PtestG(f,x,n,p0): Testing pos. integer n according to the primality test where the p-th term mod p should be p0,`):
print(`inspired by the rational generating function f of x`):
print(`PtestG(GenPerrin([0,1,1],s),s,1001,0);`):



elif nargs=1 and args[1]=GenPerrinPatternSearch then
print(`GenPerrinPatternSearch(max_deg, max_coeff): Search for Perrin-style primality tests with explicit families of pseudoprimes. `):
print(`Try:`):
print(`GenPerrinPatternSearch(2, 2);`):
print(`GenPerrinPatternSearch(3, 2);`):
print(`To see pre-cached results, see ezra(PatternDB).`):

elif nargs=1 and args[1]=PDB then
print(`PDB(x): A pre-computed data-base of fourteen pairs [p,f] where p is 2 or 3 and f is a rational function of x (of degree 2 or 3) that`):
print(`conjecturally have {p^i} as pseudo-primes. Try:`):
print(`PDB(x);`):

elif nargs=1 and args[1]=PatternDB then
print(`PatternDB(): Pre-computed database of Perrin-style primality tests with explicit families of pseudoprimes. `):
print(`Format is [p, L], where p^2, p^3, p^4, ... are the conjecture pseudoprimes for the test L.`):
print(`Try:`):
print(`PatternDB()[1]`):



elif nargs=1 and args[1]=MatC then
print(`MatC(C): The matrix corresponding to the recurrence C=[INI,Rec]. Try: `):
print(`MatC([[1,1],[1,1]]);`):


elif nargs=1 and args[1]=RoC then
print(`RtoC(R,t): Given a rational function R(t) finds the`):
print(`C-finite description of its sequence of coefficients`):
print(`For example, try:`):
print(`RtoC(1/(1-t-t^2),t);`):


elif nargs=1 and args[1]=SeqFromRec then
print(`SeqFromRec(S,N): Inputs S=[INI,ope]`):
print(`where INI is the list of initial conditions, a ope a list of`):
print(`size L, say, and a recurrence operator ope, codes a list of`):
print(`size L, finds the first N0 terms of the sequence satisfying`):
print(`the recurrence f(n)=ope[1]*f(n-1)+...ope[L]*f(n-L).`):
print(`For example, for the first 20 Fibonacci numbers,`):
print(`try: SeqFromRec([[1,1],[1,1]],20);`):


elif nargs=1 and args[1]=ValCp then
print(`ValCp(C,n,p): The n-th term of the C-finite sequence mod p, given by the recurrence C=[INI,REC]. Should be the same (but faster) than`):
print(`ValC(C,n) mod p.`):
print(`ValCp([[1,1],[1,1]],30,101);`):

else
 print(`There is no such thing as`, args):

fi:


end:
 



#####################################################


with(linalg):
with(numtheory):


GenPerrinPatternSearch := proc(max_deg, max_coeff)
    local res, d, sets, k, cs, f, start, finish, time_sum, count, next_time, pseudos, p:
    kernelopts(printbytes=false):
    res := []:
    for d from 2 to max_deg do
        sets := seq([seq(-max_coeff..max_coeff)], k=1..d):

        k := 1:
        time_sum := 0:
        # First one can't be negative, last one can't be zero.
        count := (max_coeff + 1) * (2 * max_coeff) * (2 * max_coeff + 1)^(d - 2):

        # Smallest multiple of 5 greater than the elapsed time.
        next_time := 5:

        printf("Starting d = %d\n", d):
        printf("%d iterations\n", count):

        for cs in Iterator[CartesianProduct](sets) do
            start := time():

            cs := convert(cs, list):

            if cs[..d-1] = [0 $ d-1] then
                next:
            fi:

            if cs[1] < 0 or cs[-1] = 0 then
                next:
            fi:

            f := GenPerrin(cs, x):
            pseudos := pseudoPrimes(f, x, 5^8 + 1):
            for p in [2, 3, 5] do
                if andmap(k -> member(p^k, pseudos), [seq(2..8)]) then
                    res := [op(res), [p, cs]]:
                fi:
            od:

            finish := time():
            time_sum := time_sum + finish - start:

            if time_sum > next_time then
                next_time := 5 * (floor(round(time_sum) / 5) + 1):
                printf("Iteration %d/%d\n", k, count):
                printf("Est. %ds remaining\n", round((count - k) * (time_sum / k))):
            fi:

            k := k + 1:
        od:
    od:

    res:
end:



#PDB(x): A pre-computed data-base of fourteen pairs [p,f] where p is 2 or 3 and f is a rational function of x (of degree 2 or 3) that
#conjecturally have {p^i} as pseudo-primes. Try:
#PDB(x);
PDB:=proc(x):
[
[2, (2-2*x)/(-x^2-2*x+1)], 
[2, (2-x)/(-2*x^2-x+1)],  
[3, (2-2*x)/(-2*x^2-2*x+1)], 
[2, (2*x^2+3)/(2*x^3+2*x^2+1)], 
[2, (-2*x^2+3)/(2*x^3-2*x^2+1)], 
[2, (-2*x^2+3)/(x^3-2*x^2+1)], 
[2, (2*x^2+3)/(x^3+2*x^2+1)], 
[3, (-2*x^2-2*x+3)/(-x^3-2*x^2-x+1)], 
[2, (2*x^2+3)/(-x^3+2*x^2+1)], 
[2, (2*x^2+3)/(-2*x^3+2*x^2+1)], 
[2, (-x^2-4*x+3)/(-2*x^3-x^2-2*x+1)], 
[3, (-x^2-2*x+3)/(-2*x^3-x^2-x+1)], 
[2, (-2*x^2+3)/(-2*x^3-2*x^2+1)], 
[2, (-2*x^2-2*x+3)/(-2*x^3-2*x^2-x+1)]
]:
end:

PatternDB := proc()
[[2, [2, 1]], [2, [1, 2]], [3, [2, 2]], [2, [0, -2, -2]], [2, [0, 2, -2]],
[2, [0, 2, -1]], [2, [0, -2, -1]], [3, [1, 2, 1]], [2, [0, 2, 1]],
[3, [1, -1, 1]], [5, [1, -1, 1]], [5, [2, -2, 1]], [2, [0, -2, 1]],
[2, [0, -2, 2]], [2, [1, -2, 2]], [3, [1, -2, 2]], [5, [1, -2, 2]],
[2, [2, -1, 2]], [2, [1, 0, 2]], [2, [2, 1, 2]], [3, [1, 1, 2]],
[2, [0, 2, 2]], [2, [1, 2, 2]]]:
end:


#FindPP1(f,s,k,K): All the pseudo-primes for the test inspired by the rational function f of s that happens to be L[n],of the form p^k, where p is a one of the first K primes.
#Try:
FindPP1:=proc(f,s,k,K) local f1,i,p,gu,mu,mu1,p0:

f1:=taylor(f,s=0,ithprime(K)^(k-1)+10):


mu:=[seq(coeff(f1,s,i),i=1..ithprime(K)^(k-1)+9)]:



mu1:={seq(mu[ithprime(i)] mod ithprime(i),i=10..20)}:


if nops(mu1)<>1 then
 print(`Not a primality test`):
 RETURN(FAIL):
fi:

p0:=mu1[1]:


gu:=[]:
for i from 1 to K do
p:=ithprime(i):
if mu[p^(k-1)] mod p^k=p0 then
 gu:=[op(gu),p^k]:
fi:
od:
gu:
end:




#GuessRec1(L,d): inputs a sequence L and tries to guess
#a recurrence operator with constant cofficients of order d
#satisfying it. It returns the initial d values and the operator
# as a list. For example try:
#GuessRec1([1,1,1,1,1,1],1);
GuessRec1:=proc(L,d) local eq,var,a,i,n:

if nops(L)<=2*d+2 then
 print(`The list must be of size >=`, 2*d+3 ):
 RETURN(FAIL):
fi:

var:={seq(a[i],i=1..d)}:

eq:={seq(L[n]-add(a[i]*L[n-i],i=1..d),n=d+1..nops(L))}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
else
 RETURN([[op(1..d,L)],[seq(subs(var,a[i]),i=1..d)]]):
fi:

end:




#GuessRec(L): inputs a sequence L and tries to guess
#a recurrence operator with constant cofficients 
#satisfying it. It returns the initial  values and the operator
# as a list. For example try:
#GuessRec([1,1,1,1,1,1]);
GuessRec:=proc(L) local gu,d:

for d from 1 to trunc(nops(L)/2)-2 do
 gu:=GuessRec1(L,d):
 if gu<>FAIL then
   RETURN(gu):
 fi:
od:
FAIL:

end:

#SeqFromRec(S,N): Inputs S=[INI,ope]
#where INI is the list of initial conditions, a ope a list of
#size L, say, and a recurrence operator ope, codes a list of
#size L, finds the first N0 terms of the sequence satisfying
#the recurrence f(n)=ope[1]*f(n-1)+...ope[L]*f(n-L).
#For example, for the first 20 Fibonacci numbers,
#try: SeqFromRec([[1,1],[1,1]],20);
SeqFromRec:=proc(S,N) local gu,L,n,i,INI,ope:
INI:=S[1]:ope:=S[2]:
if not type(INI,list) or not type(ope,list) then
 print(`The first two arguments must be lists `):
 RETURN(FAIL):
fi:
L:=nops(INI):
if nops(ope)<>L then
 print(`The first two arguments must be lists of the same size`):
RETURN(FAIL):
fi:

if not type(N,integer) then
 print(`The third argument must be an integer`, L):
 RETURN(FAIL):
fi:

if N<L then
 RETURN([op(1..N,INI)]):
fi:

gu:=INI:

for n from 1 to N-L do
 gu:=[op(gu),expand(add(gu[nops(gu)-i+1]*ope[i],i=1..L))]:
od:

gu:
end:



#RtoC(R,t): Given a rational function R(t) finds the
#C-finite description of its sequence of coefficients
#For example, try:
#RtoC(1/(1-t-t^2),t);
RtoC:=proc(R,t) local S1,S2,D1,R1,d,f1,i,a:
R1:=normal(R):
D1:=denom(R1):
d:=degree(D1,t):
R1:=rem(numer(R1),D1,t)/D1:
a:=coeff(D1,t,0):
S2:=[seq(-coeff(D1,t,i)/a,i=1..degree(D1,t))]:
f1:=taylor(R1,t=0,d+2):
S1:=[seq(coeff(f1,t,i),i=1..d)]:
[S1,S2]:
end:

# Evaluate the nth term of the generalized lucas sequence with respect to p(x)
# in terms of the coefficients of p(x).
# See: https://en.wikipedia.org/wiki/Newton%27s_identities
newton:=proc(p, x, n)
    option remember:
    local d, k:
    d := degree(p, x):

    if n = 0 then
        return d:
    fi:

    -n * coeff(p, x, d - n) - expand(add(coeff(p, x, d - n + k) * newton(p, x, k), k=1..n-1)):
end:

CtoR:=proc(S,t) local D1,i,N1,L1,f,f1,L:
if not (type(S,list) and  nops(S)=2 and type(S[1],list) and type(S[2],list)
        and nops(S[1])=nops(S[2]) and type(t, symbol) ) then
   print(`Bad input`):
   RETURN(FAIL):
fi:

D1:=1-add(S[2][i]*t^i,i=1..nops(S[2])):
N1:=add(S[1][i]*t^(i-1),i=1..nops(S[1])):
L1:=expand(D1*N1):


L1:=add(coeff(L1,t,i)*t^i,i=0..nops(S[1])-1):
f:=L1/D1:
L:=degree(D1,t)+10:
f1:=taylor(f,t=0,L+1):
if expand([seq(coeff(f1,t,i),i=0..L)])<>expand(SeqFromRec(S,L+1)) then
print([seq(coeff(f1,t,i),i=0..L)],SeqFromRec(S,L+1)):
 RETURN(FAIL):
else
 RETURN(f):
fi:
end:

GenPerrin:=proc(L, x)
    local p, rec, ics, k:
    p := x^(nops(L)) - add(L[k] * x^(nops(L) - k), k=1..nops(L)):
    rec := L:
    ics := [seq(newton(p, x, k), k=0..nops(L)-1)]:
    CtoR([ics, rec], x):
end:




#PPgOld(f,s,L): All the pseudo-primes <=L in the primality set induced by the rational function f of s
#PPgOld((1+s^2)/(1-s-s^2),x,10000);
PPgOld:=proc(f,s,L) local gu,i,mu,lu,D1,f1,k0:
if L<=1000 then
 print(` Make `, L, `bigger `):
 RETURN(FAIL):
fi:

for k0 from 1 while ithprime(k0)<L do od:
k0:=k0-1:

D1:=denom(f):
f1:=rem(numer(f),D1,s):
f1:=f1/D1:
f1:=taylor(f1,s=0,L+5):
gu:=[seq(coeff(f1,s,i),i=1..L)]:

lu:={seq(gu[ithprime(i)] mod ithprime(i),i=max(1,k0-4)..k0)} :

if nops(lu)<>1 then
 RETURN(FAIL):
fi:

lu:=lu[1]:

mu:=[]:

for i from 4 to L do
 if not isprime(i) and gu[i] mod i=lu then
  mu:=[op(mu),i]:
 fi:
od:
mu:


end:





#Pow1(M,k,p): M^k mod p
Pow1:=proc(M,k,p) local d,M1,i,j:
option remember:
d:=nops(M):

if k=0 then
 RETURN([seq([0$(i-1),1,0$(d-i)],i=1..d)]):
fi:

if k mod 2=1 then
 M1:=convert(multiply(Pow1(M,k-1,p),M),list):
 RETURN([seq([seq(M1[i][j] mod p,j=1..d)],i=1..d)]):
else
 M1:=Pow1(M,k/2,p):
 M1:=convert(multiply(M1,M1),list):
 RETURN([seq([seq(M1[i][j] mod p,j=1..d)],i=1..d)]):
fi:


end:


#MatC(C): The matrix corresponding to the recurrence C=[INI,Rec]. Try
#MatC([[1,1],[1,1]]);
MatC:=proc(C) local INI1,k,i,M:
 INI1:=C[1]:
 INI1:=[seq(INI1[nops(INI1)+1-i],i=1..nops(INI1))]:
 k:=nops(C[2]):
 M:=[C[2],seq([0$(i-1),1,0$(k-i)],i=1..k-1)]:
M,INI1:
end:


#ValCp(C,n,p): The n-th term of the C-finite sequence mod p, given by the recurrence C=[INI,REC]. Should be the same (but faster) than
#ValCp(C,n) mod p.
#ValCp([[1,1],[1,1]],30,101);
ValCp:=proc(C,n,p) local gu,M,k,i1,Ini:
gu:=MatC(C):
M:=gu[1]:
k:=nops(M):
Ini:=gu[2]:
Ini:=[seq([Ini[i1]],i1=1..k)]:
multiply(Pow1(M,n-k,p),Ini)[1,1] mod p:
end:

#PtestG(f,x,n,p0): Testing pos. integer n according to the primality test where the p-th term mod p should be p0,
#inspired by the rational generating function f of x
PtestG:=proc(f,x,n,p0) local C:

C:=RtoC(f,x):
evalb(ValCp(C,n,n)=p0):
end:
 

#Comp(k,N): The set of vectors [a1,...,a[k]] with non-negative integers that add-up to N. Try:
#Comp(2,3);
Comp:=proc(k,N) local gu,ak,lu,lu1:
option remember:
if k=0 then
 if N=0 then
  RETURN({[]}):
else
 RETURN({}):
fi:
fi:

if N<0 then
 RETURN({}):
fi:

gu:={}:

for ak from 0 to N do
 lu:=Comp(k-1,N-ak):
 gu:=gu union {seq([op(lu1),ak],lu1 in lu)}:
od:
gu:

end:





#PT(k,N,K1,K2,ma,x): All the compositions C of N into k parts such that  f:=GenPerrin(C,x) has <=ma pseudo-primes <=K1. It returns the
#set
#{[C,f,Pseudoprimes<=K1, Pseuprimes<=K2]}. Try:
#PT(3,2,1000,10000,3,x);
PT:=proc(k,N,K1,K2,ma,x) local lu,f,C,gu,mu:
lu:=Comp(k,N):
gu:={}:

for C in lu do
 if C[-1]<>0 then
 f:=GenPerrin(C,x):
  mu:=PPg(f,x,K1):
 if mu<>FAIL then
 if nops(mu)<=ma then
  gu:=gu union {[C,f,mu, PPg(f,x,K2)]}:
 fi:
fi:
fi:
od:
gu:
end:


#PTall(k,N,K1,K2,ma,x): inputs positive intgeres k, N, K1,K2,ma, x and applies PT(k1,N1,K1,K2,ma,x) for all 2<=k1<k, 1<=N1<=N. Try:
#PTall(3,3,1000,10000,1,x);
PTall:=proc(k,N,K1,K2,ma,x) local k1,N1,lu,gu:

gu:=[]:

for k1 from 2 to k do
for N1 from 2 to N do
 lu:=PT(k1,N1,K1,K2,ma,x):
 if lu<>{} then
  gu:=[op(gu),lu]:
 fi:
od:
od:
gu:


end:

#Champions1(x): The pre-computed output of PTall(5,5,1000,1000000,1,x): In other words all the primality tests that
#have at most one pseudo-prime <=1000, followed by those less than a million, inspired by compositions with at most
#five parts that add0-up to <=5. Try:
#Champions1(x);
Champions1:=proc(x):
[{[[1, 1], (2-x)/(-x^2-x+1), [705], [705, 2465, 2737, 3745, 4181, 5777, 6721, 
10877, 13201, 15251, 24465, 29281, 34561, 35785, 51841, 54705, 64079, 64681, 
67861, 68251, 75077, 80189, 90061, 96049, 97921, 100065, 100127, 105281, 113573
, 118441, 146611, 161027, 162133, 163081, 179697, 186961, 194833, 197209, 
209665, 219781, 228241, 229445, 231703, 252601, 254321, 257761, 268801, 272611,
283361, 302101, 303101, 327313, 330929, 399001, 430127, 433621, 438751, 447145,
455961, 489601, 490841, 497761, 512461, 520801, 530611, 556421, 597793, 618449,
635627, 636641, 638189, 639539, 655201, 667589, 687169, 697137, 722261, 741751,
851927, 852841, 853469, 920577, 925681, 930097, 993345, 999941]]}, {[[0, 1, 1],
(-x^2+3)/(-x^3-x^2+1), [], [271441, 904631]]}, {[[0, 1, 2], (-x^2+3)/(-2*x^3-x^
2+1), [42], [42, 1729, 6790, 15457, 19045, 26866, 32041, 214367]], [[1, 1, 1],
(-x^2-2*x+3)/(-x^3-x^2-x+1), [182], [182, 25201, 54289, 63618, 194390, 750890,
804055]]}, {[[0, 3, 1], (-3*x^2+3)/(-x^3-3*x^2+1), [114], [114, 1474, 5833, 
19345, 41154, 74670, 96139, 146081, 146611, 282133, 286903, 294409, 488881, 
653333, 658801, 859951, 876961, 929633]], [[2, 0, 2], (3-4*x)/(-2*x^3-2*x+1), [
245], [245, 1634, 19966, 96178]]}, {[[0, 1, 4], (-x^2+3)/(-4*x^3-x^2+1), [], [
3481, 11881, 113399]]}, {[[0, 1, 1, 1], (-x^3-2*x^2+4)/(-x^4-x^3-x^2+1), [308],
[308, 1155, 49196, 552599]]}, {[[1, 2, 0, 1], (-4*x^2-3*x+4)/(-x^4-2*x^2-x+1),
[4], [4, 5012, 49395, 135916, 599451]]}, {[[0, 1, 1, 3], (-x^3-2*x^2+4)/(-3*x^4
-x^3-x^2+1), [28], [28, 1204, 1540, 5502, 54642, 92833, 106925, 170625, 184459,
354061, 518714, 999635]], [[0, 1, 2, 2], (-2*x^3-2*x^2+4)/(-2*x^4-2*x^3-x^2+1),
[646], [646, 1729, 3961, 4066, 5329, 20595]], [[0, 3, 1, 1], (-x^3-6*x^2+4)/(-x
^4-x^3-3*x^2+1), [45], [45, 4035, 4927, 5915, 14266, 27945, 32805, 37249, 43965
, 101387, 161098, 295245, 317385, 594315]], [[1, 1, 2, 1], (-2*x^3-2*x^2-3*x+4)
/(-x^4-2*x^3-x^2-x+1), [9], [9, 3819, 343441]], [[2, 0, 0, 3], (4-6*x)/(-3*x^4-\
2*x+1), [49], [49, 15878]]}, {[[0, 1, 1, 0, 1], (-2*x^3-3*x^2+5)/(-x^5-x^3-x^2+
1), [], [2877, 12609, 142546, 227837]], [[1, 1, 0, 0, 1], (-3*x^2-4*x+5)/(-x^5-
x^2-x+1), [45], [45, 1241, 3721, 8932, 10318, 18497, 19686, 98421, 682892]], [[
2, 0, 0, 0, 1], (5-8*x)/(-x^5-2*x+1), [77], [77, 2809, 3097, 32761, 46134, 
53734, 469682, 772593]]}, {[[0, 1, 1, 0, 2], (-2*x^3-3*x^2+5)/(-2*x^5-x^3-x^2+1
), [], [1909, 45485, 675455]], [[1, 1, 0, 0, 2], (-3*x^2-4*x+5)/(-2*x^5-x^2-x+1
), [], [5083, 6517, 36481, 697461, 888501]], [[2, 0, 0, 0, 2], (5-8*x)/(-2*x^5-\
2*x+1), [841], [841, 1742, 7205, 445439]], [[2, 0, 0, 1, 1], (-x^4-8*x+5)/(-x^5
-x^4-2*x+1), [731], [731, 1369, 4687, 5041, 8927, 9409, 297733, 586861, 736291]
]}, {[[0, 1, 1, 0, 3], (-2*x^3-3*x^2+5)/(-3*x^5-x^3-x^2+1), [], [1121, 5170, 
357850, 574262]], [[0, 1, 1, 1, 2], (-x^4-2*x^3-3*x^2+5)/(-2*x^5-x^4-x^3-x^2+1)
, [], [24649, 91970, 197537, 258531]], [[0, 1, 1, 2, 1], (-2*x^4-2*x^3-3*x^2+5)
/(-x^5-2*x^4-x^3-x^2+1), [705], [705, 2465, 2737, 3745, 4181, 5777, 6721, 10877
, 13201, 15251, 24465, 29281, 34561, 35785, 51841, 54705, 64079, 64681, 67861,
68251, 75077, 80189, 90061, 96049, 97921, 100065, 100127, 105281, 113573, 
118441, 146611, 161027, 162133, 163081, 179697, 186961, 194833, 197209, 209665,
219781, 228241, 229445, 231703, 252601, 254321, 257761, 268801, 272611, 283361,
302101, 303101, 327313, 330929, 399001, 430127, 433621, 438751, 447145, 455961,
489601, 490841, 497761, 512461, 520801, 530611, 556421, 597793, 618449, 635627,
636641, 638189, 639539, 655201, 667589, 687169, 697137, 722261, 741751, 851927,
852841, 853469, 920577]], [[0, 1, 2, 0, 2], (-4*x^3-3*x^2+5)/(-2*x^5-2*x^3-x^2+
1), [169], [169, 1271, 2809, 2831, 2834, 5157, 10830, 23933, 54289, 421029]], [
[0, 1, 3, 0, 1], (-6*x^3-3*x^2+5)/(-x^5-3*x^3-x^2+1), [9], [9, 1342, 21350, 
568763]], [[1, 1, 1, 0, 2], (-2*x^3-3*x^2-4*x+5)/(-2*x^5-x^3-x^2-x+1), [49], [
49, 2695, 13958, 28273, 151585, 569427, 864973]], [[2, 0, 0, 0, 3], (5-8*x)/(-3
*x^5-2*x+1), [121], [121, 1778, 2209, 41354, 664287]]}]:
end:


#HitParade1(x): The ordered list according to the number of pseudo-primes coming out of Champions1(x) (q.v.) according to the number of pseudo-primes<=100000.
#As you can see, the original Perrin test, with 2 pseudo-primes  is the champion. Try:
#HitParade1(x);
HitParade1:=proc(x) local gu,i,mu,mu1,T:
gu:=Champions1(x):
gu:=[seq(op(gu[i]),i=1..nops(gu))]:
gu:=[seq([gu[i][2],gu[i][4]],i=1..nops(gu))]:
mu:={seq(nops(gu[i][2]),i=1..nops(gu))}:

for mu1 in mu do
T[mu1]:={}:
od:

for i from 1 to nops(gu) do
T[nops(gu[i][2])]:=T[nops(gu[i][2])] union {gu[i]}:
od:
mu:=sort(convert(mu,list)):

[seq(op(T[mu[i]]),i=1..nops(mu))]:

end:


#DBZ(x):The Dougherty-Bliss-Zeilberger primaility test
DBZ:=proc(x):(-3*x^4-5*x^2-6*x+7)/(-4*x^7-x^4-x^2-x+1):end:


#DBK(x):The Dougherty-Bliss-Kauers primaility test
DBK:=proc(x): GenPerrin([2, 0, 0, 8, 9, 3], x):end:

#start proving infinite families

#Binet(f,x,n): the explicit formula that handles the n-th coefficient of the Taylor expansion of f if it is a primaility test. Try:
#Binet( (2-2*x)/(1-2*x-x^2),x,n):
Binet:=proc(f,x,n) local f1, top,bot,gu,L,a,i,S,bot1,ka,mu,L1:
f1:=taylor(f,x=0,1001):
L:=[seq(coeff(f1,x,i),i=1..1000)]:

top:=numer(f):
bot:=denom(f):

gu:=[solve(bot,x)]:


bot1:=simplify(expand(mul(1-x/gu[i],i=1..nops(gu)))):


ka:=normal(bot1/bot):

if not type(ka,numeric) then
print(ka):

 RETURN(FAIL):
fi:

top:=ka*top:
bot:=bot1:


for i from 1 to nops(gu) do
mu[i]:=simplify(-1/gu[i]*subs(x=gu[i],top)/subs(x=gu[i],diff(bot,x))):
od:

a:=add(mu[i]/gu[i]^n,i=1..nops(gu)):

L1:=simplify([seq(subs(n=i,a),i=1..10)]):

if L1<>[op(1..10,L)] then
 RETURN(FAIL):
fi:

S:={seq(L[ithprime(i)] mod ithprime(i),i=3..10)}:

if nops(S)<>1 then
 RETURN(FAIL):
fi:

a,S[1]:
end:






#Binet1(f,x,n): the explicit formula that handles the n-th coefficient of the Taylor expansion of f if it is a primaility test. 
#outputs as a set if pairs {[c,alpha]} and the expression is the sum of c*alpha^n
#Try:
#Binet1( (2-2*x)/(1-2*x-x^2),x,n):
Binet1:=proc(f,x) local f1,top,bot,gu,L,a,i,S,bot1,ka,mu,L1,n:

f1:=taylor(f,x=0,1001):
L:=[seq(coeff(f1,x,i),i=1..1000)]:

top:=numer(f):
bot:=denom(f):

gu:=[solve(bot,x)]:


bot1:=simplify(expand(mul(1-x/gu[i],i=1..nops(gu)))):


ka:=normal(bot1/bot):

if not type(ka,numeric) then
print(ka):

 RETURN(FAIL):
fi:

top:=ka*top:
bot:=bot1:


for i from 1 to nops(gu) do
mu[i]:=simplify(-1/gu[i]*subs(x=gu[i],top)/subs(x=gu[i],diff(bot,x))):
od:

a:=add(mu[i]/gu[i]^n,i=1..nops(gu)):

L1:=simplify([seq(subs(n=i,a),i=1..10)]):

if L1<>[op(1..10,L)] then
 RETURN(FAIL):
fi:

S:={seq(L[ithprime(i)] mod ithprime(i),i=3..10)}:

if nops(S)<>1 then
 RETURN(FAIL):
fi:

{seq([mu[i],1/gu[i]],i=1..nops(gu))},S[1]:
end:


#Find22(f,x): if a(n) has the form a1^n+a2^n then it finds an expression for a(2*n) in terms of other things. Try:
#Find22(PDB(x)[3][2],x);
Find22:=proc(f,x) local lu,c,a1,a2,i,n:
lu:=Binet1(f,x)[1]:

if nops(lu)<>2 or {seq(lu[i][1],i=1..nops(lu))}<>{1} then
RETURN(FAIL):
fi:

a1:=lu[1][2]:
a2:=lu[2][2]:
c:=simplify(a1*a2):

if not type(c,integer) then
RETURN(FAIL):
fi:
a(n)^2-2*c^n:
end:


#Find23(f,x): if a(n) has the form a1^n+a2^n then it finds an expression for a(3*n) in terms of other things. Try:
#Find23(PDB(x)[3][2],x);
Find23:=proc(f,x) local lu,c,a1,a2,i,n:
lu:=Binet1(f,x)[1]:

if nops(lu)<>2 or {seq(lu[i][1],i=1..nops(lu))}<>{1} then
RETURN(FAIL):
fi:

a1:=lu[1][2]:
a2:=lu[2][2]:
c:=simplify(a1*a2):

if not type(c,integer) then
RETURN(FAIL):
fi:
a(n)^3-3*c^n*a(n):
end:



#end proving infinite families








#FindPP2(f,x,K): All the pseudo-primes for the test inspired by f(x)  p*q, where p,q<=ithprime(K)  are  
#FindPP2( (x^2-2*x+3)/(-4*x^3+x^2-x+1),x,1000);
FindPP2:=proc(f,x,K) local i,j,p,q,gu,mu,p0,mu1,f1:


f1:=taylor(f,x=0,ithprime(K+10)+10):
mu:=[seq(coeff(f1,x,i),i=1..ithprime(K+10)+1)]:

mu1:={seq(mu[ithprime(i)] mod ithprime(i),i=10..20)}:



mu1:={seq(mu[ithprime(i)] mod ithprime(i),i=10..K+10)}:

if nops(mu1)<>1 then
 print(`Not a primality test`):
 RETURN(FAIL):
fi:
p0:=mu1[1]:

gu:=[]:

for i from 2 to K do
p:=ithprime(i):
for j from  i+1 to K do
q:=ithprime(j):
if mu[p]+mu[q] mod p*q=2*p0 then
 if not (p*q<=nops(mu) and mu[p*q] mod p*q<>p0) then
 gu:=[op(gu),p*q]:
 fi:
fi:
od:
od:

sort(gu):
end:




#FindPP3(f,x,K): All the pseudo-primes for the test inspired by f(x) of the form p*q*r, where p,q,r<=ithprime(K) 
#FindPP3( (-x^3-2*x^2+4)/(-x^4-x^3-x^2+1),x,1000);
FindPP3:=proc(f,x,K) local f1,i,j,p,q,gu,mu,k,r,mu1,p0:


f1:=taylor(f,x=0,ithprime(K+10)^2+10):

mu:=[seq(coeff(f1,x,i),i=1..ithprime(K+10)^2+1)]:

mu1:={seq(mu[ithprime(i)] mod ithprime(i),i=10..20)}:


mu1:={seq(mu[ithprime(i)] mod ithprime(i),i=10..20)}:



mu1:={seq(mu[ithprime(i)] mod ithprime(i),i=10..K+10)}:

if nops(mu1)<>1 then
 print(`Not a primality test`):
 RETURN(FAIL):
fi:
p0:=mu1[1]:

if p0<>0 then
 print(`Not yet implemented`):
 RETURN(FAIL):
fi:


gu:=[]:

for i from 1 to K do
p:=ithprime(i):
for j from  i+1 to K do
q:=ithprime(j):
for k from j+1 to K do
r:=ithprime(k):
if mu[p*q]+mu[q*r]+mu[p*r]-mu[p]-mu[q]-mu[r] mod p*q*r=0 then
gu:=[op(gu),p*q*r]:
 fi:

od:
od:
od:
sort(gu):
end:




#PPg(f,s,L): Memory efficient version of PPg(f,s,L). Try:
#PPg((1+s^2)/(1-s-s^2),x,1000);
PPg:=proc(f,s,L) local k,C,gu,CU,R,mu,INI,lu,i,i1,p0,i0:

p0:=coeff(taylor(f,s=0,4),s,1):


if p0<0 then
 print(`coeff. of`, s,  `in the denominator of f should be negative`):
 RETURN(FAIL):
fi:

for i from 1 while ithprime(i)<=p0 do od:
i0:=i:

C:=RtoC(f,s):

gu:=SeqFromRec(C,ithprime(i0+20)):



lu:={seq(gu[ithprime(i)] mod ithprime(i),i=i0..i0+20)} :

if lu<>{p0} then
 RETURN(FAIL):
fi:


INI:=C[1]:
R:=C[2]:
mu:=[]:


CU:=INI:

k:=nops(INI):

for i from k+1 to L do
 CU:=[op(2..k,CU),add(R[i1]*CU[k+1-i1],i1=1..k)]:

 
 if not isprime(i) and CU[k] mod i=p0 then
  mu:=[op(mu),i]:
 fi:
od:
mu:


end:




#FindPP2g(f,s,K): Given a rational function f of s and a list of two integers K=[K1,K2] finds
#all the pseudo-primes that are of the form ithprime(i1)*ithprime(i2) with i<=K1 and i2<=K2. Try:
#FindPP2g(GenPerrin([0,2,3,1],s),s,[10,20]);
FindPP2g:=proc(f,s,K) local L,C,gu,lu,i0,i1,p1,p2,i, i2, p0:

p0:=coeff(taylor(f,s=0,4),s,1):


if p0<0 then
 print(`coeff. of`, s,  `in the denominator of f should be negative`):
 RETURN(FAIL):
fi:

for i from 1 while ithprime(i)<=p0 do od:
i0:=i:

C:=RtoC(f,s):

gu:=SeqFromRec(C,ithprime(i0+20)):



lu:={seq(gu[ithprime(i)] mod ithprime(i),i=i0..i0+20)} :

if lu<>{p0} then
 RETURN(FAIL):
fi:



L:=mul(ithprime(K[i1]),i1=2..nops(K))+1:


gu:=SeqFromRec(C,L):
lu:=[]:
for i1 from 1 to K[1] do
  p1:=ithprime(i1):
 for i2 from i1+1 to K[2] do
  p2:=ithprime(i2):

  if (gu[p1]+gu[p2]-2*p0) mod p1*p2=0 then
    lu:=[op(lu),p1*p2]:
  fi:
 od:
od:
sort(lu):
end:



#FindPP3g(f,s,K): Given a rational function f of s and a list of three integers K=[K1,K2,K3] finds
#all the pseudo-primes that are of the form ithprime(i1)*ithprime(i2)*ithprime(i3) with i<=K1 and i2<=K2,i3<=K3. Try:
#FindPP3g(GenPerrin([0,2,3,1],s),s,[10,20,30]);
FindPP3g:=proc(f,s,K) local L,C,gu,lu,i0,i1,p1,p2,p3,i, i2, i3, p0:

p0:=coeff(taylor(f,s=0,4),s,1):


if p0<0 then
 print(`coeff. of`, s,  `in the denominator of f should be negative`):
 RETURN(FAIL):
fi:

for i from 1 while ithprime(i)<=p0 do od:
i0:=i:

C:=RtoC(f,s):

gu:=SeqFromRec(C,ithprime(i0+20)):



lu:={seq(gu[ithprime(i)] mod ithprime(i),i=i0..i0+20)} :

if lu<>{p0} then
 RETURN(FAIL):
fi:



L:=mul(ithprime(K[i1]),i1=2..nops(K))+1:

gu:=SeqFromRec(C,L):
lu:=[]:
for i1 from 1 to K[1] do
  p1:=ithprime(i1):
 for i2 from i1+1 to K[2] do
  p2:=ithprime(i2):
   for i3 from i2+1 to K[3] do
  p3:=ithprime(i3):
  if gu[p1*p2]+gu[p1*p3]+gu[p2*p3]-gu[p1]-gu[p2]-gu[p3] mod p1*p2*p3=0 then
    lu:=[op(lu),p1*p2*p3]:
  fi:
 od:
od:
od:
sort(lu):
end:




#FindPPg4(f,s,K): Given a rational function f of s and a list of four integers K=[K1,K2,K3,K4] finds
#all the pseudo-primes that are of the form ithprime(i1)*ithprime(i2)*ithprime(i3)*ithprime(i4) with i<=K1 and i2<=K2, i3<=K3, i4<=K4. Try:
#FindPP4g(GenPerrin([0,2,3,1],s),s,[5,10,10,10]);
FindPP4g:=proc(f,s,K) local L,C,gu,lu,i0,i1,p1,p2,p3,p4,i, i2, i3, i4,p0:

p0:=coeff(taylor(f,s=0,4),s,1):


if p0<0 then
 print(`coeff. of`, s,  `in the denominator of f should be negative`):
 RETURN(FAIL):
fi:

for i from 1 while ithprime(i)<=p0 do od:
i0:=i:

C:=RtoC(f,s):

gu:=SeqFromRec(C,ithprime(i0+20)):



lu:={seq(gu[ithprime(i)] mod ithprime(i),i=i0..i0+20)} :

if lu<>{p0} then
 RETURN(FAIL):
fi:



L:=mul(ithprime(K[i1]),i1=2..nops(K))+1:

gu:=SeqFromRec(C,L):
lu:=[]:
for i1 from 1 to K[1] do
  p1:=ithprime(i1):
 for i2 from i1+1 to K[2] do
  p2:=ithprime(i2):
   for i3 from i2+1 to K[3] do
  p3:=ithprime(i3):
  for i4 from i3+1 to K[4] do
  p4:=ithprime(i4):
  if 
gu[p1*p2*p3]+gu[p1*p2*p4]+gu[p1*p3*p4]+gu[p2*p3*p4]
-(gu[p1*p2]+gu[p1*p3]+gu[p1*p4]+gu[p2*p3]+gu[p2*p4]+gu[p3*p4])
+gu[p1]+gu[p2]+gu[p3]+gu[p4]
-2*p0 mod p1*p2*p3*p4 =0 then
    lu:=[op(lu),p1*p2*p3*p4]:
  fi:
 od:
od:
od:
od:

sort(lu):
end: