#Procedure rec2ortho (version 1b3) for Maple V, release 4 or 5.
#Copyright 1996, 1997, 1998, 2005 by Tom H. Koornwinder and Rene' F. Swarttouw.
#email thk@science.uva.nl
#Input: coefficients of three term recurrence relation.
#Returns: explicit system of orthogonal polynomials satisfying these relations,
#provided the system is in the Askey scheme up to two-parameter level and has
#positive orthogonality measure, otherwise FAIL
#Last modified: 19 March 2005.
#
`rec2ortho/addprop`:=proc() local i,sequence;
option `Copyright 1996,1997,1998, 2005 by Tom H Koornwinder and Rene' F Swarttouw`;
sequence:=NULL:
for i from 1 to nargs/2
do
  if is(args[2*i-1],args[2*i])
  then
  elif not coulditbe(args[2*i-1],args[2*i])
  then
    RETURN(FAIL)
  else
    sequence:=sequence,args[2*i-1],args[2*i]
  fi
od:
if nops([sequence])>0
then
  traperror(additionally(sequence)):
  if assigned(lasterror)
  then
    print(lasterror);
    RETURN(FAIL)
  fi
fi:
RETURN(true)
end:

rec2ortho:=proc(an,bn,cn,A0,B0,C1)
local Anmin1,C,Fn,Gn,Gn0,Gnd,H,Hn,K,M,N,P,a,a0,ak,alpha,aminb,aplusb,
b,b0,beta,c1,cc,d,dFn,dGn,discri,dfn,dgn,fn,gn,k,kn,lambda,mu0,mun,nu1,nun,phi,
pp,prefn,procfail,s,sol,t,u,v; global n,x;
option `Copyright 1996,1997,1998, 2005 by Tom H Koornwinder and Rene' F Swarttouw`;

if nargs>6
then
  print(`Warning:`):
  print(`only the first six arguments of rec2ortho will be taken into account`)
fi:
procfail:=`not in Askey scheme or no positive orthogonality measure`:
if nargs>=4 and has(A0,n)
then
  print(`fourth argument should not depend on n`);
  RETURN(FAIL)
fi:
if nargs>=5 and has(B0,n)
then
  print(`fifth argument should not depend on n`);
  RETURN(FAIL)
fi:
if nargs>=6 and has(C1,n)
then
  print(`sixth argument should not depend on n`); RETURN(FAIL)
fi:
if nargs<6
then
  c1:=simplify(subs(n=1,simplify(cn)))
else
  c1:=C1
fi:
if nargs<5
then
  b0:=simplify(subs(n=0,simplify(bn)))
else
  b0:=B0
fi:
if nargs<4
then
  a0:=simplify(subs(n=0,simplify(an)))
else
  a0:=A0
fi:
if cn = 0 or c1=0
then
  print(procfail); RETURN(FAIL)
fi:
if an=0 or a0=0
then
  RETURN(`p(n,x) not of degree n`)
fi:
ak := subs(n=k,an);
kn :=1/product(ak,'k'=1..n-1)/a0;
#ak := subs(n=n-k,an);
#kn :=simplify(1/product(ak,'k'=1..n-1)/a0);
Anmin1 := subs(n=n-1,an);
nun := simplify(cn*Anmin1);
nu1:=c1*a0;
mun := bn;
mu0:=b0;
if not type(nun,ratpoly(anything,n)) or not type(mun,ratpoly(anything,n))
then
  print(procfail); RETURN(FAIL)
fi:
if interface(showassumed)<>2
then
  interface(showassumed=2)
fi:
Fn := collect(numer(nun),n);
Gn := collect(denom(nun),n);
dFn := degree(Fn,n);
dGn := degree(Gn,n);
fn := collect(numer(mun),n);
gn := collect(denom(mun),n);
dfn := degree(fn,n);
dgn := degree(gn,n);
if dGn=0
then
  nun:=collect(nun,n)
fi:
if dgn=0
then
  mun:=collect(mun,n)
fi:
if dFn-dGn < 0 or dFn-dGn=3 or dFn-dGn>4
then
  print(procfail); RETURN(FAIL)
fi:
if dFn-dGn = 0
then                      #Jacobi
  if dFn=1 or dFn=3 or dFn>4
  then
    print(procfail); RETURN(FAIL)
  fi:
  if simplify(subs(n=1,nun)-nu1)<>0 and simplify(subs(n=0,mun)-mu0)<>0
  then
    print(procfail); RETURN(FAIL)
  fi:
  if dFn=0
  then
    if has(mun,n)
    then
      print(procfail); RETURN(FAIL)
    fi:
    sol:=solve(x^2=1/nun/4,x):
    u:=sol[1]:
    v:=-u*mun;
    if nu1-nun=0 and mu0-mun=0
    then
      alpha:=1/2;
      beta:=1/2;
    elif nu1-nun=0 and mu0-mun<>0
    then
      alpha:=-v-u*mu0;
      if alpha-1/2<>0 and alpha+1/2<>0
      then
        print(procfail); RETURN(FAIL)
      fi:
      beta:=-alpha
    elif nu1-2*nun=0
    then
      alpha:=-1/2;
      beta:=-1/2
    else
      print(procfail); RETURN(FAIL)
    fi
  else            #dFn=2 or 4
    if dFn=2
    then
      s:=simplify(coeff(Fn,n,1)/coeff(Fn,n,2));
      if simplify(coeff(Gn,n,1)/coeff(Gn,n,2)-s)<>0
      then
        print(procfail); RETURN(FAIL)
      fi:
      Gn0:=simplify(coeff(Gn,n,0)/coeff(Gn,n,2));
      if simplify(coeff(Fn,n,0))=0
      then
        Gnd:=simplify(s^2-4*Gn0);
        if Gnd<>0 and Gnd-1<>0
        then
          print(procfail); RETURN(FAIL)
        fi:
        a:=(s+1-Gnd)/2;
        b:=(s-1+Gnd)/2
      else
        if (s+1<>0 and s<>0 and s-1<>0) or simplify(s^2-2*Gn0-1/2)<>0
        then
          print(procfail); RETURN(FAIL)
        fi:
        t:=coeff(Fn,n,0)/coeff(Fn,n,2);
        sol:=solve(x^2-s*x+t,x);
        a:=sol[1]: b:=sol[2]:
      fi:
      if is(a+1<=0) or is(b+1<=0)
      then
        print(procfail); RETURN(FAIL)
      fi:
      sol:=solve(x^2=simplify(coeff(Gn,n,2)/coeff(Fn,n,2)/4),x):
      u:=sol[1]:
    else                                      #dFn=4
      if coeff(Fn,n,0)<>0 or coeff(Fn,n,3)=0
      then
        print(procfail); RETURN(FAIL)
      fi:
      s:=coeff(Fn,n,3)/coeff(Fn,n,4);
      if simplify(subs(n=-s/2,Fn))<>0
      then
        print(procfail); RETURN(FAIL)
      fi:
      t:=2*coeff(Fn,n,1)/coeff(Fn,n,4)/s;
      sol:=solve(x^2-s*x/2+t,x);
      a:=sol[1]: b:=sol[2]:
      if is(a<=-1) or is(b<=-1) or is(a+b<=-2)
      then
        print(procfail); RETURN(FAIL)
      fi:
      sol:=solve(x^2=simplify(4*n*(n+a)*(n+b)*(n+a+b)/
        nun/(2*n+a+b-1)/(2*n+a+b)^2/(2*n+a+b+1)),x):
      u:=sol[1]:
      if depends(u,n)
      then
        print(procfail); RETURN(FAIL)
      fi
    fi:
    Hn:=-simplify(u*mun*(2*n+a+b)*(2*n+a+b+2));
    if not type(Hn,polynom(anything,n))
    then
      print(procfail); RETURN(FAIL)
    fi:
    Hn:=collect(Hn,n);
    if degree(Hn,n)>2
    then
      print(procfail); RETURN(FAIL)
    fi:
    v:=coeff(Hn,n,2)/4;
    d:=simplify(Hn-v*(2*n+a+b)*(2*n+a+b+2));
    if depends(d,n)
    then
      print(procfail); RETURN(FAIL)
    fi:
    if simplify(a+b)=0 and d=0
    then
      alpha:=simplify(-v-u*mu0);
      beta:=-alpha;
      if simplify(alpha^2-a^2)<>0
      then
        print(procfail); RETURN(FAIL)
      fi
    elif simplify(d-a^2+b^2)=0
    then
      alpha:=a;
      beta:=b
    elif simplify(d-b^2+a^2)=0
    then
      alpha:=b;
      beta:=a
    else
      print(procfail); RETURN(FAIL)
    fi:
    if simplify(alpha+beta+1)=0
    then
      if simplify(nu1+2*(alpha+1)*alpha/u^2)<>0
      then
        print(procfail); RETURN(FAIL)
      fi
    elif simplify(subs(n=1,nun)-nu1)<>0
    then
      print(procfail); RETURN(FAIL)
    fi:
    if simplify(alpha+beta)<>0 and simplify(subs(n=0,mun)-mu0)<>0
    then
      print(procfail); RETURN(FAIL)
    fi:
  fi:
  prefn:=simplify(kn*n!*2^n/pochhammer(n+alpha+beta+1,n)/u^n);
  aplusb:=simplify(alpha+beta):
  aminb:=simplify(alpha-beta):
  if is(aplusb,constant)
  then
    if not is(aplusb>-2)
    then
      print(procfail); RETURN(FAIL)
    fi
  fi:
  if is(aminb,constant)
  then
    if not is(aminb,real)
    then
      print(procfail); RETURN(FAIL)
    fi
  fi:
  if `rec2ortho/addprop`(alpha,RealRange(Open(-1),infinity),
    beta,RealRange(Open(-1),infinity),
    u,AndProp(real,Non(0)),v,real)
  then
  else
    print(procfail); RETURN(FAIL)
  fi:
  alpha:=eval(alpha):
  beta:=eval(beta):
  u:=eval(u):
  v:=eval(v):
  print(`Jacobi `);
  RETURN(prefn * P(n,u*x+v,alpha,beta))
fi:
if dFn-dGn = 1
then           #Hermite or Charlier
  if simplify(subs(n=1,nun)-nu1)<>0 or simplify(subs(n=0,mun)-mu0)<>0
  then
    print(procfail); RETURN(FAIL)
  fi:
  if dGn > 0 
  then
    print(procfail); RETURN(FAIL)
  fi:
  if coeff(nun,n,0) <> 0 or not type(mun,polynom(anything,n))
  then
    print(procfail); RETURN(FAIL)
  fi:
  if is(coeff(nun,n,1)<0)
  then
    print(procfail); RETURN(FAIL)
  fi:
  if degree(mun,n) > 1
  then
    print(procfail); RETURN(FAIL)
  fi:
  if coeff(mun,n,1) = 0
  then #Hermite
    sol:=solve(x^2=simplify(1/(2*coeff(nun,n,1))),x):
    u:=sol[1]:
    v:=simplify(-coeff(mun,n,0)*u):
    prefn := simplify(kn/(2*u)^n);
    if `rec2ortho/addprop`(u,AndProp(real,Non(0)),v,real)
    then
    else
      print(procfail); RETURN(FAIL)
    fi:
    u:=eval(u); v:=eval(v);
    print(`Hermite`);
    RETURN(prefn*H(n,u*x+v))
  else #Charlier
    u:=simplify(1/coeff(mun,n,1)):
    a:=simplify(coeff(nun,n,1)/(coeff(mun,n,1))^2);
    v:=simplify(-coeff(mun,n,0)/coeff(mun,n,1)+
      coeff(nun,n,1)/(coeff(mun,n,1))^2);
    prefn := simplify(kn*(-a/u)^n);
    if `rec2ortho/addprop`(u,AndProp(real,Non(0)),a,RealRange(Open(0),infinity),
      v,real)
    then
    else
      print(procfail); RETURN(FAIL)
    fi:
    u:=eval(u): v:=eval(v); a:=eval(a);
    print(`Charlier`);
    RETURN(prefn*C(n,u*x+v,a))
  fi
fi:
if dFn-dGn = 2 #Hahn, Continuous Hahn, Laguerre, Meixner-Pollaczek, Meixner or 
               #Krawtchouk
then
  if dGn>0        #Hahn or Continuous Hahn
  then
    RETURN(`possibly Hahn or Continuous Hahn`)
  fi:
#dGn=0
  if coeff(nun,n,0) <> 0
  then
    print(`possibly Hahn or Continuous Hahn`); RETURN(FAIL)
  fi:      
#coeff(nun,n,0)=0  
  if dgn > 0 
  then
    print(procfail); RETURN(FAIL)
  fi:
  if dfn>1
  then
    print(procfail); RETURN(FAIL)
  fi:
#dgn=0, dfn=0 or 1
#Meixner-Pollaczek, Lagurerre, Meixner or Krawtchouk
  if simplify(subs(n=1,nun)-nu1)<>0 or simplify(subs(n=0,mun)-mu0)<>0
  then
    print(procfail); RETURN(FAIL)
  fi:
  discri:=simplify((coeff(mun,n,1))^2-4*coeff(nun,n,2));
  if not is(discri,real)
  then
    print(procfail); RETURN(FAIL)
  elif discri = 0
  then          #Laguerre
    a:=simplify(4*coeff(nun,n,1)/(coeff(mun,n,1))^2):
    u:=simplify(2/coeff(mun,n,1)):
    v:=simplify(4*coeff(nun,n,1)/(coeff(mun,n,1))^2+1-
    2*coeff(mun,n,0)/coeff(mun,n,1)):
    prefn := simplify(kn*(-1)^n*n!/u^n);
    if `rec2ortho/addprop`(a,RealRange(Open(-1),infinity),u,AndProp(real,Non(0)),
      v,real)
    then
    else
      print(procfail); RETURN(FAIL)
    fi:
    a:=eval(a):
    u:=eval(u):
    v:=eval(v):
    print(`Laguerre`);
    RETURN(prefn * L(n,u*x+v,a))
  elif is(discri<0)
  then          #Meixner-Pollaczek
    lambda := simplify((coeff(nun,n,1)+coeff(nun,n,2))/2/coeff(nun,n,2));
    sol:=solve(x^2=simplify(1/((4*coeff(nun,n,2)-(coeff(mun,n,1))^2))),x):
    u:=sol[1]:
    phi := simplify(arccot(-u*coeff(mun,n,1)));
    v := simplify(-u*coeff(mun,n,0)+lambda*u*coeff(mun,n,1));
    prefn := simplify(kn*n!/(2^n)/(sin(phi))^n*u^(-n));
    if `rec2ortho/addprop`(u,AndProp(real,Non(0)),lambda,RealRange(Open(0),
      infinity),
    phi,RealRange(Open(0),Open(Pi)),v,real)
    then
    else
      print(procfail); RETURN(FAIL)
    fi:
    u:=eval(u):
    lambda:=eval(lambda):
    phi:=eval(phi):
    v:=eval(v):
    print(`Meixner-Pollaczek`):
    RETURN(prefn*P(n,u*x+v, lambda,phi))
  elif is(discri>0)
  then          #Meixner or Krawtchouk      
    beta :=simplify((coeff(nun,n,1)+coeff(nun,n,2))/coeff(nun,n,2));
    if beta=0
    then
      print(procfail); RETURN(FAIL)
    elif not is(beta,real)
    then
      print(procfail); RETURN(FAIL)
    elif is(beta<0)
    then            #Krawtchouk
      N:=-beta;
      sol:=solve(x^2=simplify(1/((coeff(mun,n,1))^2-4*coeff(nun,n,2))),x);
      u:=sol[1]:
      pp := simplify(1/2 - u*coeff(mun,n,1)/2);
      v := simplify(-u*coeff(mun,n,0)+pp*N);
      prefn := simplify(kn*(-N)*product(k-N,'k'=1..n-1)*pp^n/u^n);
      if `rec2ortho/addprop`(u,AndProp(real,Non(0)),N,integer,
        pp,RealRange(Open(0),Open(1)),v,real)
      then
      else
        print(procfail); RETURN(FAIL)
      fi:
      u:=eval(u):
      N:=eval(N):
      pp:=eval(pp):
      v:=eval(v):
      if is(pp<=0) or is(pp>=1)
      then
        print(procfail); RETURN(FAIL)
      fi:
      print(`Krawtchouk`);
      RETURN(prefn*K(n,u*x+v,pp,N))
    elif is(beta>0)
    then            #Meixner
      sol:=solve(x^2=simplify(1/((coeff(mun,n,1))^2-4*coeff(nun,n,2))),x):
      u:=sol[1]:
      cc := simplify((u*coeff(mun,n,1)-1)/(u*coeff(mun,n,1)+1));
      v := simplify(-u*coeff(mun,n,0)+beta*(u*coeff(mun,n,1)-1)/2);
      prefn := simplify(kn*pochhammer(beta,n)*(cc/(cc-1))^n/u^n);
      if `rec2ortho/addprop`(u,AndProp(real,Non(0)),v,real,
        cc,OrProp(RealRange(Open(0),Open(1)),RealRange(Open(1),infinity)))
      then
      else
        print(procfail); RETURN(FAIL)
      fi:
      u:=eval(u):
      cc:=eval(cc):
      v:=eval(v):
      print(`Meixner`):
      RETURN(prefn*M(n,u*x+v,beta,cc))
    fi:
    print(`possibly Meixner or Krawtchouk`, `assume more about the sign of`);
    print(beta); RETURN(FAIL)
  fi:
  print(`possibly Meixner or Krawtchouk or Meixner-Pollaczek`,
 `assume more about the sign of`);
  print(discri); RETURN(FAIL)
fi:
if dFn-dGn = 4
then
  print(`possibly Hahn, Continuous Dual Hahn, Racah, Wilson`);; RETURN(FAIL)
fi
end;