######################################################################
##Condorcet.txt: Save this file as  Condorcet.txt                    #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read `Condorcet.txt` <Enter>                       #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created:  Nov. 4, 2016
 
print(`Created:  Nov. 4, 2016`):
print(` This is Condorcet.txt `):
print(`It is the Maple package that accompanyies the article `):
print(` "The Probability of a Condorcet Paradox Occuring   tends to 0.087739828045910... as the (odd) number of voters tends to infinity`):
print(`(and each voter chooses (independently) one of the six possible rankings uniformly at random)"`):

print(`by Shalosh B. Ekhad and Doron Zeilberger`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the procedures from the Maple package CompIneq.txt, type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):

Digits:=30:


ezraCompIneq:=proc()

if args=NULL then
 print(` The procedures from the Maple package CompIneq.txt are:  CE, GFtE, GFtEv,nuCE, nuCEmD, SnuCE, SnuCEm `):
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: CU, CGseq, EmpAsy, GuessQP, GuessRec, MBC, Zinn `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: Condorcet, CondorcetPodd,  CondorcetPoddPC, CondorcetQP,  CPseq, MCcon `):
 print(` `):



elif nops([args])=1 and op(1,[args])=CE then
print(`CE(VarList,SetOfIneqalities,n): inputs a list of discrete variables, a set of affine linear-expressions in the variables`):
print(`of VarList, and a non-negative integer n, outputs the SET of compositions,C, of length nops(VarList) that add up to n`):
print(`and such that if you replace  VarList[i] by C[i] in SeqOfInequalties, there are all non-negative`):
print(`For example, to get the set of integer partitions of 5 do`):
print(`CE([x,y,z], {x-y,y-z},5); `) ;
print(`For example, to get the set of all Condorcet scenariors with 3 candidates with 1->2->3->1, and 10 people voting, type`):
print(`CE([p123,p132,p213,p231,p312,p321],{p123+p132+p312-p213-p231-p321-1, p123+p213+p231-p132-p312-p321-1,p312+p321+p231-p132-p123-p213-1},10);`):

elif nops([args])=1 and op(1,[args])=CGseq then
print(`CGseq(x123,x132,x213,x231,x312,x321,N): the first N terms of the sequence of Condorcet scenariors`. Try):
print(`CGseq(x123,x132,x213,x231,x312,x321,10);`):

elif nops([args])=1 and op(1,[args])=CPseq then
print(`CPseq(N): the sequence of probabilities (times 6^n) of the Condorcet scenario 1->2->3->1 from n=1 to n=N. Try:`):
print(`CPseq(20);`):

elif nops([args])=1 and op(1,[args])=Condorcet then
print(`Condorcet(t):An article with the rational generating function for (strict) Condorect scenarios`):


elif nops([args])=1 and op(1,[args])=CondorcetPodd then
print(`CondorcetPodd(n,N,K,w):  Inputs (i)n: a discrete variable  name (ii) N: a letter denoting its shift operator`):
print(` (iii)K: a large positive integer,  (iv)w : a small positive integer `):
print(` the story of the probability of Condorect's paradox for an odd number of voters `):
print(`where each of the six rankings is equally likely. `):
print(`tells the story of the probability of Condorect's paradox for an odd number of voters`):
print(`where each of the six rankings is equally likely. Try:`):
print(`CondorcetPodd(n,N,1000,3);`):

elif nops([args])=1 and op(1,[args])=CondorcetPoddPC then
print(`CondorcetPoddPC(n,N,K,w):  faster version of CondorcetPodd(n,N,K,w), using pre-computed values. Try:`):
print(`CondorcetPoddPC(n,N,1000,3);`):

elif nops([args])=1 and op(1,[args])=CondorcetQP then
print(`CondorcetQP(n):The story of Condorcet, displaying the expressions for the number of occurence, using the variable n. Try:`):
print(`CondorcetQP(n);`):


elif nops([args])=1 and op(1,[args])=CtoR then
print(`CtoR(S,t), the rational function, in t, whose coefficients`):
print(`are the members of the C-finite sequence S. For example, try:`);
print(`CtoR([[1,1],[1,1]],t);`):

elif nops([args])=1 and op(1,[args])=CU then
print(`CU(S) given a set of expressions kicks out all non-negative  integers. Try:`):
print(` CU({x-y,5,6,z-w}); `)


elif nops([args])=1 and op(1,[args])=EmpAsy then
print(`EmpAsy(L,n0,n): Empirical asympototics of a sequence that converges to a constant in terms of 1/n.`):
print(`Inputs a list of numbers L of size k, say, and  a positive integer n0, where the L=[f(n0),..., f(n0+k-1)]`):
print(`returns an expression of the form a0+a1/n+..a[k-1]/n^(k-1), let's call it g(n) such that `):
print(`f(n)=g(n) for n=n0..n0+k01. Try`):
print(`EmpAsy([seq(4+5/i+7/i^2,i=1000..1002)],1000,n);`):

elif nops([args])=1 and op(1,[args])=GFtE then
print(`GFtE(V,S,t,N): The rational generating function whose coefficient of t^n is the number of compositions in the set`):
print(`of variables V that sum up to n, that satisfy the inequalities in V. Try:`):
print(`using N data points. Try:`):
print(`GFtE([a,b,c],{a-b,b-c},t,30);`):

elif nops([args])=1 and op(1,[args])=GFtEv then
print(`GFtEv(V,S,t,N): Verbose version of GFtE(V,S,t,N) (q.v.). Try: `):
print(`GFtEv([a,b,c],{a-b,b-c},t,30);`):

elif nops([args])=1 and op(1,[args])=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 nops([args])=1 and op(1,[args])=GuessQP then
print(`GuessQP(L,R,n,d): Inputs a list L , positive integer R,  a variable n, and a positive inteter d, and outputs a list of polynomials`):
print(` of degree <=d, of length r<=R, let's call it M, such that L[i+m*r]=M[i](i+m*r) for`):
print(` i=1..r. If nothing found, it returns FAIL. For example try: `):
print(` GuessQP([1,-1,1,-1,1,-1,1,-1,1,-1],4,n,0); `):

elif nops([args])=1 and op(1,[args])=MBC then
print(`MBC(L): The multi-binomial coefficient of the composition L`):

elif nops([args])=1 and op(1,[args])=MCcon then
print(`MCcon(n,N): an estimate for the prob. of the Condorcet scenario using a simulation with n voters with N elections`):
print(`Try: `):
print(`MCon(1001,1000); `):

elif nops([args])=1 and op(1,[args])=nuCE then
print(`nuCE(VarList,SetOfIneqalities,n): inputs a list of discrete variables, a set of affine linear-expressions in the variables`):
print(`of VarList, and a non-negative integer n, outputs the NUMBER of compositions,C, of length nops(VarList) that add up to n`):
print(`and such that if you replace  VarList[i] by C[i] in SeqOfInequalties, there are all non-negative`):
print(`For example, to get the number of integer partitions of 5 do`):
print(`nuCE([x,y,z], {x-y,y-z},5); `) ;
print(`For example, to get number of all Condorcet scenariors with 3 candidates with 1->2->3->1, and 10 people voting, type`):
print(`nuCE([p123,p132,p213,p231,p312,p321],{p123+p132+p312-p213-p231-p321-1, p123+p213+p231-p132-p312-p321-1,p312+p321+p231-p132-p123-p213-1},10);`):

elif nops([args])=1 and op(1,[args])=nuCEmD then
print(`nuCEmD(VarList,SetOfIneqalities,n): inputs a list of discrete variables, a set of affine linear-expressions in the variables`):
print(`of VarList, and a non-negative integer n, outputs the SUM  of the Multi-binomial coefficients of the`):
print(`compositions,C, of length nops(VarList) that add up to n`):
print(`and such that if you replace  VarList[i] by C[i] in SeqOfInequalties, there are all non-negative.`):
print(`It is done directly, by generating the set (using CE(V,S,N)) and then summing over MBC(c) over all c in it.`):
print(`For example, to get the set of integer partitions of 5 do`):
print(`nuCEmD([x,y,z], {x-y,y-z},5); `):
print(`For example, to get the probability of the Condorcet scenariors with 3 candidates with 1->2->3->1, and 10 people voting, type`):
print(`nuCEmD([p123,p132,p213,p231,p312,p321],{p123+p132+p312-p213-p231-p321-1, p123+p213+p231-p132-p312-p321-1,p312+p321+p231-p132-p123-p213-1},10)/6^10;`):

elif nops([args])=1 and op(1,[args])=SnuCE then
print(`SnuCE(V,S,N):  The first N terms of the sequence nuCE(V,S,n) (q.v.).Try:`):
print(`SnuCE([p123,p132,p213,p231,p312,p321],{p123+p132+p312-p213-p231-p321-1, p123+p213+p231-p132-p312-p321-1,p312+p321+p231-p132-p123-p213-1},14);`):

elif nops([args])=1 and op(1,[args])=SnuCEm then
print(`SnuCEm(V,S,N):  The first N terms of the sequence nuCEm(V,S,n) (q.v.).Try:`):
print(`SnuCEm([p123,p132,p213,p231,p312,p321],{p123+p132+p312-p213-p231-p321-1, p123+p213+p231-p132-p312-p321-1,p312+p321+p231-p132-p123-p213-1},14);`):

elif nops([args])=1 and op(1,[args])=Zinn then
print(`Zinn(resh): Zinn-Justin's method to estimate`):
print(`the C,mu, and theta such that`):
print(`resh[i] is appx. Const*mu^i*i^theta`):
print(`For example, try:`):
print(`Zinn([seq(5*i*2^i,i=1..30)]);`):

else
print(`There is no ezra for`,args):
fi:
 
end:


#From Findrec
ezraFindrec:=proc()
if args=NULL then

print(` FindRec: A Maple package for empirically guessing partial recurrence`):
print(`equations satisfied by Discrete Functions of ONE Variable`):
print():
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(`Contains procedures:  `):
print(`  findrec, Findrec, FindrecF`):
print():

elif nargs=1 and args[1]=findrec then
print(`findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating`):
print(`the sequence f of degree DEGREE and order ORDER.`):
print(`For example, try: findrec([seq(i,i=1..10)],0,2,n,N);`):

elif nargs=1 and args[1]=Findrec then
print(`Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with`):
print(`poly coffs. ope(n,N), where n is the discrete variable, and N is the shift operator `):
print(`of maximum DEGREE+ORDER<=MaxC`):
print(`e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);`):

elif nargs=1 and args[1]=FindrecF then
print(`FindrecF(f,n,N): Given a function f of a single variable tries to find a linear recurrence equation with`):
print(`poly coffs. .g. try FindrecF(i->i!,n,N);`):


elif nargs=1 and args[1]=SeqFromRecF then
print(`SeqFromRecF(ope,n,N,Ini,K): Given the first L-1`):
print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`):
print(`satisfied by the recurrence ope(n,N)f(n)=0`):
print(`extends it to the first K values`):
print(`For example, try:`):
print(`SeqFromRecF(N-n-1,n,N,[1],10);`):

elif nargs=1 and args[1]=SeqFromRecFw then
print(`SeqFromRecFw(ope,n,N,Ini,K,w): Given the first L-1`):
print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`):
print(`satisfied by the recurrence ope(n,N)f(n)=0`):
print(`extends outputs the values from K-w+1 to K`):
print(`For example, try:`):
print(`SeqFromRecFw(N-n-1,n,N,[1],10,2);`):


fi:
 
end:


###Findrec 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrecVerbose:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
print(`There is some slack, there are `, nops(ope)):
print(ope):
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
option remember:


if not findrecEx(f,DEGREE,ORDER,ithprime(20)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(80)) then
 RETURN(FAIL):
fi:


if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 


Yafe:=proc(ope,N) local i,ope1,coe1,L: 
if ope=0 then
 RETURN(1,0):
fi:
ope1:=expand(ope):
L:=degree(ope1,N):
coe1:=coeff(ope1,N,L):
ope1:=normal(ope1/coe1):
ope1:=normal(ope1):
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),ope1:
end:


#Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with
#poly coffs.
#of maximum DEGREE+ORDER<=MaxC
#e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);
Findrec:=proc(f,n,N,MaxC)
local DEGREE, ORDER,ope,L:

for L from 0 to MaxC do
for ORDER from 0 to L do
 DEGREE:=L-ORDER:
if (2+DEGREE)*(1+ORDER)+4>=nops(f) then
 print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
 RETURN(FAIL):
fi:
 ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):
     if ope<>FAIL then
       RETURN(ope):
      fi:
 od:
od:
FAIL:

end:





 
#SeqFromRecF(ope,n,N,Ini,K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#extends it to the first K values
SeqFromRecF:=proc(ope,n,N,Ini,K)
local ope1,gu,L,n1,j1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

ope1:=expand(subs(n=n-L,ope1)/N^L):

gu:=Ini:

for n1 from nops(Ini)+1 to K do
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)]:
od:

gu:

end:
 
 
#SeqFromRecFw(ope,n,N,Ini,K,w): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#outputs the values from K to K+w-1
SeqFromRecFw:=proc(ope,n,N,Ini,K,w)
local ope1,gu,L,n1,j1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

ope1:=expand(subs(n=n-L,ope1)/N^L):

gu:=Ini:

for n1 from nops(Ini)+1 to w do
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)]:
od:

for n1 from w+1 to K+w-1 do
 gu:=[op(2..nops(gu),gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),j1=1..L)]:
od:

gu:

end:
 
#end Findrec


with(linalg):

#findrecEx(f,DEGREE,ORDER,m1): Explores whether thre
#is a good chance that there is a recurrence of degree DEGREE
#and order ORDER, using the prime m1
#For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003);
findrecEx:=proc(f,DEGREE,ORDER,m1)
local ope,var,eq,i,j,n0,eq1,a,A1,
D1,E1,Eq,Var,f1,n,N:
option remember:
f1:=f mod m1:
if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f1) mod m1:
  od:
 
   eq:= eq union {eq1}:
od:


Eq:= convert(eq,list):
Var:= convert(var,list):


D1:=nops(Var):
E1:=nops(Eq):
if E1<D1 then
  RETURN(true):
fi:

A1:=matrix(D1,D1):

for i from 1 to D1-1 do
for j from 1 to D1 do
  A1[i,j]:=coeff(Eq[i],Var[j]):
od:
od:

for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1],Var[j]):
od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:

if E1-D1>=1 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+1],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

if E1-D1>=2 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+2],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

true:

end:
#End From Findrec

##########From C-finite

sn:=proc(resh,n1):
-1/log(op(n1+1,resh)*op(n1-1,resh)/op(n1,resh)^2):
end:
 
#Zinn(resh): Zinn-Justin's method to estimate
#the C,mu, and theta such that
#resh[i] is appx. Const*mu^i*i^theta
#For example, try:
#Zinn([seq(5*i*2^i,i=1..30)]);
Zinn:=proc(resh)
local s1,s2,theta,mu,n1,i:
if nops({seq(sign(resh[i]),i=1..nops(resh))})<>1 then
 RETURN(FAIL):
fi:

n1:=nops(resh)-1:
s1:=sn(resh,n1):
s2:=sn(resh,n1-1):
theta:=evalf(2*(s1+s2)/(s1-s2)^2):
mu:=evalf(sqrt(op(n1+1,resh)/op(n1-1,resh))*exp(-(s1+s2)/((s1-s2)*s1))):
[theta,mu]:
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:

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

#CtoR(S,t), the rational function, in t, whose coefficients
#are the members of the C-finite sequence S. For example, try:
#CtoR([[1,1],[1,1]],t);

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(factor(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:


###ene from Cfinite

##from GuessQuasiPol

#GuessQPold(L,R,n): Inputs a list L and outputs a list of polynomials
#of length r<=R, let's call it M, such that L[i+m*r]=M[i](i+m*r) for
#i=1..r. If nothing found, it returns FAIL. For example try:
#GuessQP([1,-1,1,-1,1,-1,1,-1,1,-1],4,n);


GuessQPold:=proc(L,R,n) local r,gu:
 for r from 1 to R do
  gu:=GuessQP1(L,r,n):
  if gu<>FAIL then
    RETURN(gu):
  fi:
 od:
FAIL:
end:


#GuessQP1(L,r,n): Inputs a list L and outputs a list of polynomials
#of length r, let's call it M, such that L[i+m*r]=M[i](i+m*r) for
#i=1..r, For example try:
#GuessQP1:=proc([1,-1,1,-1,1,-1,1,-1,1,-1],2,n);
GuessQP1:=proc(L,r,n) local M,i,L1,gu,m:
M:=[]:
for i from 1 to r do
L1:=[seq(L[i+(m-1)*r],m=1..trunc((nops(L)-i)/r)+1)]:
gu:=GuessPol(L1,m):

if gu=FAIL then
  RETURN(FAIL):
fi:
gu:=expand(subs(m=(n-i)/r+1,gu)):
M:=[op(M),gu]:
end:
M:
end:

#GuessPol1(L,d,n): guesses a polynomial of degree d in n for
# the list L, such that P(i)=L[i+1] for i=1..nops(L)
#For example, try: 
#GuessPol1([seq(i,i=1..10)],1,n);
GuessPol1:=proc(L,d,n) local P,i,a,eq,var:
if d>nops(L)-3 then
 ERROR(`the list is too small`):
fi:

P:=add(a[i]*n^i,i=0..d):
var:={seq(a[i],i=0..d)}:
eq:={seq(subs(n=i,P)-L[i],i=1..d+1)}:


var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
fi:

P:=subs(var,P):
if {seq(subs(n=i,P)-L[i],i=d+1..nops(L))}<>{0} then
 FAIL:
 else
P:
fi:


end:

#GuessPol(L,n): guesses a polynomial of degree d in n for
# the list L, such that P(i)=L[i] for i=1..nops(L) for example, try: 
#GuessPol([seq(i,i=1..10)],n);
GuessPol:=proc(L,n) local d,gu:

for d from 0 to nops(L)-3 do

 gu:=GuessPol1(L,d,n):
 if gu<>FAIL then
    RETURN(gu):
 fi:
od:

FAIL:

end:
##End from GuessQuasiPol

#CU(S) given a set of expressions kicks out all non-negative  integers
CU:=proc(S) local s,S1:
S1:=S:

for s in S do
if type(s,integer) and s>=0 then
  S1:=S1 minus {s}
fi:
od:
S1:
end:


#CE(VarList,SetOfIneqalities,n): inputs a list of discrete variables, a set of affine linear-expressions in the variables
#of VarList, and a non-negative integer n, outputs the SET of compositions,C, of length nops(VarList) that add up to n
#and such that if you replace  VarList[i] by C[i] in SeqOfInequalties, there are all non-negative
#For example, to get the set of integer partitions of 5 do
#CE([x,y,z], {x-y,y-z},5)
#For example, to get the set of all Condorcet scenariors with 3 candidates with 1->2->3->1, and 10 people voting, type
#CE([p123,p132,p213,p231,p312,p321],{p123+p132+p312-p213-p231-p321-1, p123+p213+p231-p132-p312-p321-1,p312+p321+p231-p132-p123-p213-1},10);
CE:=proc(V,S,n) local K,V1,s,gu,i,gu1,gu11,S1:
option remember:

for s in S do
 if type(s,integer) and s<0 then
   RETURN({}):
 fi:
od:


if V=[] then
 if n=0 then
  RETURN({[]}):
  else
   RETURN({}):
 fi:
fi:


K:=nops(V):
V1:=[op(1..K-1,V)]:
gu:={}:

for i from 0 to n do
 S1:=CU(subs(V[K]=i,S)):
 gu1:=CE(V1,S1,n-i):
 gu:=gu union {seq([op(gu11),i],gu11 in gu1)}:
od:
gu:
end:


#nuCE(VarList,SetOfIneqalities,n): inputs a list of discrete variables, a set of affine linear-expressions in the variables
#of VarList, and a non-negative integer n, outputs the NUMBER of compositions,C, of length nops(VarList) that add up to n
#and such that if you replace  VarList[i] by C[i] in SeqOfInequalties, there are all non-negative
#For example, to get the NUMBER of integer partitions of 5 do
#nuCE([x,y,z], {x-y,y-z},5)
#For example, to get the NUMBER of all Condorcet scenariors with 3 candidates with 1->2->3->1, and 10 people voting, type
#nuCE([p123,p132,p213,p231,p312,p321],{p123+p132+p312-p213-p231-p321-1, p123+p213+p231-p132-p312-p321-1,p312+p321+p231-p132-p123-p213-1},10);
nuCE:=proc(V,S,n) local K,V1,s,gu,i,S1:
option remember:

for s in S do
 if type(s,integer) and s<0 then
   RETURN(0):
 fi:
od:


if V=[] then
 if n=0 then
  RETURN(1):
  else
   RETURN(0):
 fi:
fi:


K:=nops(V):
V1:=[op(1..K-1,V)]:
gu:=0:

for i from 0 to n do
 S1:=CU(subs(V[K]=i,S)):
 gu:=gu+nuCE(V1,S1,n-i):
od:
gu:
end:



#SnuCE(V,S,N):  The first N terms of the sequence nuCE(V,S,n) (q.v.).
SnuCE:=proc(V,S,N)  local n:
[seq(nuCE(V,S,n),n=1..N)]:
end:

#SnuCE0(V,S,N):  The first N terms of the sequence nuCE(V,S,n) (q.v.).
SnuCE0:=proc(V,S,N)  local n:
[seq(nuCE(V,S,n),n=0..N)]:
end:


#CondorcetQP(n):The story of Condorcet, displaying the expressions for the number of occurence, using the variable n. Try:
#CondorcetQP(n);
CondorcetQP:=proc(n) local gu,mu,t0,lu1, p123,p132,p213,p231,p312,p321,i:
t0:=time():
print(`How likely is Condorcet's paradox?`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Suppose that there three candidates and n voters, the sequence for the number of ways (out of binomial(n+5,5)) possible ways`):
print(`in which the votes can be cast such that there is a cycle 1>2>3>1 is`):
gu:=SnuCE([p123,p132,p213,p231,p312,p321],{p123+p132+p312-p213-p231-p321-1, p123+p213+p231-p132-p312-p321-1,p312+p321+p231-p132-p123-p213-1},16):
print(gu):
print(`hence the number of ways in which there is a cycle 1>2>3>1 or 1>3>2>1, is twice that`):
gu:=2*gu:
print(``):
print(` namely `):
print(``):
print(gu):
print(``):	
mu:=factor(GuessQP(gu,2,n,5));
print(``):
print(`-------------------------------------------------------------------------------------`):
print(``):

print(`Case I: n is odd`):
print(``):
print(`The number of ways of getting a Condorcet scenario when n is odd is`):
print(``):
print(mu[1]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(mu[1]):
print(``):
print(`dividing by the total, binomial(n+5,5), the probabilty is`):
lu1:=factor(mu[1]/expand(binomial(n+5,5))):
print(lu1):
print(``):	
print(`Hence the probabilty that there is NO Condorcet scenario, i.e. there is a clear-cut ranking is`):
print(``):
lu1:=factor(normal(1-lu1)):
print(``):
print(lu1):
print(``):
print(`and in Maple notation`):
print(``):
lprint(lu1):
print(``):
print(`This agrees with William Gehrlein and Peter Fishburn's result in "Condorcet's paradox and anonymous preferance profiles"`):
print(`published in Public Choice v. 26 (1976), 1-18`):
print(`and also quoted in William Gehrlein's article "The Expected Probability of the Condorceet Paradox"`):
print(`published in Economic Letters v. 7 (1981). 33-37 `):

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

print(`Case II: n is even`):
print(``):
print(`The number of ways of getting a Condorcet scenario when n is odd is`):
print(``):
print(mu[2]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(mu[2]):
print(``):
print(`dividing by the total, binomial(n+5,5), the probabilty is`):
lu1:=factor(mu[2]/expand(binomial(n+5,5))):
print(lu1):
print(``):	
print(`Hence the probabilty that there is NO Condorcet scenario, i.e. there is a clear-cut ranking is`):
print(``):
lu1:=factor(normal(1-lu1)):
print(``):
print(lu1):
print(``):
print(`and in Maple notation`):
print(``):
print(lu1):
print(``):
print(`To conclude, for the sake of the OEIS, here are the first 100 terms of the number of ways Condorcet can occur, using the above quasi-polynomial`):
print(``):
print([seq(op([subs(n=2*i-1,mu[1]),  subs(n=2*i,mu[2])]),i=1..50)]):
print(``):
print(`------------------------------------------------------------------`):
print(``):
print(`This concludes this article that took`, time()-t0, `seconds to generate.`):
end:



#GuessQP1d(L,r,n,d): Inputs a list L and outputs a list of polynomials of degree d 
#of length <=r, let's call it M, such that L[i+m*r]=M[i](i+m*r) for
#i=1..r, For example try:
#GuessQP1d:=proc([1,-1,1,-1,1,-1,1,-1,1,-1],2,n,0);
GuessQP1d:=proc(L,r,n,d) local gu,r1:

for r1 from 1 to r do

 gu:=GuessQP1rd(L,r1,n,d):

 if gu<>FAIL then

       RETURN(gu):
 fi:
od:
FAIL:
end:

#GuessQP1rd(L,r,n,d): Inputs a list L and outputs a list of polynomials of degree d 
#of length r, let's call it M, such that L[i+m*r]=M[i](i+m*r) for
#i=1..r, For example try:
#GuessQP1rd:=proc([1,-1,1,-1,1,-1,1,-1,1,-1],2,n,0);
GuessQP1rd:=proc(L,r,n,d) local M,i,L1,gu,m:
M:=[]:
for i from 1 to r do
L1:=[seq(L[i+(m-1)*r],m=1..trunc((nops(L)-i)/r)+1)]:

gu:=GuessPol1(L1,d,m):

if gu=FAIL then
  RETURN(FAIL):
fi:
gu:=expand(subs(m=(n-i)/r+1,gu)):
M:=[op(M),gu]:
end:
M:
end:

#GuessQPd(L,R,n,d): Inputs a list L and outputs a list of polynomials of degree d
#of length r<=R, let's call it M, such that L[i+m*r]=M[i](i+m*r) for
#i=1..r. If nothing found, it returns FAIL. For example try:
#GuessQPd([1,-1,1,-1,1,-1,1,-1,1,-1],4,n,0);


GuessQPd:=proc(L,R,n,d) local r,gu:
 for r from 1 to R do
  gu:=GuessQP1d(L,r,n,d):
  if gu<>FAIL then
    RETURN(gu):
  fi:
 od:
FAIL:
end:


#GuessQP(L,R,n,d): Inputs a list L and outputs a list of polynomials
#of length r<=R, let's call it M, such that L[i+m*r]=M[i](i+m*r) for
#i=1..r. If nothing found, it returns FAIL. For example try:
#GuessQP([1,-1,1,-1,1,-1,1,-1,1,-1],4,n);
GuessQP:=proc(L,R,n,d) local gu,d1:

 for d1 from 0 to d do
  gu:=GuessQPd(L,R,n,d1):
  if gu<>FAIL then
    RETURN(gu):
  fi:
 od:
FAIL:
end:




#Condorcet(t):An article with the rational generating function for (strict) Condorect scenarios
Condorcet:=proc(t) local gu,mu,t0, p123,p132,p213,p231,p312,p321,n:
t0:=time():
print(`How Many  Condorcet Scenarios Are There?`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Suppose that there three candidates and n voters, the sequence for the number of ways (out of binomial(n+5,5)) possible ways`):
print(`in which the votes can be cast such that there is a cycle 1>2>3>1 is`):
gu:=SnuCE([p123,p132,p213,p231,p312,p321],{p123+p132+p312-p213-p231-p321-1, p123+p213+p231-p132-p312-p321-1,p312+p321+p231-p132-p123-p213-1},26):
print(gu):
print(`hence the number of ways in which there is a cycle 1>2>3>1 or 1>3>2>1, is twice that`):
gu:=2*gu:
gu:=[0,op(gu)]:
mu:=CtoR(GuessRec(gu),t):
print(``):
print(` Let a(n) be the number of Condorcet scenarios, then`):
print(``):
print(Sum(a(n)*t^n,n=0..infinity)=mu):
print(``):	
print(`------------------------------------------------`):
print(``):

print(`This concludes this article that took`, time()-t0, `seconds to generate.`):
end:

#GFtE(V,S,t,N): The rational generating function whose coefficient of t^n is the number of compositions in the set
#of variables V that sum up to n, that satisfy the inequalities in V. Try:
#using N data points. Try:
#GFtE([a,b,c],{a-b,b-c},t,30);
GFtE:=proc(V,S,t,N) local gu:
gu:=GuessRec(SnuCE0(V,S,N)):  

if gu=FAIL then
  RETURN(FAIL):
else
 RETURN(CtoR(gu,t)):
fi:

end:


#MBC(L): The multi-binomial coefficient of the composition L
MBC:=proc(L) local i: convert(L,`+`)!/mul(L[i]!,i=1..nops(L)):end:


#nuCEmD(VarList,SetOfIneqalities,n): inputs a list of discrete variables, a set of affine linear-expressions in the variables
#of VarList, and a non-negative integer n, outputs the SUM  of the Multi-binomial coefficients of the
#compositions,C, of length nops(VarList) that add up to n
#and such that if you replace  VarList[i] by C[i] in SeqOfInequalties, there are all non-negative.
#It is done directly, by generating the set (using CE(V,S,N)) and then summing over MBC(c) over all c in it.
#It is done DIRECTLY
#For example, to get the set of integer partitions of 5 do
#nuCEmD([x,y,z], {x-y,y-z},5);
#For example, to get the probability of the Condorcet scenariors with 3 candidates with 1->2->3->1, and 10 people voting, type
#nuCEmD([p123,p132,p213,p231,p312,p321],{p123+p132+p312-p213-p231-p321-1, p123+p213+p231-p132-p312-p321-1,p312+p321+p231-p132-p123-p213-1},10)/6^10;
nuCEmD:=proc(V,S,n) local gu,g:

gu:=CE(V,S,n):

add(MBC(g), g in gu):

end:





#nuCEm(VarList,SetOfIneqalities,n): inputs a list of discrete variables, a set of affine linear-expressions in the variables
#of VarList, and a non-negative integer n, outputs the sum of the Multi-Binomial Coefficient of compositions,C, of length nops(VarList) that add up to n
#and such that if you replace  VarList[i] by C[i] in SeqOfInequalties, there are all non-negative
#For example, to get the set of integer partitions of 5 do
#nuCEm([x,y,z], {x-y,y-z},5)
#For example, to get the set of all Condorcet scenariors with 3 candidates with 1->2->3->1, and 10 people voting, type
#nuCEm([p123,p132,p213,p231,p312,p321],{p123+p132+p312-p213-p231-p321-1, p123+p213+p231-p132-p312-p321-1,p312+p321+p231-p132-p123-p213-1},10);
nuCEm:=proc(V,S,n) local K,V1,s,gu,i,S1:
option remember:

for s in S do
 if type(s,integer) and s<0 then
   RETURN(0):
 fi:
od:


if V=[] then
 if n=0 then
  RETURN(1):
  else
   RETURN(0):
 fi:
fi:


K:=nops(V):
V1:=[op(1..K-1,V)]:
gu:=0:

for i from 0 to n do
 S1:=CU(subs(V[K]=i,S)):
 gu:=gu+binomial(n,i)*nuCEm(V1,S1,n-i):
od:
gu:
end:

#SnuCEm(V,S,N):  The first N terms of the sequence numCE(V,S,n) (q.v.).
SnuCEm:=proc(V,S,N)  local n:
[seq(nuCEm(V,S,n),n=1..N)]:
end:


#CPseq(N): the sequence of probabilities (times 6^n) of the Condorcet scenario 1->2->3->1 from n=1 to n=N
CPseq:=proc(N) local p123,p132,p213,p231,p312,p321:
SnuCEm([p123,p132,p213,p231,p312,p321],{p123+p132+p312-p213-p231-p321-1, p123+p213+p231-p132-p312-p321-1,p312+p321+p231-p132-p123-p213-1},N):
end:



#CGseqSlow(x123,x132,x213,x231,x312,x321,N): the first N terms of the sequence of Condorcet scenariors
CGseq:=proc(x123,x132,x213,x231,x312,x321,N) local n,mu,t:
mu:=(t^3*x123*x231*x312+1)*t^3*x312*x123*x231/(t^2*x123*x321-1)/(t^2*x213*x312-1)/(t^2*x312*x123-1)/(t^2*x132*x231-1)/(t^2*x123*x231-1)/(t^2*x231*x312-1):
mu:=taylor(mu,t=0,N+1):
[seq(coeff(mu,t,n),n=1..N)]:
end:


#CT1(F,var): the constant term of the Laurent polynomial F in the list of variables var
CT1:=proc(F,var) local i,F1:
F1:=F:

  for i from 1 to nops(var) do
   F1:=coeff(F1,var[i],0):
  od:
F1:
end:


#CPseq(N): the sequence of probabilities (times 6^n) of the Condorcet scenario 1->2->3->1 from n=1 to n=N done Fast
CPseq:=proc(N) local x123,x132,x213,x231,x312,x321,i,gu:
gu:=expand(CGseq(x123,x132,x213,x231,x312,x321,N)):

[seq(CT1(gu[i]*(1/x123+1/x132+1/x213+1/x231+1/x312+ 1/x321)^i,[x123,x132,x213,x231,x312,x321]),i=1..N)]:
end:



#GFtEv(V,S,t,N): Verbose version of GFtE(V,S,t,N) (q.v.)
#using N data points. Try:
#GFtEv([a,b,c],{a-b,b-c},t,30);
GFtEv:=proc(V,S,t,N) local gu,t0,i,n:
t0:=time():
gu:=GFtE(V,S,t,N):

if gu=FAIL then
print(`You need to make`, N, `bigger `):
  RETURN(FAIL):
fi:

print(`The Ordinary Generating Function for the Number of Compositions Satisfying a Certain set of inequalities`):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):
print(`Theorem: Let a(n) be the number of solutions, in non-negative integers, of `):
print(``):
print(add(V[i],i=1..nops(V))=n):
print(``):
print(`subject to the inequalities`):
print(``):
for i from 1 to nops(S) do
 print(S[i]>=0):
od:
print(``):
print(`Then the ordinary generating function of a(n) is given by the rational function`):
print(``):
print(Sum(a(n)*t^n,n=0..infinity)=gu):
print(``):
print(`and in Maple format it is`):
print(``):
lprint(gu):
print(``):
print(`--------------------------------`):
print(``):
print(`This ends this article that took`, time()-t0, `seconds to generate `):
print(``):

end:




#EmpAsy(L,n0,n): Empirical asympototics of a sequence that converges to a constant in terms of 1/n.
#Inputs a list of numbers L of size k, say, and  a positive integer n0, where the L=[f(n0),..., f(n0+k-1)]
#returns an expression of the form a0+a1/n+..a[k-1]/n^(k-1), let's call it g(n) such that 
#f(n)=g(n) for n=n0..n0+k01. Try
#EmpAsy([seq(4+5/i+7/i^2,i=1000..1002)],1000,n);
EmpAsy:=proc(L,n0,n) local k,var,eq,a,i,g:
k:=nops(L):
var:={seq(a[i],i=0..k-1)}:
g:=add(a[i]/n^i,i=0..k-1):
eq:={seq(subs(n=n0+i-1,g)-L[i],i=1..k)}:
var:=solve(eq,var):
subs(var,g):
end:


#CondorcetPodd(n,N,K,w):  Inputs (i)n: a discrete variable  name (ii) N: a letter denoting its shift operator
#(iii)K: a large positive integer,  (iv)w : a small positive integer
#the story of the probability of Condorect's paradox for an odd number of voters
#where each of the six rankings is equally likely. Try:
CondorcetPodd:=proc(n,N,K,w) local t0,gu,ope,lu,i,ku1,ku2,mu1,mu2,sfarot,i1,ku2a:
t0:=time():
print(`The probability that a Condorcet Scenario will occur with an odd number of voters and 3 Candidates`):
print(`Where each voter chooses one of the six rankings with the same probability (i.e. 1/6)`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

gu:=CPseq(60):

gu:=[seq(gu[2*i-1],i=1..30)]:

ope:=Findrec(gu,n,N,10):

if gu=FAIL then
  RETURN(FAIL):
fi:

print(`Theorem: Suppose three candidates, 1,2,3, are running for office, and there 2n-1 voters, each of them choosing one of the 3!=6`):
print(`possible ranking with the same probability (1/6). The votes are counted and it turns out that in the vote of `):
print(` 1 vs. 2, 1 won`):
print(` 2 vs. 3, 2 won`):
print(` 1 vs. 3, 3 won`):
print(``):
print(`In other words, there is a cycle 1->2->3->1 `):
print(`Let a(n) be the probability of that happening times 6^(2n-1) (i.e. the numerator)`):
print(`The sequence a(n) satisfies the following recurrence `):
print(``):
print(add(factor(coeff(ope,N,i))*a(n+i),i=0..degree(ope,N))=0):
print(``):
print(`and in Maple format`):
print(``):
lprint(add(factor(coeff(ope,N,i))*a(n+i),i=0..degree(ope,N))=0):
print(``):
print(`Using this recurrence it is easy to compute many values, for example, the probability for 1999 voters is`):

lu:=evalf(SeqFromRecFw(ope,n,N,[seq(gu[i],i=1..3)],1000,3)[1]/6^(2*1000-1),20):
print(evalf(lu,30)):


mu1:=evalf(SeqFromRecFw(ope,n,N,[seq(gu[i],i=1..3)],K,w)):
mu1:=[seq(mu1[i]/6^(2*(K+i-1)-1),i=1..nops(mu1))]:
ku1:=EmpAsy(mu1,K,n):

mu2:=evalf(SeqFromRecFw(ope,n,N,[seq(gu[i],i=1..3)],2*K,w)):
mu2:=[seq(mu2[i]/6^(2*(2*K+i-1)-1),i=1..nops(mu2))]:
ku2:=EmpAsy(mu2,2*K,n):

ku1:=ku2-ku1:
print(`ku1 is`, ku1):

sfarot:=[seq(trunc(-log[10](abs(coeff(ku1,n,-i1))))-2	 ,i1=0..-ldegree(ku1,n))]:
print(`sfarot is`, sfarot):
	      


ku2a:=0:

for i1 from 0 to -ldegree(ku2,n) do
 if sfarot[i1+1]>=5 then
  ku2a:=ku2a+evalf(coeff(ku2,n,-i1),sfarot[i1+1])/n^i1:
 fi:
od:

ku2:=ku2a:
print(`The asymptotic expression is estimated to be`):

print(ku2):
print(``):
print(`and in Maple notation`):
print(``):
lprint(ku2):
print(``):
print(`But it is also possible to have 1->3->2->1, so the probability of a Condorcet scenario is`):
print(``):
ku2:=2*ku2:
print(ku2):
print(``):
print(`----------------------------------------`):
print(``):
print(`This ends this article, that took`, time()-t0, `seconds to generate `):
RETURN([ope, ku2]):

end:

#CondorcetPoddPC(n,N,K,w):  Inputs (i)n: a discrete variable  name (ii) N: a letter denoting its shift operator
#(iii)K: a large positive integer,  (iv)w : a small positive integer
#the story of the probability of Condorect's paradox for an odd number of voters
#where each of the six rankings is equally likely. Try:
CondorcetPoddPC:=proc(n,N,K,w) local t0,gu,ope,lu,i,ku1,ku2,mu1,mu2,sfarot,i1,ku2a:
t0:=time():
print(`The probability that a Condorcet Scenario will occur with an odd number of voters and 3 Candidates`):
print(`Where each voter chooses one of the six rankings with the same probability (i.e. 1/6)`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):


gu:=[0, 0, 6, 0, 270, 90, 10500, 6720, 392910, 355950, 14478156, 16439346, 
529404876, 706810104, 19272783816, 29111692896, 699677161326, 1165747995126, 
25353979229436, 45774091141350, 917549732994180, 1771818950042856, 
33174260647148088, 67845990480870498, 1198569964080659820, 2576217011384276520,
43279979034548697720, 97172396306861238720, 1562151084311426575320, 
3645494682112649782440, 56364895696677322208400, 136156526172748961935200, 
2033162581408415428229550, 5066507265019626353202630, 
73322069395323981081793500, 187939363676310774342456750, 
2643706560575331509749420500, 6952888937871643394146357800, 
95306151219566977999148535000, 256632941320578733939956369270, 
3435329091855352094619142017060, 9453458631038672249283002481720, 
123812379936714903942117773215080, 347623412550354359085613206886560, 
4461851663701935271912084788763560, 12763153196675804005956747381942360, 
160778258838191766507231740405560560, 467965621751281288493359486772101410, 
5793026200170719131181221470879740460, 17137271130912909801809606767015420200,
208715023775792055004589499371039269656, 
626897384888398383786501373897438885536, 
7519265559117988709964896895641635381416, 
22910043002293523401246081641910619333304, 
270877893397621394036184987041567724515376, 
836510093539168123146589556095498183292256, 
9757771639085931004530711120359767623040536, 
30518848286150870253293452807846195757365656, 
351486666796749948852186133697196078418141296, 
1112623225434768542071755478271892713782039176]:
#gu:=CPseq(60):

gu:=[seq(gu[2*i-1],i=1..30)]:

ope:=Findrec(gu,n,N,10):

if gu=FAIL then
  RETURN(FAIL):
fi:

print(`Theorem: Suppose three candidates, 1,2,3, are running for office, and there 2n-1 voters, each of them choosing one of the 3!=6`):
print(`possible ranking with the same probability (1/6). The votes are counted and it turns out that in the vote of `):
print(` 1 vs. 2, 1 won`):
print(` 2 vs. 3, 2 won`):
print(` 1 vs. 3, 3 won`):
print(``):
print(`In other words, there is a cycle 1->2->3->1 `):
print(`Let a(n) be the probability of that happening times 6^(2n-1) (i.e. the numerator)`):
print(`The sequence a(n) satisfies the following recurrence `):
print(``):
print(add(factor(coeff(ope,N,i))*a(n+i),i=0..degree(ope,N))=0):
print(``):
print(`and in Maple format`):
print(``):
lprint(add(factor(coeff(ope,N,i))*a(n+i),i=0..degree(ope,N))=0):
print(``):
print(`Using this recurrence it is easy to compute many values, for example, the probability for 1999 voters is`):

lu:=evalf(SeqFromRecFw(ope,n,N,[seq(gu[i],i=1..3)],1000,3)[1]/6^(2*1000-1),20):
print(evalf(lu,30)):


mu1:=evalf(SeqFromRecFw(ope,n,N,[seq(gu[i],i=1..3)],K,w)):
mu1:=[seq(mu1[i]/6^(2*(K+i-1)-1),i=1..nops(mu1))]:
ku1:=EmpAsy(mu1,K,n):

mu2:=evalf(SeqFromRecFw(ope,n,N,[seq(gu[i],i=1..3)],2*K,w)):
mu2:=[seq(mu2[i]/6^(2*(2*K+i-1)-1),i=1..nops(mu2))]:
ku2:=EmpAsy(mu2,2*K,n):

ku1:=ku2-ku1:
print(`ku1 is`, ku1):

sfarot:=[seq(trunc(-log[10](abs(coeff(ku1,n,-i1))))-2	 ,i1=0..-ldegree(ku1,n))]:



ku2a:=0:

for i1 from 0 to -ldegree(ku2,n) do
 if sfarot[i1+1]>=5 then
  ku2a:=ku2a+evalf(coeff(ku2,n,-i1),sfarot[i1+1])/n^i1:
 fi:
od:

ku2:=ku2a:
print(`The asymptotic expression is estimated to be`):

print(ku2):
print(``):
print(`and in Maple notation`):
print(``):
lprint(ku2):
print(``):
print(`But it is also possible to have 1->3->2->1, so the probability of a Condorcet scenario is`):
ku2:=2*ku2:
print(``):
print(ku2):
print(``):

print(`And in Maple format`):
lprint(ku2):
print(``):

print(`----------------------------------------`):
print(``):	
print(`This ends this article, that took`, time()-t0, `seconds to generate `):
RETURN([ope, ku2]):

end:
 



###start simulation

    #Inputs a composition,L of n into non-neg. integers of size 6
    #[123,132,213,231,312,321], decides whether 1->2->3->1
    Con1:=proc(L)
     evalb(
       #1->2
        (L[1]+L[2]+L[5])-(L[3]+L[4]+L[6])>0
        and
        #2->3
         (L[1]+L[3]+L[4])- (L[2]+L[5]+L[6])>0
        and
        #3->1
         (L[4]+L[5]+L[6])- (L[1]+L[2]+L[3])>0
       ):
end:

 #Inputs a composition,L of n into non-neg. integers of size 6
    #[123,132,213,231,312,321], decides whether 1->3->2->1
    Con2:=proc(L)
     evalb(
       #1->2
        (L[1]+L[2]+L[5])-(L[3]+L[4]+L[6])<0
        and
        #2->3
         (L[1]+L[3]+L[4])- (L[2]+L[5]+L[6])<0
        and
        #3->1
         (L[4]+L[5]+L[6])- (L[1]+L[2]+L[3])<0
       ):
end:

  #Con(L): Is L a Condorect paradox example?  

 Con:=proc(L):    Con1(L) or Con2(L)     : end:


  #RV(): generates one of [1,0$5], ..., [0$5,1] each with prob. 1/6
  RV:=proc() local i:
  i:=rand(1..6)():
  [0$(i-1),1,0$(6-i)]:
  end:

  Rcom:=proc(n) local i: add(RV(),i=1..n): end:
  
  #an estimate for the prob. of the Condorcet scenario using a simulation with n voters with N elections
  MCcon:=proc(n,N) local L,c,i:
  c:=0:
   for i from 1 to N do
    L:=Rcom(n):
    if Con(L) then
      c:=c+1:
    fi:
   od:
    
   evalf(c/N):

  end:

###end simulation