######################################################################
## PerrinVV.txt Save this file as  PerrinVV.  to use it,             #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read PerrinVV.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 PerrinVV.txt, one of the  Maple packages`):
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(`----------------------------`):
print(`It finds Perrin-style primality tests inspired by Vince Vatter's beautiful combinatorial proof`):
print(`of the original Perrin test, and the Edlin-Zeilberger version of the Goulden-Jackons algorithm`):
print(`----------------------------`):
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 Generalized functions type: ezraG();`):
print():
print(`For help with a specific procedure, type:`):
print(`"ezra(procedure_name);"`):
 print(`------------------------------------`):


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



print():
print(`For  help with procedures taken from CGJ `):
print(` type "ezraCGJ();". For specific help type "ezra(procedure_name);" `):
print(`For a list of the supporting functions type: ezraCGJ(procedure_name);`):
print():
 

with(linalg):
with(numtheory):


ezraG:=proc()
if args=NULL then

print(`The General procedures are`):
print(`   FindPP1G, GenPerrin3, GenPerrin, PPg,  PtestG,`):
print(`   GenPerrinPatternSearch, PatternDB`):

else
ezra(args):

fi:


end:

ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are`):
print(`  CheckSeq, FindPP1, FindPP2, FindPP2gen, FindPP3,  GenLucas, IncVecs, MatC, Pow1, PPlucas, RecPat, SeqLucas, ValC, ValCp, ValPatP `):

else
ezra(args):

fi:


end:



ezraS:=proc()
if args=NULL then

print(`The STORY procedures are`):
print(`  Info1, Paper1, Paper2, Paper3 `):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

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

print(`The MAIN procedures are`):
print(` FindPP, PP, Ptest, SeqPat, ValPat `):



elif nargs=1 and args[1]=CheckSeq then
print(`CheckSeq(L): Given a sequece of integers L,checks that the quantity`):
print(`Sum(L[n/d]*mobius(d), d in Divisors(n)). if true for n from 1 to nops(L). Try:`):
print(` CheckSeq(SeqPat({0,1},{[1,1]},100));`):

elif nargs=1 and args[1]=FindPP then
print(`FindPP(Alph,Pats,A,K): Inputs an alphabet Alph and a set of patterns Pats, and a template A=[a1,a2,...,ak] of positive integers , outputs`):
print(`of pseudo-primes for the primality tests induced by [Alph,Pats] of the form p1^a1*...*pk^ak with p1<..<pk primes among the first K primes and <=ithprime(K)^(sum(A)-1). Try:`):
print(`FindPP({0,1},{[1,1]},[1,1],30);`):

elif nargs=1 and args[1]=FindPP1 then
print(`FindPP1(Alph,Pats,k,K): All the pseudo-primes for the test inspired by Alph and Pats, of the form p^k, where p is a one of the first K primes.`):
print(`Try:`):
print(`FindPP1({0,1},{[1,1],[0,0,0]},2,1000);`):

elif nargs=1 and args[1]=FindPP1G then
print(`FindPP1G(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(`FindPP1G(CGJsG({0,1},{[1,1],[0,0,0]},s,0),s,2,1000);`):

elif nargs=1 and args[1]=FindPP2 then
print(`FindPP2(Alph,Pats,K): All the pseudo-primes for the test inspired by Alph and Pats, of the form p*q, where p and q are  one of the first K primes.`):
print(`Try:`):
print(`FindPP2({0,1},{[1,1],[0,0,0]},1000);`):

elif nargs=1 and args[1]=FindPP2gen then
print(`FindPP2gen(Alph,Pats,a1,a2,K): All the pseudo-primes for the test inspired by Alph and Pats, of the form p^a1*q^a2, where p,q<=ithprime(K)  are  `):
print(`FindPP2gen({0,1},{[1,1],[0,0,0]},1,1,1000);`):

elif nargs=1 and args[1]=FindPP3 then
print(`FindPP3(Alph,Pats,K): All the pseudo-primes for the test inspired by Alph and Pats, of the form p*q*r, where p,q,r<=ithprime(K) `):
print(`FindPP3({0,1},{[1,1],[0,0,0]},100);`):

elif nargs=1 and args[1]=GenLucas then
print(`GenLucas(x,y): The generating function of the General Lucas test`):

elif nargs=1 and args[1]=GenPerrin3 then
print(`GenPerrin3 (a1, a2, a3, x): The general Perrin-style rational function with cubic denominator. `):
print(`To get the original Perrin generating function.Try:`):
print(`GenPerrin3(0,-1,1,x);`):

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]=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]=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]=IncVecs then
print(`IncVecs(k,K): The set of all increasing vectors of length k with entries<=K, Try:`):
print(`IncVecs(4,5);`):


elif nargs=1 and args[1]=Info1 then
print(`Info1(Alph,Pats,K1,K2,K3,K4): Given an alphabet Alph and a set of patterns Pats, and positive integers K1,K2,K3, outputs`):
print(`* The generating function in s for the sequence enumerating circular words in the alphabet Alph, avoiding  the patterns in Pats (that induce this primality test)`):
print(`*for the sake of OEIS, the first K1 terms of the sequence`):
print(`*All the pseudo-primes <=K2`):
print(`*All the presudoprimes that are perfect squares of prime  <=ithprime(K3)`):
print(`*All the presudoprimes that are  products of two distinct primes <=ithprime(K3). For the Lucas test, Try:`):
print(`Info1({0,1},{[1,1]},30,10000,1000,100);`):

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]=Paper1 then
print(`Paper1(A,k1,s,K): Finds all the primality tests inspired by an alphabet of size from 2 to A with set of forbidden patterns (all constant words) of size 1 where the patterns have length from 2 to to k1`):
print(`that are constant words, and all their pseudo-primes <=K, it also determines the one with the smallest number of pseu-primes <=K. Try:`):
print(`Paper1(2,5,s,1000);`):


elif nargs=1 and args[1]=Paper2 then
print(`Paper2(A,k1,s,K): Finds all the primality tests inspired by an alphabet of size from 2 to A with set of forbidden patterns (all constant words) of size 2 where the patterns have length from 2 to to k1`):
print(`and the patterns are constant words (i.e. only have one letter) and all their pseudo-primes <=K, it also determines the one with the smallest number of pseu-primes <=K. Try:`):
print(`Paper2(2,5,s,1000);`):



elif nargs=1 and args[1]=Paper3 then
print(`Paper3(A,k1,s,K): Finds all the primality tests inspired by an alphabet of size from 3 to A with set of forbidden patterns (all constant words) of size 3 where the patterns have length from 2 to to k1`):
print(`and the patterns are constant words (i.e. only have one letter) and all their pseudo-primes <=K, it also determines the one with the smallest number of pseu-primes <=K. Try:`):
print(`Paper3(3,5,s,1000);`):


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]=PP then
print(`PP(Alph,Pats,L): All the pseudo-primes <=L in the primality set induced by the alphabet Alph and the set of patterns Pats. Try:`):
print(`PP({0,1},{[1,1],[0,0,0]},1000);`):

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((1+s^2)/(1-s-s^2),s,1000);`):

elif nargs=1 and args[1]=PPlucas then
print(`PPlucas(y0,L): All the pseudo-primes <=L in the generalized Lucas test with y=y0. Try:`):
print(`PPlucas(1,1000);`):

elif nargs=1 and args[1]=Ptest then
print(`Ptest(Alph,Pats,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 set of patterns Pats in the alphabet Alph. For the Perrin test try:`):
print(`Ptest({0,1},{[1,1],[0,0,0]},1001,0); For the Lucas test try:`):
print(`Ptest({0,1},{[1,1]},1001,1); `):

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(CGJs({0,1},{[1,1],[0,0,0]},s),s,1001,0);`):

elif nargs=1 and args[1]=SeqLucas then
print(`SeqLucas(y0,L): Given an alphabet Alph, a finite set of patterns Pats, and a positive integer L, outputs the first L terms of`):
print(`the number of the generating function for the generalized Lucas numbers with y=y0.`):
print(`SeqLucas(1,00);`):

elif nargs=1 and args[1]=RecPat then
print(`RecPat(Alph,Pats): The C-finite recurrence satisfied by the set of circular words in the alphabet Alph avoiding the set of patterns Pats. For the Perrim sequence, type:`):
print(`RecPat({0,1},{[1,1],[0,0,0]});`):


elif nargs=1 and args[1]=SeqPat then
print(`SeqPat(Alph,Pats,L): Given an alphabet Alph, a finite set of patterns Pats, and a positive integer L, outputs the first L terms of`):
print(`the number of circular words (with  a clasp) in the alphabet Alph avoding the set of Patterns Pats. For the first 100 terms of the Perrin sequence (except the begining) try:`):
print(`SeqPat({0,1},{[1,1],[0,0,0]},100);`):

elif nargs=1 and args[1]=ValC then
print(`ValC(C,n): The n-th term of the C-finite sequence given by the recurrence C=[INI,REC], using Matrix multiplication. Should be the same as`):
print(`SeqFromRec(C,n)[n], but hopefully faster. Try:`):
print(`ValC([[1,1],[1,1]],30);`):

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);`):

elif nargs=1 and args[1]=ValPat then
print(`ValPat(Alph,Pats,n): The number of circular words of length n in the alphabet Alph avoiding as consecutive subwords the patterns in the set Pats`):
print(`Same, but faster as SeqPat(Alph,Pats,n)[n]. Try:`):
print(`ValPat({0,1},{[1,1],[0,0,0]},1000);`):


elif nargs=1 and args[1]=ValPatP then
print(`ValPatP(Alph,Pats,n,p): The number of circular words of length n in the alphabet Alph avoiding as consecutive subwords the patterns in the set Pats mod p`):
print(`Same, but faster as ValPat(Alph,Pats,n) mod p. Try:`):
print(`ValPatP({0,1},{[1,1],[0,0,0]},1000,1003);`):


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

fi:


end:
 
 

## START FROM Cfinite.txt
#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:


##END FROM Cfinite.txt

###FROM CGJ
ezraCGJ:=proc()
if args=NULL then
print(`CGJ`):
 
print(` Contains Procedures: `):
print(` bdokC, CGJs , CGJsG, ProveBW, SetOfCorrectWordsC,WrapAround `):
print(``):
print(` For specific help type "ezra(procedure_name)" `):
fi:
 
 
if nops([args])=1 and op(1,[args])=bdokC then
print(`bdokC(alphabet,MISTAKES,L): checks the validity `):
print(`of the Edlin-Zeilberger method for the given alphabet`):
print(` and set of MISTAKES up to words of length L`):
print(` For example bdokC({1,2},{[1,1],[2,2]},10); `):
fi:
 
if nops([args])=1 and op(1,[args])=ProveBW then
print(`ProveBW(): proves the Burstein-Wilf formula`):
fi:
 
if nops([args])=1 and op(1,[args])=WrapAround then
print(`WrapAround(MISTAKES,s): The contibution from Case III clusters`):
print(`in the Edlin-Zeilberger cyclic words analog of Goulden-Jackosn`):
print(`index by the LIST of MISTAKES, and whose entries give all possible`):
print(`overlaps, phrased in terms  of the variable s`):
fi:
 
 
if nops([args])=1 and op(1,[args])=`CGJs` then
print(`CGJs(alphabet,MISTAKES,s) finds the generating function`):
print(`g(s), is the g.f. sum(a[i]*s^i,i=0..infty)  `):
print(`where a[i]:=number of cyclic words of length i without any mistakes`):
print(`from MISTAKES`):
print(`For example to get the g.f. for a_n:=number of cyclic `):
print(` words of length n`):
print(`in the alphabet {1,2,3}, avoiding 123,231,312, type`):
print(` type CGJs({1,2,3},{[1,2,3],[2,3,1],[3,1,2]},s)`):
fi:
 

if nops([args])=1 and op(1,[args])=`CGJsG` then
print(`CGJsG(alphabet,MISTAKES,s,t) finds the generating function`):
print(`in s of the sequence of weight-enumerators  of the words of length n according to the weight t^NumverOfOccurrencesOfAmistake `):
print(`For example to get the g.f. for a_n:=number of cyclic `):
print(` words of length n`):
print(`in the alphabet {1,2,3}, avoiding 123,231,312, type`):
print(` type CGJsG({1,2,3},{[1,2,3],[2,3,1],[3,1,2]},s,t)`):
fi:
 
 
 
end:
 
with(linalg):
CGJs:=proc(alphabet,MISTAKES1,s)
local v,eq, var,i,lu,C,MISTAKES,mu:
 
MISTAKES:=Hakten(MISTAKES1):
eq:={}:
var:={}:
 
for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 eq:= eq union {findeqz(v,MISTAKES,C,s)}:
 var:=var union {C[op(v)]}:
od:
 
var:=solve(eq,var):
 
lu:=1:
 mu:=0:
for i from 1 to nops(alphabet) do
 lu:=lu-s:
od:
 
for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 lu:=lu-subs(var,C[op(v)]):
 mu:=mu+subs(var,C[op(v)]):
od:
 
normal((1/lu)*(1+s*diff(mu,s)-mu)+WrapAround(convert(MISTAKES,list),s)):
 
end:
 
 
#findeqz sets up the equ C[v]= s+t*Sum_u overlap(u,v,x) *C[u]
 
findeqz:=proc(v,MISTAKES,C,s)
local eq,i,u:
eq:=-1:
 
for i from 1 to nops(v) do
 eq:=eq*s:
od:
 
for i from 1 to nops(MISTAKES) do
 u:=op(i,MISTAKES):
 eq:=eq-overlapz(u,v,s)*C[op(u)]:
od:
 
C[op(v)]-eq=0:
 
end:
 
 
 
#overlapz is a procedure that given two words u and v, and a variable s
#computes the weight-enumerator of all v\suffix(u),
#for all suffixes of u that are prefixes of v, but with uniform weight s
overlapz:=proc(u,v,s)
local i,j,lu,gug:
 
lu:=0:
 
for i from 2 to nops(u) do
  for j from i to nops(u) while (j-i+1<=nops(v) and op(j,u)=op(j-i+1,v))
         do
   :
  od:   
 
if j-i=nops(v) and u<>v then
 ERROR(v,`is a subword of`,u,`illegal input`):
 fi:
 
 
 
 if j=nops(u)+1 and (i>1 or j>2) then
  gug:=1:
 
   gug:=gug*s^(nops(v)-(j-i)):
 
  lu:=lu+gug:
 
fi:
 
od:
 
lu:
 
end:
 
 
 
#WrapAround(MISTAKES,s): The contibution from Case I clusters
#in the Edlin-Zeilberger cyclic words analog of Goulden-Jackosn
#index by the LIST of MISTAKES, and whose entries give all possible
#overlaps, phrased in terms  of the variable s
WrapAround:=proc(MISTAKES,s)
local n,A,B,C,i,j,D1,lu,gu,gu1,r:
with(linalg):
n:=nops(MISTAKES):
A:=matrix(n,n):
B:=matrix(n,n):
C:=matrix(n,n):
n:=nops(MISTAKES):
for i from 1 to n do
for j from 1 to n do
A[i,j]:=-overlapz(MISTAKES[i],MISTAKES[j],s):
B[i,j]:=s*diff(A[i,j],s):
if i=j then
C[i,j]:=1-A[i,j]:
else
C[i,j]:=-A[i,j]:
fi:
od:
od:
D1:=multiply(multiply(A,inverse(C)),B):
gu:=0:
for i from 1 to n do
gu1:=D1[i,i]:
lu:=taylor(gu1,s=0,nops(MISTAKES[i])+1):
for r from 0 to nops(MISTAKES[i])-1 do
gu1:=gu1-coeff(lu,s,r)*s^r:
od:
gu:=normal(gu+gu1):
od:
gu:
end:
 
 
 
 
 
hakten:=proc(B)
local w,i:
 
for i from 1 to nops(B) do
w:=op(i,B):
 
if superflous(B,w)=1 then
 RETURN(B minus {w}):
fi:
od:
B:
end:
 
#Hakten(B) removes all superflous words
 
Hakten:=proc(B)
local B1,B2:
 
B1:=B:
 
B2:=hakten(B1):
 
while B1<>B2 do
B1:=B2:
B2:=hakten(B2):
od:
 
B2:
end:
 
 
#BW(k,w,x): The Burstein-Wilf generating function
BW:=proc(k,w,x)
normal(1+
(1-x^w)/(1-x)*
(k*x+(k-1)*x*(
                 (w+1-w*k*x)/(1-k*x+(k-1)*x^(w+1))- (w+1)/(1-x^(w+1))
             )
)):
end:
 
 
#BWb(k,w): Checking the Edlin-Zeilberger Analog of Goulden-Jackson
#for cyclic words against the Burstein-Wilf generating function
BWb:=proc(k,w)
local al,MI,i,j,x,gu:
al:={seq(i,i=1..k)}:
 
MI:={seq([seq(j,i=1..w+1)],j=1..k)}:
 
gu:=CGJs(al,MI,x):
gu:=gu-BW(k,w,x):
normal(gu):
 
end:
 
#bdokC(alphabet,MISTAKES,L): checks the validity of
#of the Edlin-Zeilberger method
bdokC:=proc(alphabet,MISTAKES,L)
local lu,gu,s,i:
 
lu:=EmpiricalGJseriesC(alphabet,MISTAKES,L):
gu:=taylor(CGJs(alphabet,MISTAKES,s),s=0,L+1):
gu:=[seq(coeff(gu,s,i),i=0..L)]:
evalb(lu=gu):
end:
 
#EmpiricalGJSeriesC(alphabet,MISTAKES1,L): Does what
#GJseriesC(nops(alphabet),MISTAKES,L) does, but
#empirically. Don't make L too large, unless there 
#are lots of mistakes. It is meant to be a check
#on the validity of the (usually) much much faster
#and should not be used except for curiosity
EmpiricalGJseriesC:=proc(alphabet,MISTAKES,L)
local n,gu:
gu:=[]:
 
for n from 0 to L do
gu:=[op(gu),nops(SetOfCorrectWordsC(alphabet,MISTAKES,n))]:
od:
gu:
end:
 
 
#SetOfCorrectWords(Alphabet,MISTAKES,n): given an alphabet, Alphabet
#and a set of mistakes, MISTAKES, finds, empirically, the
#set of words that that do not have factors that belong to MISTAKES
SetOfCorrectWords:=proc(Alphabet,MISTAKES,n)
local MISTAKES1,gu,mu,i,j,cand:
option remember:
MISTAKES1:=Hakten(MISTAKES):
 
if n=0 then
RETURN({[]}):
fi:
 
mu:=SetOfCorrectWords(Alphabet,MISTAKES,n-1):
 
gu:={}:
 
for i from 1 to nops(mu) do
 for j from 1 to nops(Alphabet) do
   cand:=[op(op(i,mu)),op(j,Alphabet)]:
 
     if not isbad(cand,MISTAKES1) then
          gu:=gu union {cand}:
     fi:
 od:
od:
 
gu:
 
end:
 
 
#SetOfCorrectWordsC(Alphabet,MISTAKES,n): given an alphabet, Alphabet
#and a set of mistakes, MISTAKES, finds, empirically, the
#set of CYCLIC words (with a marked beginning)
#that that do not have factors that belong to MISTAKES
SetOfCorrectWordsC:=proc(Alphabet,MISTAKES,n)
local gu,mila,mu,i:
mu:=SetOfCorrectWords(Alphabet,MISTAKES,n):
gu:={}:
 
for i from 1 to nops(mu) do
mila:=op(i,mu):
if not isbadC(mila,MISTAKES) then
 gu:=gu union {mila}:
fi:
od:
gu:
end:
 
 
 
#isbadC(mila,MISTAKES): returns true if one of the suffices of
#one of the CYCLIC shifts of mila belongs to the set of MISTAKES
isbadC:=proc(mila,MISTAKES)
local i,mila1,n:
n:=nops(mila):
for i from 1 to n do
mila1:=[op(i..n,mila),op(1..i-1,mila)]:
 
if isbad(mila1,MISTAKES) then
RETURN(true):
fi:
od:
 
false:
end:
 
 
#isbad(mila,MISTAKES): returns true if one of the suffices of
#mila belongs to the set of MISTAKES
isbad:=proc(mila,MISTAKES)
local i,shegi:
 
for i from 1 to nops(MISTAKES) do
shegi:=op(i,MISTAKES):
 
if nops(mila)>=nops(shegi) and 
 [op(nops(mila)-nops(shegi)+1..nops(mila),mila)]=shegi then
 RETURN(true):
fi:
 
od:
 
false:
 
end:
 
 
SBW:=proc(k,w,s)
local i,L,gu,mu,A,matkhil:
A:=-sum(s^i,i=1..w):
L:=-k*s^(w+1)/(1-A):
matkhil:=sum((i-1)*s^i,i=1..w):
gu:=k*normal(s*diff(A,s)*A/(1-A))-k*matkhil:
mu:=normal(1/(1-k*s-L)):
normal((s*diff(L,s)-L+1)*mu+gu):
end:
 
#ProveBW(): proves the Burstein-Wilf formula
ProveBW:=proc()
local i,L,gu,mu,A,matkhil,k,w,s:
A:=-sum(s^i,i=1..w):
L:=-k*s^(w+1)/(1-A):
matkhil:=sum((i-1)*s^i,i=1..w):
gu:=k*normal(s*diff(A,s)*A/(1-A))-k*matkhil:
mu:=normal(1/(1-k*s-L)):
evalb(simplify(normal((s*diff(L,s)-L+1)*mu+gu-BW(k,w,s)))=0):
end:
 

findeqzG:=proc(v,MISTAKES,C,s,t)
local eq,i,u:
eq:=-1:
 
for i from 1 to nops(v) do
 eq:=eq*s:
od:
 
for i from 1 to nops(MISTAKES) do
 u:=op(i,MISTAKES):
 eq:=eq+(t-1)*overlapz(u,v,s)*C[op(u)]:
od:
 
C[op(v)]-eq=0:
 
end:

CGJsG:=proc(alphabet,MISTAKES1,s,t)
local v,eq, var,i,lu,C,MISTAKES,mu:
 
MISTAKES:=Hakten(MISTAKES1):
eq:={}:
var:={}:
 
for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 eq:= eq union {findeqzG(v,MISTAKES,C,s,t)}:
 var:=var union {C[op(v)]}:
od:
 
var:=solve(eq,var):
 
lu:=1:
 mu:=0:
for i from 1 to nops(alphabet) do
 lu:=lu-s:
od:
 
for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 lu:=lu-subs(var,C[op(v)]):
 mu:=mu+subs(var,C[op(v)]):
od:
 
normal((1/lu)*(1+s*diff(mu,s)-mu)+WrapAround(convert(MISTAKES,list),s)):
 
end:

###END FROM CGJ


#SeqPatOri(Alph,Pats,L): Given an alphabet Alph, a finite set of patterns Pats, and a positive integer L, outputs the first L terms of
#the number of circular words (with  a clasp) in the alphabet Alph avoding the set of Patterns Pats. For the first 100 terms of the Perrin sequence (except the begining) try:
#SeqPatOri({0,1},{[1,1],[0,0,0]},100);
SeqPatOri:=proc(Alph,Pats,L) local f,s,f1,i:

f:=CGJs(Alph,Pats,s):
f1:=taylor(f,s=0,L+5):
[seq(coeff(f1,s,i),i=1..L)]:
end:


#SeqPat(Alph,Pats,L): Given an alphabet Alph, a finite set of patterns Pats, and a positive integer L, outputs the first L terms of
#the number of circular words (with  a clasp) in the alphabet Alph avoding the set of Patterns Pats. For the first 100 terms of the Perrin sequence (except the begining) try:
#SeqPat({0,1},{[1,1],[0,0,0]},100);
SeqPat:=proc(Alph,Pats,L) local f,s,f1,i,D1:

f:=CGJs(Alph,Pats,s):

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

end:



#SeqPatSlow(Alph,Pats,L): Given an alphabet Alph, a finite set of patterns Pats, and a positive integer L, outputs the first L terms of
#the number of circular words (with  a clasp) in the alphabet Alph avoding the set of Patterns Pats. For the first 100 terms of the Perrin sequence (except the begining) try:
#SeqPatSlow({0,1},{[1,1],[0,0,0]},100);
SeqPatSlow:=proc(Alph,Pats,L) local f,s,C:
f:=CGJs(Alph,Pats,s):
C:=RtoC(f,s):
SeqFromRec(C,L):
end:



#PP(Alph,Pats,L): All the pseudo-primes <=L in the primality set induced by the alphabet Alph and the set of patterns Pats. Try:
#PP({0,1},{[1,1],[0,0,0]},1000);
PP:=proc(Alph,Pats,L) local gu,i,mu,lu:
if L<=100 then
 print(` Make `, L, `bigger `):
 RETURN(FAIL):
fi:

gu:=SeqPat(Alph,Pats,L):

lu:={seq(gu[ithprime(i)] mod ithprime(i),i=4..10)}:

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:


#ValC(C,n): The n-th term of the C-finite sequence given by the recurrence C=[INI,REC], using Matrix multiplication. Should be the same as
#SeqFromRec(C,n)[n], but hopefully faster. Try:
#ValC([[1,1],[1,1]],30);
ValC:=proc(C,n) 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(M&^(n-k),Ini)[1,1]:
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:



#RecPat(Alph,Pats): The C-finite recurrence satisfied by the set of circular words in the alphabet Alph avoiding the set of patterns Pats. For the Perrim sequence, type:
#RecPat({0,1},{[1,1],[0,0,0]});
RecPat:=proc(Alph,Pats) local gu,s:
option remember:
gu:=CGJs(Alph,Pats,s):
RtoC(gu,s):
end:


#ValPat(Alph,Pats,n): The number of circular words of length n in the alphabet Alph avoiding as consecutive subwords the patterns in the set Pats
#Same, but faster as SeqPat(Alph,Pats,n)[n]. Try:
#ValPat({0,1},{[1,1],[0,0,0]},1000);
ValPat:=proc(Alph,Pats,n) ValC(RecPat(Alph,Pats),n):end:



#ValPatP(Alph,Pats,n,p): The number of circular words of length n in the alphabet Alph avoiding as consecutive subwords the patterns in the set Pats mod p
#Same, but faster as ValPat(Alph,Pats,n) mod p. Try:
#ValPatP({0,1},{[1,1],[0,0,0]},1000,101);
ValPatP:=proc(Alph,Pats,n,p) ValCp(RecPat(Alph,Pats),n,p):end:


#Ptest(Alph,Pats,n,p0): Testing pos. integer n according to the primality test where the p-th term mod p should be p0,
#inspired by the set of patterns Pats in the alphabet Alph. For the Perrin test try:
#Ptest({0,1},{[1,1],[0,0,0]},1001,0); For the Lucas test try:
#Ptest({0,1},{[1,1]},1001,1); 
Ptest:=proc(Alph,Pats,n,p0)
evalb(ValPatP(Alph,Pats,n,n)=p0):
end:



#FindPP1slow(Alph,Pats,k,K): All the pseudo-primes for the test inspired by Alph and Pats, of the form p^k, where p is a one of the first K primes.
#Try:
#FindPP1slow({0,1},{[1,1],[0,0,0]},2,1000);
FindPP1slow:=proc(Alph,Pats,k,K) local i,p,gu:

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



#FindPP1(Alph,Pats,k,K): All the pseudo-primes for the test inspired by Alph and Pats, that happens to be L[n],of the form p^k, where p is a one of the first K primes.
#Try:
#FindPP1({0,1},{[1,1],[0,0,0]},2,1000);
FindPP1:=proc(Alph,Pats,k,K) local i,p,gu,mu,mu1,p0:

mu:=SeqPat(Alph,Pats,ithprime(K)^(k-1)+10):



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:



#FindPP2(Alph,Pats,K): All the pseudo-primes for the test inspired by Alph and Pats, of the form p*q, where p,q<=ithprime(K)  are  
#FindPP2({0,1},{[1,1],[0,0,0]},1000);
FindPP2:=proc(Alph,Pats,K) local i,j,p,q,gu,mu,p0,mu1:


mu:=SeqPat(Alph,Pats,ithprime(K+10)+1):

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(Alph,Pats,K): All the pseudo-primes for the test inspired by Alph and Pats, of the form p*q*r, where p,q,r<=ithprime(K) 
#FindPP3({0,1},{[1,1],[0,0,0]},1000);
FindPP3:=proc(Alph,Pats,K) local i,j,p,q,gu,mu,k,r:


mu:=SeqPat(Alph,Pats,ithprime(K)^2+1):


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:



#CheckSeq(L,n): Given a sequence of integers L, and a positive integer n, checks that the
#Sum(L[n/d]*mobius(d), d in Divisors(n)) is always 0. Try:
#L:=SeqPat({0,1},{[1,1]},1000): CheckSeq(L);
CheckSeq:=proc(L) local D1,d1,i,n:
if not (type(L,list) and {seq(type(L[i],integer),i=1..nops(L))}={true}) then
 print(`bad input`):
RETURN(FAIL):
fi:
for n from 1 to nops(L) do

D1:=divisors(n):
 if add(L[n/d1]*mobius(d1), d1 in D1) mod n<>0 then
  print(`does not work for n=`, n):
  RETURN(false):
 fi:

od:
true:
end:



#IncVecs(k,K): The set of all increasing vectors of length k with entries<=K, Try
#IncVecs(4,5);
IncVecs:=proc(k,K) local gu,mu,i,mu1:
option remember:
if k=1 then
 RETURN({seq([i],i=1..K)}):
fi:
gu:={}:

for i from 1 to K do
 mu:=IncVecs(k-1,i-1):
 gu:=gu union {seq([op(mu1),i],mu1 in mu)}:
od:
gu:
end:



#FindPP(Alph,Pats,A,K): Inputs an alphabet Alph and a set of patterns Pats, and a template A=[a1,a2,...,ak] of positive integers , outputs
#of pseudo-primes for the primality tests induced by [Alph,Pats] of the form p1^a1*...*pk^ak with p1<..<pk primes among the first K primes and <=ithprime(K)^(sum(A)-1). Try:
#FindPP({0,1},{[1,1]},[1,1],30);
FindPP:=proc(Alph,Pats,A,K) local mu,S,k,gu,n,s,N,D1,d1,i:
N:=ithprime(K)^(convert(A,`+`)-1):
mu:=SeqPat(Alph,Pats,N):

k:=nops(A):
S:=IncVecs(k,K):
gu:=[]:

for s in S do

 n:=mul(ithprime(s[i])^A[i],i=1..nops(A)):

D1:=divisors(n) minus {n}:

if max(D1)<=N then
 if (add(mu[d1]*mobius(n/d1), d1 in D1)+mu[1]) mod n=0 then
  gu:=[op(gu),n]:
 fi:
fi:
od:
sort(gu):

end:



#FindPP2gen(Alph,Pats,a1,a2,K): All the pseudo-primes for the test inspired by Alph and Pats, of the form p^a1*q^a2, where p,q<=ithprime(K)  are  
#FindPP2gen({0,1},{[1,1],[0,0,0]},1,1,1000);
FindPP2gen:=proc(Alph,Pats,a1,a2,K) local i,j,p,q,gu,mu,N:


N:=ithprime(K)^(a1+a2-1)+10:
mu:=SeqPat(Alph,Pats,N):


gu:=[]:

for i from 1 to K do
p:=ithprime(i):
for j from  i+1 to K do
q:=ithprime(j):
if max(p^(a1-1)*q^a2, p^(a1-1)*q^a2)<=N then
if mu[p^(a1-1)*q^a2]+mu[p^a1*q^(a2-1)]-mu[p^(a1-1)*q^(a2-1)]-mu[1] mod p^a1*q^a2=0 then
 gu:=[op(gu),p^a1*q^a2]:
fi:
fi:
od:
od:

sort(gu):
end:





#Info1(Alph,Pats,K1,K2,K3,K4): Given an alphabet Alph and a set of patterns Pats, a symbol s, and positive integers K1,K2,K3, outputs
#* The generating function in s for the sequence enumerating circular words in the alphabet Alph, avoiding  the patterns in Pats (that induce this primality test)
#*for the sake of OEIS, the first K1 terms of the sequence
#*All the pseudo-primes <=K2
#*All the presudoprimes that are perfect squares of prime  <=ithprime(K3)
#*All the presudoprimes that are  products of two distinct primes <=ithprime(K3). For the Lucas test, Try:
#Info1({0,1},{[1,1]},s,30,10000,1000,100);

Info1:=proc(Alph,Pats,K1,K2,K3,K4) local s,f,D1,L,n,f1,f1a,C,a,i,j,q,lu,lu1,t0,lu2:
t0:=time():
f:=CGJs(Alph,Pats,s):

D1:=denom(f):
f1:=rem(numer(f),D1,s):
f1:=f1/D1:

f1a:=taylor(f1,s=0,K2+5):
L:=[seq(coeff(f1a,s,i),i=1..K2)]:

print(`A Primality Test inspired by circular words in`, Alph, `avoiding Consecutive subwords in`, Pats, `and their first few Pseudoprimes`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`The following primality test is generated by the sequence enumerating  words of length n in the alphabet`, Alph, `avoiding as consecutive words the patterns in`, Pats):
print(``):
print(`Theorem: Let a(i) be the sequence given by the generating function`):
print(``):
print(Sum(a(n)*s^n,n=0..infinity)= f1):
print(``):
C:=RtoC(f,s):
print(`Or equivalenty the sequence satisfying the recurrence`):
print(``):
print(a(n)=add(C[2][j]*a(n-j),j=1..nops(C[2]))):
print(``):
print(`with initial conditions `):
print(``):
print(seq(a[j]=C[1][j],j=1..nops(C[1]))):
print(``):
print(`then if p is a prime, then a(p) is always congruent mod p to`, L[1]):
print(``):
print(`Proof: Using the Edlin-Zeilberger variation of the Goulden-Jackson method,  Adv. in Appl. Math 25(2000), 228-232,  https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/cgj.html`):
print(`the above generating function (except for the first few terms) enumerates circular words (with a clasp)`):
print(`Obviously unless  all the letters are the same, all the p rotations are distinct, hence this set of non-constant words can be partitioned into equivalence classes under rotation, each of them`):
print(`with p members`):
print(``):
print(`---------------------`):
print(``):

print(``):
print(`For the sake of the OEIS the first`, min(K1,K2), `terms, starting at n=1 are `):
print(``):
print(op(1..min(K1,K2), L)):
print(``):
lu:=PP(Alph,Pats,K2):
print(``):
print(`There are`, nops(lu), `pseduo-primes <=`, K2, `here there are `):
print(``):
print(lu):
print(``):
lu1:=FindPP1(Alph,Pats,2,K3):
print(``):
print(`There are`, nops(lu1), `pseduo-primes that are squares of primes, among the first`, K3, `primes. Here there are  `):
print(lu1):

print(``):
lu2:=FindPP2(Alph,Pats,K4):
print(``):
print(`There are`, nops(lu2), `pseduo-primes that are products of two distinct primes, among the first`, K4, `primes. Here there are  `):
print(lu2):
print(``):

print(`-------------------------------`):
print(``):
print(`This ends this article that took`, time()-t0, `seconds. `):
print(``):
[lu,lu1,lu2]:

end:



#Paper1(A,k1,s,K): Finds all the primality tests inspired by an alphabet of size from 2 to A, with set of forbidden patterns (all constant words) of size 1 where the patterns have length from 2 to to k1
#that are constant words and all their pseudo-primes <=K, it also determines the one with the smallest number of pseudo-primes <=K. Try:
#Paper1(2,5,s,1000);
Paper1:=proc(A,k1,s,K) local a,gu, aluf, si,Alph,Pats,f,lu,t0,D1,f1,i,mu:

print(`All the Primality Tests Inspired by an alphabet of 2 to `, A, `letters and sets of forbidden words of size 1 where the words are constants and of length from 2 to`, k1):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):

t0:=time():
si:=K:

for a from 2 to A do

Alph:={seq(i,i=1..a)}:

gu:=[]:

for i from 2 to k1 do
 Pats:={[1$i]}:

f:=CGJs(Alph,Pats,s):

D1:=denom(f):
f1:=rem(numer(f),D1,s):
f1:=f1/D1:

lu:=PP(Alph,Pats,K):

mu:=[[Alph,Pats],f1, lu]:

if nops(lu)<si then
si:=nops(lu):
aluf:=mu:
fi:


lprint(mu):
print(``):
gu:=[op(gu),mu]:
od:
od:
print(`The set of patterns with, followed by the generating function, with the least number of pseudoprimes less than`, K , `is: `):
print(``):
lprint(aluf):
print(``):
print(`---------------------`):
print(``):
print(`This took`, time()-t0, `seconds. `):
print(``):
end:





#Paper2(A,k1,s,K): Finds all the primality tests inspired by an alphabet of size from 2 to A, with set of forbidden patterns (all constant words) of size 2 where the patterns have length from 2 to to k1
#that are constant words and all their pseudo-primes <=K, it also determines the one with the smallest number of pseudo-primes <=K. Try:
#Paper2(2,5,s,1000);
Paper2:=proc(A,k1,s,K) local a,gu, aluf, si,Alph,Pats,f,lu,t0,D1,f1,i1,i2,mu,i:

print(`All the Primality Tests Inspired by an alphabet of 2 to `, A, `letters and sets of forbidden words of size 1 where the words are constants and of length from 2 to`, k1):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):

t0:=time():
si:=K:

for a from 2 to A do

Alph:={seq(i,i=1..a)}:

gu:=[]:

for i1 from 2 to k1 do
 for i2 from i1+1 to k1 do
 Pats:={[1$i1],[2$i2]}:

f:=CGJs(Alph,Pats,s):

D1:=denom(f):
f1:=rem(numer(f),D1,s):
f1:=f1/D1:

lu:=PP(Alph,Pats,K):

mu:=[[Alph,Pats],f1, lu]:

if nops(lu)<si then
si:=nops(lu):
aluf:=mu:
fi:


lprint(mu):
print(``):
gu:=[op(gu),mu]:
od:
od:
od:

print(`The set of patterns with, followed by the generating function, with the least number of pseudoprimes less than`, K , `is: `):
print(``):
lprint(aluf):
print(``):
print(`---------------------`):
print(``):
print(`This took`, time()-t0, `seconds. `):
print(``):
end:




#Paper3(A,k1,s,K): Finds all the primality tests inspired by an alphabet of size from 3 to A, with set of forbidden patterns (all constant words) of size 3 where the patterns have length from 2 to to k1
#that are constant words and all their pseudo-primes <=K, it also determines the one with the smallest number of pseudo-primes <=K. Try:
#Paper3(2,5,s,1000);
Paper3:=proc(A,k1,s,K) local a,gu, aluf, si,Alph,Pats,f,lu,t0,D1,f1,i1,i2,i3,mu,i:

if A<3 then
 print(`A must be at least 3`):
 RETURN(FAIL):
fi:

print(`All the Primality Tests Inspired by an alphabet of 2 to `, A, `letters and sets of forbidden words of size 1 where the words are constants and of length from 2 to`, k1):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):

t0:=time():
si:=K:

for a from 3 to A do

Alph:={seq(i,i=1..a)}:

gu:=[]:

for i1 from 2 to k1 do
 for i2 from i1+1 to k1 do
  for i3 from i2+1 to k1 do
 Pats:={[1$i1],[2$i2],[3$i3]}:

f:=CGJs(Alph,Pats,s):

D1:=denom(f):
f1:=rem(numer(f),D1,s):
f1:=f1/D1:

lu:=PP(Alph,Pats,K):

mu:=[[Alph,Pats],f1, lu]:

if nops(lu)<si then
si:=nops(lu):
aluf:=mu:
fi:


lprint(mu):
print(``):
gu:=[op(gu),mu]:
od:
od:
od:
od:

print(`The set of patterns with, followed by the generating function, with the least number of pseudoprimes less than`, K , `is: `):
print(``):
lprint(aluf):
print(``):
print(`---------------------`):
print(``):
print(`This took`, time()-t0, `seconds. `):
print(``):
end:






#GenLucas(x,y): The generating function of the General Lucas test
GenLucas:=proc(x,y)
normal((1+(1-y)*x+x^2*y^2+(1-y)*x^3*y-(1-y)^2*x^4*y)/((1-x*y)*(1-x^2*y)-x^3*y)):
end:


#SeqLucas(y0,L): Given an alphabet Alph, a finite set of patterns Pats, and a positive integer L, outputs the first L terms of
#the number of the generating function for the generalized Lucas numbers with y=y0.
#SeqLucas(1,00);
SeqLucas:=proc(y0,L) local f,x,i,D1,f1:

f:=GenLucas(x,y0):

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

end:


#PPlucas(y0,L): All the pseudo-primes <=L in the generalized Lucas test with y=y0. Try:
#PPlucas(1,1000);
PPlucas:=proc(y0,L) local gu,i,mu,lu:

gu:=SeqLucas(y0,L):

lu:={seq(gu[ithprime(i)] mod ithprime(i),i=2*y0..2*y0+10)}:

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

lu:=lu[1]:

mu:=[]:

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


end:

###start general



#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(CGJs({0,1},{[1,1],[0,0,0]},s),s,1001,0);
PtestG:=proc(f,x,n,p0) local C:
C:=RtoC(f,x):
evalb(ValCp(C,n,n)=p0):
end:

 



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


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=4..10)}:

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:



#FindPP1G(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:
#FindPP1G(CGJsG({0,1},{[1,1],[0,0,0]},s,0),s,2,1000);
FindPP1G:=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:

#GenPerrin3 (a1, a2, a3, x): The general Perrin-style rational function with cubic denominator. 
#To get the original Perrin generating function.Try:
#GenPerrin3(0,-1,1,x);
GenPerrin3:=proc (a1, a2, a3, x): (3-2*a1*x+a2*x^2)/(1-a1*x+a2*x^2-a3*x^3): 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:

#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:


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:

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:

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:

####