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


print(`First Written: Sept. 2020: tested for Maple 2018 `):
print(`Version  : Sept. 2020  `):
print():
print(`This is  JumpModel.txt, A Maple package`):
print(`accompanying the two papers  `):
print(` Katriel H, Mahanaimi U, Koutschan C, Zeilberger D, Steel M, and Snir S, Gene Order modeling via generating functions.`):
print(` Katriel H, Mahanaimi U, Szklarczyk D, Koutschan C, Zeilberger D, Steel M, and Snir S, Gene Order Phylogenetics is Consistent under a Markov Process`):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/XXX .`):
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(`For a list of the supporting functions type: ezra1();`):
print():
 



ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are`):
print(``):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(`The MAIN procedures are: `):
print(` pijt, pijtF, qkt, qktF, , qktFF, qktPE, qktPF `):

elif nargs=1 and args[1]=pijt then
print(`pijt(i,j): The rational function p_{ij}(t) given by Eq. (15) in the paper, using Eq. (15) `):
print(`For example, try: `):
print(`pijt(10,20,t); `):


elif nargs=1 and args[1]=pijtF then
print(`pijtF(i,j,t): A FAST version of pijt(i,j,t) (q.v.) using the second-order recurrence supplied by the Zeilberger algorithm.`):
print(`For example, try: `):
print(`pijtF(10,20,t); `):

elif nargs=1 and args[1]=qkt then
print(`qkt(k,t):q_k(t) defined in Eq. (16) of the paper, using the original definition of pij(i,j,t). SLOW. Try:`):
print(`qkt(10,t); `):

elif nargs=1 and args[1]=qktF then
print(`qktF(k,t):same as qkt(k,t) but faster, using pijtF(i,j,t) rather than the original pijt(i,j,t). Try:`):
print(`qktF(10,t); `):


elif nargs=1 and args[1]=qktFF then
print(`qktFF(k,t):same as qkt(k,t) but even faster, using  the third-order recurrence. Try:`):
print(`qktFF(10,t); `):

elif nargs=1 and args[1]=qktPE then
print(`qktPE(k,t): q_k'(t) acording to the explicit expression in the paper`):


elif nargs=1 and args[1]=qktPF then
print(`qktPF(k,t): q_k'(t) using the recurrence found by the Zeilberger algorithm`):


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

fi:


end:

ez:=proc(): print(` pijt(i,j,t), pijtF(i,j,t) , qkt(k,t), qktF(k,t), qktPE(k,t) , qktPF(k,t)  `): end:


#opeI:=(2*i*t^2-3*t^2-i+j+1)/t/(t+1)/(i-1)/Si-(t-1)*(i-2)/(t+1)/(i-1)/Si^2:
#opeJ:=(2*j*t^2-3*t^2+i-j+1)/t/(t+1)/(j-1)/Sj

#pijt(i,j,t): p_{ij}(t): given in Eq. 15 of Snir et. al.
pijt:=proc(i,j,t) local L:

normal(
add((i+j-L-1)!/(i-L)!/(j-L)!/(L-1)!*(1-t^2)^(L-1)*t^(i+j-2*L),L=1..min(i,j))/
(1+t)^(i+j-1)):
end:


#pijtF(i,j,t): Same as p_{ij}(t) but using the recurrence given by the Zeilberger algorithm.
#Try
#pijFt(10,10,t);
pijtF:=proc(i,j,t) option remember:
if j=0 then
 0:
elif j=1 then
if i=0 then
 0:
else
t^(i-1)/(1+t)^i:
fi:

else
normal((2*j*t^2-3*t^2+i-j+1)/t/(t+1)/(j-1)*pijtF(i,j-1,t)-(t-1)*(j-2)/(t+1)/(j-1)*pijtF(i,j-2,t)):
fi:
end:




#qkt(k,t):q_k(t) defined in Eq. (16) of Snir et. al.
qkt:=proc(k,t) local i,j:
normal(add(add(pijt(i,j,t),i=1..k),j=1..k)/k):
end:


#qktF(k,t):like q_k(t) but using pijtF rather than pijt
qktF:=proc(k,t) local i,j:
normal(add(add(pijtF(i,j,t),i=1..k),j=1..k)/k):
end:


#qktFF(k,t): q_k(t) using the third-order recurrence:

qktFF:=proc(k,t) option remember:
if k<=3 then
 qktF(k,t):
else
normal(
(3*k*t^2+2*k*t-3*t^2+3*k-3)/(1+t)^2/k*qktFF(k-1,t)
-(3*k*t^2-2*k*t-6*t^2+3*k+6*t-6)/(1+t)^2/k*qktFF(k-2,t)
+(k-3)*(t-1)^2/(1+t)^2/k*qktFF(k-3,t)):
fi:
end:


#qktPE(k,t): q_k'(t) acording to the explicit expression in the paper
qktPE:=proc(k,t) local m:
normal(-add(binomial(k-1,m)*binomial(k,m)*t^(2*m),m=0..k-1)/(t+1)^(2*k)):
end:

#qktPF(k,t): q_k'(t) using the recurrence
qktPF:=proc(k,t) option remember:
if k<=2 then
normal(diff(qkt(k,t),t)):
else
normal(2*(2*k^2*t^2-4*k*t^2+2*k^2+2*t^2-4*k+1)/(t+1)^2/k/(2*k-3)*qktPF(k-1,t)
-(t-1)^2*(2*k-1)*(k-2)/(t+1)^2/k/(2*k-3)*qktPF(k-2,t)):
fi:

end: