print(``):
print(`This Maple program accompanies the paper`):
print(`"WZ-Style Certification and Sister Celine's Technique`):
print(`for Abel-Type Sums"`):
print(`by John Majewicz`);

print(`For a hypergeometric term F(n,k), the program "celine" returns an`):
print(`operator  of order ORDER, with coefficients of degree DEG,`):
print(`which annihilates the Abel-type`): 
print(`term F(n,k,r,s):=F(n,k)r(r+k)^{k-1}(n-k+s)^{n-k}.  If ORDER or DEG`):
print(`are too low, celine will return an appropriate message`):



celine:=proc(F,ORDER,DEG)

  local  testr, tests, var, EXPR, ope, JOE1, JOE2, eq, temp_eq, temp_var,

         temp_ope, opefull, jj, ll, i, j, ii:


testr:=2:  tests:=1:  JOE1:=ORDER:  JOE2:=DEG:

var:={a[0,0]}:

EXPR:=a[0,0]:  ope:=a[0,0]:


# At this stage of the procedure, celine will do three things.  First, it will
# create the set of coefficients for an operator ope annihilating the Abel-type
# form F(n,k,r,s).  It also creates the operator ope itself as well as the
# expression EXPR which if `large enough' will have coefficients wrt k equal
# to 0.


for i from 1 to JOE1 do

  for j from 0 to JOE2 do

    var:=var union {a[i,j]}:

    EXPR:=EXPR+a[i,j]*simplify( subs({n=n-i,k=k-j},F)/F )*(r+j)*(n-k+s)^(-i+j)
/r/(k+r)^j:

    ope:=ope+a[i,j]*N^i*K^j*S^(-i+j)*R^(-j):

  od:

od:  

EXPR:=normal(EXPR):  EXPR:=numer(EXPR):  EXPR:=expand(EXPR):

# eq will be the set of coefficients

eq:={}:  print(degree(EXPR,k)):

for ii from 0 to degree(EXPR,k) do

  eq:=eq union {coeff(EXPR,k,ii)=0}:

od:


# It's much quicker for Maple to find solutions {a[i,j]} for specific values
# of r, s and n.  Moreover, if no solution exists for specific values of r, 
# s and n, then a general solution doesn't exist.  Therefore, before proceeding,
# celine checks to see if there is a particular solution.  If a `non-zero'
# solution exists, celine goes on to find the general solution.  Otherwise, it
# stops.


temp_eq:=subs({r=testr,s=tests,n=1},eq):

lprint(`temp_q has`,nops(temp_eq),`equations to be solved.`);

temp_var:=solve(temp_eq,var):

temp_ope:=subs(temp_var,ope): 

if temp_ope=0 then 

  RETURN(`No recurrence of order`,JOE1,`degree`,JOE2); 

fi:

print(`Now I'm trying to solve the general set of`,nops(eq),`equations!`);

var:=solve(eq,var);

ope:=subs(var,ope):

opefull:=ope;

ope:=subs(K=1,opefull):


# If there are infinitely many solutions a[i,j], celine will return the one 
# with all such a[i,j]=1.  That is the purpose of the next double loop.


for jj from 1 to JOE1 do

  for ll from 0 to JOE2 do
    
    ope:=subs(a[jj,ll]=1,ope):

  od:

od: 

RETURN(ope):

lprint(ope):

end: