#File elim_recurrence_oper.maple
 
#Warning; This is a pseudo-algorithm. It may give up because
#it takes too long
 
#elim:A Maple program for eliminating opertors by the naive approach
#The operators are P(N,K,n,k) and Q(N,K,n,k)
#IMPORTANT: A typical monomial n^i1 N^j1 K^j2 k^i2. Note that
#the K comes before the k but the n comes before the N
#To use it for your own opertors P and Q change the input
#part at the very bottom
 
#To run it, call this file 'elim', adjust
#the output at the bottom and type maple -q<elim>outelim
#To read the output do 'more outelim'
 
#elim(P,Q,ordN,N,n,K,k,para) finds operators A(N,n,K,k)
#B(N,n,K,k) and S(N,n) s.t. AP+BQ=S+(K-1)GARBAGE
#P and Q are the opertors, OrdNS is the expected order of N in  S(N,n)
#n and k are the "variables", k is the one to be eliminated
#N and K are the shifts in n and k respectively. para is a list
#of parameters featuring in P and Q are not n or k
 
 
#-------Description of subroutines followes
#mult is a program that multiplies two difference operators of the form
#Sum a_{i2,i3,i4}(n) N^i2 * K^i3 * k^i4
#times Sum b_{i1,i2,i3,i4} * n^i1 * N^i2 * K^i3 * k^i4
#multmtm multiples a monomial by a monomial by the rule
#a(n)*N^i*K^j*k^r *( n^a *N^b*K^c* k^d)=
#a(n)*(n+i)^a*N^(i+b)*K^(j+c)* (k-c)^r * k^d
#multopetm multiplies an operator by a monomial
#genpol(N,K,k,ordN,ordK,degk,a) generates a generic polynomial
#indexed by a[i] in N,K,k, of degrees ordN,ordK,degk respectively
#elim eliminates k mod (K-1) between P and Q
#khasel takes a set of assigenemtns {a[1]=a[1],...} and weeds out those
#that are a[1]=a[1]
 
yafe:=proc(oper,N)
local i,oper1,tosaa:
 
tosaa:=0:
 
oper1:=expand(oper):
 
for i from 0 to degree(oper1,N) do
 tosaa:=tosaa+factor(coeff(oper1,N,i))*N^i:
od:
 
tosaa:
 
end:
 
 
khasel:=proc(var)
local kv,i,mu:
 
kv:={}:
for i from 1 to nops(var) do
 mu:=op(i,var):
 
 if op(1,mu)=op(2,mu) then
  kv:=kv union {op(1,mu)}:
 fi:
 
od:
 
 
kv:
 
end:
 
 
 
 
genpol1:=proc(N,ordN,a)
local gu,iN,var,co:
gu:=1:
co:=1:
 
var:={}:
 
for iN from 1 to ordN do
   var:=var union {a[co]}:   
   gu:=gu+a[co]*N^iN:
   co:=co+1:
od:
 
gu,var:
 
end:
 
genpol3:=proc(N,K,k,ordN,ordK,degk,a)
local gu,iN,iK,ik,var,co:
gu:=0:
co:=1:
 
var:={}:
 
for iN from 0 to ordN do
 for iK from 0 to ordK do
  for ik from 0 to degk do
   var:=var union {a[co]}:   
   gu:=gu+a[co]*N^iN*K^iK*k^ik:
   co:=co+1:
  od:
 od:
od:
 
gu,var:
 
end:
 
 
multmtm:=proc(mono1,mono2,N,n,K,k)
local i,j,r,a,b,c,d,mek1,mek2:
 
i:=degree(mono1,N):
j:=degree(mono1,K):
r:=degree(mono1,k):
 
mek1:=normal(mono1/N^i/K^j/k^r):
 
a:=degree(mono2,n):
b:=degree(mono2,N):
c:=degree(mono2,K):
d:=degree(mono2,k):
 
mek2:=normal(mono2/n^a/N^b/K^c/k^d):
expand(mek1*mek2*(n+i)^a*N^(i+b)*K^(j+c)*(k-c)^r*k^d):
 
end:
 
 
 
 
multopetm:=proc(oper1,mono2,N,n,K,k)
local gu,i,mono1,oper:
 
oper:=expand(oper1):
if not type(oper,`+`) then
 RETURN(multmtm(oper,mono2,N,n,K,k)):
fi:
 
gu:=0:
 
for i from 1 to nops(oper) do
  mono1:=op(i,oper):
  gu:=gu+multmtm(mono1,mono2,N,n,K,k):
od:
 
gu:
 
end:
 
 
mult:=proc(oper1a,oper2a,N,n,K,k)
local mono2,gu,i,oper1,oper2:
 
oper1:=expand(oper1a):
oper2:=expand(oper2a):
 
if not type(oper2,`+`) then
  RETURN(multopetm(oper1,oper2,N,n,K,k)):
fi:
 
gu:=0:
 
for i from 1 to nops(oper2) do
 mono2:=op(i,oper2):
 gu:=gu+multopetm(oper1,mono2,N,n,K,k):
od:
 
gu:
 
end:
 
 
elimbound:=proc(P,Q,ordNS,ordN,ordK,degk,N,n,K,k,para)
local A,B,a,b,S,s,eq,var,gu,DEGguk,DEGguN,mu,eqspec,HALEV, 
SORRY, gu1, i, iN, ik,kv:
 
var:={}:
 
A:=genpol3(N,K,k,ordN,ordK,degk,a):
 
var:=var union A[2]:
 
A:=A[1]:
 
B:=genpol3(N,K,k,ordN,ordK,degk,b):
 
var:=var union B[2]:
 
B:=B[1]:
 
S:=genpol1(N,ordNS,s):
 
var:=var union S[2]:
 
S:=S[1]:
 
gu:=mult(A,P,N,n,K,k)+mult(B,Q,N,n,K,k)+S:
gu1:=subs(K=1,gu):
gu1:=expand(gu1):
 
eq:={}:
 
DEGguk:=degree(gu1,k):
 
 
for ik from 0 to DEGguk do
 mu:=coeff(gu1,k,ik):
 
  DEGguN:=degree(mu,N):
 
  for iN from 0 to DEGguN do
   eq:=eq union {coeff(mu,N,iN)=0}:
  od:
 
od:
 
eqspec:=subs(n=3,eq):
 
for i from 1 to nops(para) do
 eqspec:=subs(op(i,para)=i^3+2,eqspec):
od:
 
HALEV:=solve(eqspec,var):
 
 
if HALEV=NULL then
 RETURN(`SORRY`):
else
 print(`Looks good, please stand by`):
fi:
 
var:=solve(eq,var):
 
 
if var=NULL then
  RETURN(`SORRY`):
 
fi:
 
kv:=khasel(var):
 
A:=subs(var,A):
 
B:=subs(var,B):
 
S:=subs(var,S):
 
 
for i from 1 to nops(kv) do
 A:=subs(op(i,kv)=0,A):
 B:=subs(op(i,kv)=0,B):
 S:=subs(op(i,kv)=0,S):
od:
 
S:=yafe(S,N):
A:=yafe(A,N):
B:=yafe(B,N):
S,A,B:
 
end:
 
 
elim:=proc(P,Q,ordNS,N,n,K,k,para)
local ordN,ordK,degk,gu,SORRY,P1,Q1,maxordN,maxordK,maxdegk:
 
P1:=expand(P):
Q1:=expand(Q):
 
maxordN:=max(degree(P1,N),degree(Q1,N))+2:
maxordK:=max(degree(P1,K),degree(Q1,K))+2:
maxdegk:=max(degree(P1,k),degree(Q1,k)):
gu:=SORRY:
 
for degk from 0 to maxdegk do
 for ordN from 0 to maxordN do
  for ordK from 0 to maxordK do
 
   
    gu:=elimbound(P1,Q1,ordNS,ordN,ordK,degk,N,n,K,k,para):
 
     if nops([gu])>1 then
      RETURN(gu):
     fi:
  
  od:
 od:
od:
 
lprint(`No operator of order`,ordNS,`in`,N,`seems to exist`):
lprint(`try to raise the order`):
 
SORRY:
 
end:
 
#findop finds the first-order operator annihilating a given 
#closed-form function F, w.r.t. to the variable n, where N is 
#the corresponding shift in n
 
findop:=proc(F,n,N)
local gu:
 
gu:=simplify(expand(subs(n=n+1,F)/F)):
 
denom(gu)*N-numer(gu):
 
end:
 
findrec:=proc(F,n,k,N,K,ordNS,para)
local P,Q:
 
P:=findop(F,n,N):
Q:=findop(F,k,K):
Q:=K*subs(k=k-1,coeff(Q,K,1))+coeff(Q,K,0):
 
elim(P,Q,ordNS,N,n,K,k,para):
end:
 
 
#----Begin Input------------------
#This example proves Dixon's theorem
 
F:=(-1)^k/(n+k)!/(n-k)!/(a+k)!/(a-k)!/(b+k)!/(b-k)!:
ordNS:=1;
para:=[b,c];
findrec(F,n,k,N,K,ordNS,para):
lprint(");
quit;
 
#--------end program file--------------------------