# multsum.mpl
# (C) Torsten Sprenger, University of Kassel, 2003-2005
#
# Torsten Sprenger
# sprenger@mathematik.uni-kassel.de
# http://www.mathematik.uni-kassel.de/~sprenger


_EnvSmax:=10:
_EnvC:=1:
_EnvP:=1:
_EnvSolve:=SolveTools[Linear]:
_EnvFilter:=0:

########################################################### 
#                                                         #
#  'check_hyper' überprüft, ob ein eingegebener Term t    #
#  hypergeometrisch ist. Ist das nicht der Fall, so wird  #
#  ein Fehler ausgegeben, andernfalls passiert nichts.    #
#                                                         #
###########################################################

check_hyper:=proc(t,k,n)
local i,var,hyper;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  var:=[n,op(k)];
  hyper:=andmap(evalb,
                {seq(type(ratio(t,var[i]),ratpoly(anything,var[i])),
                     i=1..nops(var))
                });
  if not(hyper) 
    then ERROR(`Term is not hypergeometric`)
  end if;
end proc:

checkhyper:=proc(t,k,n)
local check;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  check:=check_hyper(t,k,n);
  if check=NULL
    then return(true) 
    else return(false)
  end if;
end proc:


########################################################### 
#                                                         #
#  'check_proper' überprüft, ob ein eingegebener Term t   #
#  zulässig (proper) ist. Ist das nicht der Fall, so      #
#  wird ein Fehler ausgegeben, andernfalls wird das       #
#  Polynom (pol) und der GAMMA-Ausdruck (gam) des PHT's   #
#  (proper hypergeometric term) zurückgegeben.            #
#                                                         #
###########################################################

check_proper:=proc(t,k,n)
local i,x,term,kk,pol_den,pol_den_,den_set,var,rules,gam_den,gam_num,gam,pol;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  check_hyper(t,k,n);
  kk:=[op(k)];
  var:=[n,op(k)];
  rules:=[seq((x::freeof(var))^var[i]=1,i=1..nops(kk)+1)];
  term:=applyrule(rules,simpcomb(t));
  gam_num:=numer(term);
  pol_den:=factor(applyrule(GAMMA(x::anything)=1,denom(term)));
  pol_den_:=1;
  if type(pol_den,`*`)
    then 
      for i from 1 to nops(pol_den) do
        if type(op(i,pol_den),freeof(var))
          then gam_num:=gam_num/op(i,pol_den); 
               pol_den_:=pol_den_*op(i,pol_den); 
        end if;
      end do;
  end if;
  pol_den:=pol_den/pol_den_;
  gam_den:=1;
  if not(type(pol_den,freeof(var))) 
    then den_set:=convert(factor(pol_den),`multiset`);
      for i from 1 to nops(den_set) do
        if degree(den_set[i,1],{op(var)})=1
          then gam_den:=gam_den*(GAMMA(den_set[i,1]+1)/GAMMA(den_set[i,1]))^(den_set[i,2])
          else ERROR(`Hypergeometric term is not proper`)
        end if 
      end do 
    else pol_den:=1;
  end if;
  gam_den:=gam_den*(denom(term)/pol_den);
  gam:=gam_num/gam_den;
  pol:=applyrule(GAMMA(x::anything)=1,gam);
  gam:=normal(gam/pol);

  return([pol,gam]);    
end proc:

checkproper:=proc(t,k,n)
local check;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  try check:=check_proper(t,k,n);
  catch: return(false);
  end try;
  return(true);
end proc:


########################################################### 
#                                                         #
#  'convertgam' extrahiert die Argumente der GAMMA-       #
#  Funktionen des GAMMA-Ausdrucks und überführt sie in    #
#  eine Liste für die Zähler-Ausdrücke (num_fac_ex) und   #
#  in eine Liste für die Nenner-Ausdrücke (den_fac_ex).   #
#                                                         #
###########################################################

convertgam:=proc(gam)
local i,j,num,den,upper,lower,num_fac_ex,den_fac_ex;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  num:=convert(numer(gam),`multiset`);
  den:=convert(denom(gam),`multiset`);
  upper:=applyrule(GAMMA(x::anything)=x,num); 
  lower:=applyrule(GAMMA(x::anything)=x,den);

  num_fac_ex:=[];
  den_fac_ex:=[];
  
  if not(type(num,freeof(GAMMA))) then 
    for i from 1 to nops(upper) do
      num_fac_ex:=[op(num_fac_ex),seq(upper[i,1],j=1..upper[i,2])]
    end do;
  end if;
  
  if not(type(den,freeof(GAMMA))) then
    for i from 1 to nops(lower) do
      den_fac_ex:=[op(den_fac_ex),seq(lower[i,1],j=1..lower[i,2])]
    end do;
  end if;
  
  return([num_fac_ex,den_fac_ex]);    
end proc:


########################################################### 
#                                                         #
#  'convertfacex' bestimmt                                #
#  die Koeffizienten der Zähler-Ausdrücke (a,b,c) und     #
#  die Koeffizienten der Nenner-Ausdrücke (u,v,w) und     #
#  gibt sie jeweils als Listen aus.                       #
#                                                         #
###########################################################

convertfacex:=proc(num_fac_ex,den_fac_ex,k,n)
local i,j,kk,pp,qq,a,b,c,u,v,w,r,var;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  kk:=[op(k)];
  var:=[n,op(k)];
  r:=nops(k);
  pp:=nops(num_fac_ex);
  qq:=nops(den_fac_ex);

  a:=[seq(coeff(num_fac_ex[i],n,1),i=1..pp)];
  b:=[seq([seq(coeff(num_fac_ex[i],var[j+1],1),i=1..pp)],j=1..r)];
  c:=[];
  for i from 1 to pp do
    if ldegree(num_fac_ex[i],var)=1 
      then c:=[op(c),0]
      else c:=[op(c),tcoeff(num_fac_ex[i],var)]
    end if
  end do;

  u:=[seq(coeff(den_fac_ex[i],n,1),i=1..qq)];
  v:=[seq([seq(coeff(den_fac_ex[i],var[j+1],1),i=1..qq)],j=1..r)];
  w:=[];
  for i from 1 to qq do
    if ldegree(den_fac_ex[i],var)=1 
      then w:=[op(w),0]
      else w:=[op(w),tcoeff(den_fac_ex[i],var)]
    end if
  end do; 
  
  if nargs>4 
    then return([a,b,u,v]) 
    else return([a,b,c,u,v,w])
  end if;
end proc:


########################################################### 
#                                                         #
#  'properterm' bestimmt die Bestandteile eines einge-    #
#  gebenen Terms t. Ausgegeben wird eine Liste mit drei   #
#  Elementen:                                             #
#  (1): Polynom (pol)                                     #
#  (2): die Koeffizienten der Zähler- und Nenner-         # 
#       Ausdrücke (a,b,c,u,v,w)                           #
#  (3): restliche Potenzen.                               #
#                                                         #
###########################################################

properterm:=proc(t,k,n)
local pol_gam,abc_uvw,var;
option remember, `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  var:=[n,op(k)];
  pol_gam:=check_proper(t,k,n);
  abc_uvw:=convertgam(pol_gam[2]);

  return([pol_gam[1],
          convertfacex(abc_uvw[1],abc_uvw[2],k,n),
          simpcomb(t/(pol_gam[1]*pol_gam[2]))]);  
end proc:

                                                     
###########################################################
#                                                         #
#  'facdiff' bestimmt die Fakultäts-Differenz             #
#  und gibt sie als [A_t,B_t] aus.                        #
#                                                         #
###########################################################

facdiff:=proc(t,k,n)
local i,j,r,propterm,abc_uvw,p,q,pp,qq,A,B;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  propterm:=properterm(t,k,n);
  r:=nops(k);
  pp:=nops(propterm[2,1]);
  qq:=nops(propterm[2,4]); 
  abc_uvw:=propterm[2];
  
  A:=0;
  for p from 1 to pp do 
    if {seq(abc_uvw[2,j,p],j=1..r)}<>{0}
      then A:=A+abc_uvw[1,p]
    end if
  end do;
  for q from 1 to qq do 
    if {seq(abc_uvw[5,j,q],j=1..r)}<>{0} 
      then A:=A-abc_uvw[4,q] 
    end if 
  end do;
   
  for i from 1 to r do
    B[i]:=0;
    for p from 1 to pp do
      B[i]:=B[i]+abc_uvw[2,i,p]
    end do; 
    for q from 1 to qq do
      B[i]:=B[i]-abc_uvw[5,i,q] 
    end do;
  end do;
   
  return([A,seq(B[i],i=1..r)]);
end proc: 


########################################################### 
#                                                         #
#  'structurefct' bestimmt die Parameter der Struktur-    #
#  funktionen eines PHT's t.                              #
#                                                         #
###########################################################

structurefct:=proc(t,k,n)
local i,j,r,propterm,abc_uvw,AB,p,q,pp,qq,gh_abc,gh_uvw,gh_AB;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  propterm:=properterm(t,k,n);
  r:=nops(k);
  pp:=nops(propterm[2,1]);
  qq:=nops(propterm[2,4]); 
  abc_uvw:=propterm[2];
  
  for p from 1 to pp do 
    gh_abc[p]:=[];
    if {seq(abc_uvw[2,j,p],j=1..r)}<>{0}
      then gh_abc[p]:=[op(gh_abc[p]),abc_uvw[1,p]/igcd(abc_uvw[1,p],seq(abc_uvw[2,j,p],j=1..r))];
        for i from 1 to r do
          gh_abc[p]:=[op(gh_abc[p]),abc_uvw[2,i,p]/igcd(abc_uvw[1,p],seq(abc_uvw[2,j,p],j=1..r))];
        end do 
    end if
  end do;
  
  for q from 1 to qq do 
    gh_uvw[q]:=[];
    if {seq(abc_uvw[5,j,q],j=1..r)}<>{0} 
      then gh_uvw[q]:=[op(gh_uvw[q]),-abc_uvw[4,q]/igcd(abc_uvw[4,q],seq(abc_uvw[5,j,q],j=1..r))];
        for i from 1 to r do
          gh_uvw[q]:=[op(gh_uvw[q]),-abc_uvw[5,i,q]/igcd(abc_uvw[4,q],seq(abc_uvw[5,j,q],j=1..r))];
        end do 
    end if
  end do;
  
  gh_AB:=NULL;
  AB:=facdiff(t,k,n);
  if {op(AB)}<>{0}
    then gh_AB:=[seq(-AB[j]/igcd(seq(AB[i],i=1..r+1)),j=1..r+1)]
  end if;
  
  return(remove(x->(x=[]),[seq(gh_abc[p],p=1..pp),
                           seq(gh_uvw[q],q=1..qq),
                           gh_AB]));
end proc:


########################################################### 
#                                                         #
#  'boundpts_Sstd' bestimmt Randpunkte der                #
#  Standard-Strukturmenge Sstd (bp) und gibt zusätzlich   #
#  die zugrunde liegenden Parameter der Struktur-         #
#  funktionen (sf) aus.                                   #
#                                                         #
###########################################################

boundpts_Sstd:=proc(t,k,n,IJ)
local i,j,num_IJ,sf,num_sf,bp;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  sf:=structurefct(t,k,n);
  num_IJ:=nops(IJ);
  num_sf:=nops(sf);
  for i from 1 to num_sf do
    for j from 1 to num_IJ do
      if type(sf[i,j],nonnegative) 
        then bp[i,j]:=IJ[j] 
        else bp[i,j]:=0 
      end if;
    end do;
  end do;

  return([sf,[seq([seq(bp[i,j],j=1..num_IJ)],i=1..num_sf)]]);
end proc:

boundptsSstd:=proc(t,k,n,IJ)
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  return(boundpts_Sstd(t,k,n,IJ)[2]);
end proc:


########################################################### 
#                                                         #
#  'boundpts_S' bestimmt Randpunkte einer                 #
#  Strukturmenge S (bp) und gibt zusätzlich               #
#  die zugrunde liegenden Parameter der Struktur-         #
#  funktionen (sf) aus.                                   #
#                                                         #
###########################################################

boundpts_S:=proc(t,k,n,S)
local i,j,num_S,sf,num_sf,bp;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  sf:=structurefct(t,k,n);
  num_S:=nops(S);
  num_sf:=nops(sf);
  
  for i from 1 to num_sf do
    bp[i]:=[S[1,1],S[1,2]];
    for j from 2 to num_S do
      if bp[i][1]*sf[i,1]+bp[i][2]*sf[i,2]<S[j,1]*sf[i,1]+S[j,2]*sf[i,2]
        then bp[i]:=[S[j,1],S[j,2]];
      end if;
    end do;
  end do;

  return([sf,[seq(bp[i],i=1..num_sf)]]);
end proc:


########################################################### 
#                                                         #
#  'inequat' erzeugt die zu den Strukturfunktionen und    #
#  Randpunkten der Standard-Strukturmenge Sstd zuge-      #
#  hörigen Ungleichungen, die die P-maximale Struktur-    #
#  menge Smax eingrenzen.                                 #
#                                                         #            
###########################################################

inequat:=proc(t,k,n,IJ,x)
local i,j,num_IJ,num_sf,sf_bp;
option remember, `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  sf_bp:=boundpts_Sstd(t,k,n,IJ);
  num_IJ:=nops(IJ);
  num_sf:=nops(sf_bp[1]);  

  return({seq(
           sum('sf_bp[1,i,j]*x[j]','j'=1..num_IJ)<=
           sum('sf_bp[1,i,j]*sf_bp[2,i,j]','j'=1..num_IJ),
          i=1..num_sf)});
end proc:


########################################################### 
#                                                         #
#  'isfinite' überprüft, ob die Menge S, die die          #
#  Standard-Strukturmenge Sstd umfasst, endlich und       #
#  somit die P-maximale Strukturmenge ist oder nicht.     #
#  Ausgabe: min_max_IJ, eine Liste der Grenzen des        #
#  kleinsten S umfassenden Rechtecks/Quaders (sofern S    #
#  endlich ist), andernfalls false.                       #
#                                                         #
###########################################################

isfinite:=proc(t,k,n,IJ)
local i,x,var,ineq,min_max_IJ;
option remember, `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  var:=[n,op(k)];
  ineq:=inequat(t,k,n,IJ,x);
  
  min_max_IJ:=[seq([ceil(subs(simplex[minimize](x[i],ineq),x[i])),
                   floor(subs(simplex[maximize](x[i],ineq),x[i]))],i=1..nops(var))];
  if not(type(min_max_IJ,freeof(x))) 
    then min_max_IJ:=false;
  end if;

  return(min_max_IJ);
end proc:


########################################################### 
#                                                         #
#  'Sstd' kann zum Erzeugen der Standard-Strukturmenge    #
#  Sstd, aber auch zur Erzeugung eines einfachen          #
#  kartesischen Produktes verwendet werden.               #
#                                                         #            
###########################################################

Sstd:=proc(IJ)
local i,num_IJ,set,Sstd_;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  num_IJ:=nops(IJ);
  if type(IJ,`listlist`)
    then set:=combinat[cartprod]([seq([$IJ[i,1]..IJ[i,2]],i=1..num_IJ)]);
    else set:=combinat[cartprod]([seq([$0..IJ[i]],i=1..num_IJ)]);
  end if;  

  Sstd_:={};
  while not set[finished] do 
    Sstd_:=Sstd_ union {set[nextvalue]()} 
  end do;

  return(Sstd_);
end proc:


########################################################### 
#                                                         #
#  'Scov' erzeugt aus einer vorgegebenen Struktur-        #
#  menge Sstd eine neue Strukturmenge Scov, in der        #
#  die P-maximale Strukturmenge Smax bezüglich Sstd       #
#  liegt. In Scov werden nur die Grenzen gespeichert.     #
#                                                         #
###########################################################

Scov:=proc(t,k,n,IJ)
local i,num_IJ,Scov_,min_max_IJ,extend;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  num_IJ:=nops(IJ);
  min_max_IJ:=isfinite(t,k,n,IJ);
  if type(min_max_IJ,`list`) 
    then Scov_:=min_max_IJ;
    else if 4<nargs then extend:=args[5] else extend:=_EnvSmax end if;
         Scov_:=[seq([-extend,extend],i=1..num_IJ)];
  end if;
  
  return(Scov_);
end proc:


########################################################### 
#                                                         #
#  'm_Smax' erzeugt die P-maximale Strukturmenge Smax,    #
#  die die Standard-Strukturmenge Sstd umfasst.           #
#  Ist Smax nicht endlich, so kann man mit einem          #
#  zusätzlichen Argument die Grösse der Menge Smax        #
#  steuern. Scov ist die P-maximale Menge umfassende      #
#  Menge, von der man sukzessive Elemente eliminiert, um  #
#  schliesslich Smax zu bekommen.                         #
#                                                         #
###########################################################

m_Smax:=proc(t,k,n,IJ)
local i,j,x,extend,m,num_IJ,num_ineq,Scov_,rules,Smax_,set,set_el,ineq;
option remember, `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  if 4<nargs then extend:=args[5] else extend:=NULL end if;
  num_IJ:=nops(IJ);
  Scov_:=Scov(t,k,n,IJ,extend);
    
  ineq:=inequat(t,k,n,IJ,x);
  num_ineq:=nops(ineq);
  set:=combinat[cartprod]([seq([$Scov_[i,1]..Scov_[i,2]],i=1..num_IJ)]);
  
  Smax_:={};
  while not set[finished] do 
    set_el:=set[nextvalue]();
    rules:={seq(x[m]=set_el[m],m=1..num_IJ)};
    if andmap(evalb,{seq(subs(rules,ineq[j]),j=1..num_ineq)}) 
      then Smax_:=Smax_ union {set_el};
    end if;
  end do;    
   
  return(Smax_);
end proc:


########################################################### 
#                                                         #
#  'drawsupport' zeichnet die nicht-verschwindenden       #
#  Ausdrücke eines zulässigen Terms als Punkte in ein     #
#  Koordinatensystem.                                     #
#                                                         #
###########################################################

binomial_:=proc(n,k)
  if type(k,negint)
    then return 0
    else return binomial(n,k)
  end if
end proc:

drawsupport:=proc(t,k,n,bd)
local i,j,m,pts,val,val_num,epsilon,undef_pts;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  pts:=[];
  undef_pts:=[];
  epsilon:=10^(-15);
  if nops(k)=1 then
    for i from 0 to bd do
      for j from -bd to bd do
        try
          val:=eval(subsindets(subs({n=i,k=j},t),'specfunc(anything,binomial)',f->subsop(0=binomial_,f)));
          val_num:=eval(subsindets(subs({n=i,k=j},numer(t)),'specfunc(anything,binomial)',f->subsop(0=binomial_,f)));
        catch:
          val:=evalf(subsindets(subs({n=i,k=j},t),'specfunc(anything,binomial)',f->subsop(0=binomial_,f)));
          val_num:=evalf(subsindets(subs({n=i,k=j},numer(t)),'specfunc(anything,binomial)',f->subsop(0=binomial_,f)));
        end try; 
        if hastype(val,undefined) or hastype(val_num,infinity)
          then undef_pts:=[op(undef_pts),[i,j]];
          else if not(hastype(val,name)) 
                 then if abs(val)>epsilon
                        then pts:=[op(pts),[i,j]]  
                      end if;
                 else pts:=[op(pts),[i,j]]
               end if;
        end if;      
      end do;
    end do;
    return(plots[display](seq(plottools[point](pts[i],symbol=circle,symbolsize=15,color=black),i=1..nops(pts)),
                          seq(plottools[point](undef_pts[i],symbol=circle,symbolsize=15,color=grey),i=1..nops(undef_pts)),
                          axes=normal,scaling=constrained,view=[0..bd,-bd..bd]));
  end if;

  if nops(k)=2 and nargs<5 then
    for i from 0 to bd do
      for j from -bd to bd do
        for m from -bd to bd do
          try
            val:=eval(subsindets(subs({n=i,k[1]=j,k[2]=m},t),'specfunc(anything,binomial)',f->subsop(0=binomial_,f)));
            val_num:=eval(subsindets(subs({n=i,k[1]=j,k[2]=m},numer(t)),'specfunc(anything,binomial)',f->subsop(0=binomial_,f)));
          catch:
            val:=evalf(subsindets(subs({n=i,k[1]=j,k[2]=m},t),'specfunc(anything,binomial)',f->subsop(0=binomial_,f)));
            val_num:=evalf(subsindets(subs({n=i,k[1]=j,k[2]=m},numer(t)),'specfunc(anything,binomial)',f->subsop(0=binomial_,f)));
          end try;
          if hastype(val,undefined) or hastype(val_num,infinity)
            then undef_pts:=[op(undef_pts),[i,j,m]];
            else if not(hastype(val,name)) 
                   then if abs(val)>epsilon
                          then pts:=[op(pts),[i,j,m]]  
                        end if;
                   else pts:=[op(pts),[i,j,m]]
                 end if;
          end if;    
        end do;
      end do;
    end do;
    return(plots[display3d](seq(plottools[point](pts[i],symbol=circle,symbolsize=15,color=black),i=1..nops(pts)),
                          seq(plottools[point](undef_pts[i],symbol=circle,symbolsize=15,color=grey),i=1..nops(undef_pts)),
                          axes=normal,scaling=constrained,view=[0..bd,-bd..bd,-bd..bd]));
  end if;

  if nops(k)=2 and nargs>4 then
    for i from -bd to bd do
      for j from -bd to bd do
        try
          val:=eval(subsindets(subs({n=args[5],k[1]=i,k[2]=j},t),'specfunc(anything,binomial)',f->subsop(0=binomial_,f)));
          val_num:=eval(subsindets(subs({n=args[5],k[1]=i,k[2]=j},numer(t)),'specfunc(anything,binomial)',f->subsop(0=binomial_,f)));
        catch: 
          val:=evalf(subsindets(subs({n=args[5],k[1]=i,k[2]=j},t),'specfunc(anything,binomial)',f->subsop(0=binomial_,f)));
          val_num:=evalf(subsindets(subs({n=args[5],k[1]=i,k[2]=j},numer(t)),'specfunc(anything,binomial)',f->subsop(0=binomial_,f)));
        end try;
        if hastype(val,undefined) or hastype(val_num,infinity)
          then undef_pts:=[op(undef_pts),[i,j]];
          else if not(hastype(val,name)) 
                 then if abs(val)>epsilon
                        then pts:=[op(pts),[i,j]]  
                      end if;
                 else pts:=[op(pts),[i,j]]
               end if;
        end if;    
      end do;
    end do;
    return(plots[display](seq(plottools[point](pts[i],symbol=circle,symbolsize=15,color=black),i=1..nops(pts)),
                          seq(plottools[point](undef_pts[i],symbol=circle,symbolsize=15,color=grey),i=1..nops(undef_pts)),
                          axes=normal,scaling=constrained,view=[-bd..bd,-bd..bd]));
  end if;

end proc:


########################################################### 
#                                                         #
#  'drawSmax' zeichnet die Punkte der P-maximalen         #
#  Strukturmenge eines PHT's in ein zwei- bzw. drei-      #
#  dimensionales Koordinatensystem.                       #
#                                                         #
###########################################################

drawSmax:=proc(t,k,n,IJ)
local i,j,Smax_pts,Sstd_pts,extend,extend_,equat,lines,Smax_,Sstd_,num_IJ,Scov_,bdpts,bdpts_,
      plotlist,planes;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  if 4<nargs then extend:=abs(args[5]) else extend:=3 end if;
  Sstd_pts:=Sstd(IJ);
  Smax_pts:=m_Smax(t,k,n,IJ,extend) minus Sstd_pts;
  bdpts:=boundpts_Sstd(t,k,n,IJ)[2];
  num_IJ:=nops(IJ);
  Scov_:=Scov(t,k,n,IJ,0);
  equat:=convert(inequat(t,k,n,IJ,x),equality);  

  if num_IJ=2 then
    lines:=plots[implicitplot](equat,
                x[1]=Scov_[1,1]-extend..Scov_[1,2]+extend,
                x[2]=Scov_[2,1]-extend..Scov_[2,2]+extend,color=black,linestyle=DOT);
    if 4<nargs then if args[5]<0 then lines:=NULL end if end if;

    Smax_:=seq(plottools[point](Smax_pts[i],
                                    symbol=circle,color=grey,symbolsize=15),
              i=1..nops(Smax_pts));
    Sstd_:=seq(seq(plottools[point](Sstd_pts[i],
                                    symbol=circle,color=black,symbolsize=15-j),
              i=1..nops(Sstd_pts)),j=0..14);
    bdpts_:=seq(plottools[point](bdpts[i],
                                    symbol=circle,color=black,symbolsize=22),
              i=1..nops(bdpts));

    return(plots[display](Sstd_,Smax_,bdpts_,lines,
                 axes=normal,scaling=constrained,
                 view=[Scov_[1,1]-extend..Scov_[1,2]+extend,
                       Scov_[2,1]-extend..Scov_[2,2]+extend]));
  end if;   
  
  if num_IJ=3 then
    for j from 1 to 3 do plotlist[j]:=[] end do;
    for i from 1 to nops(equat) do
      for j from 1 to 3 do
        if has(equat[i],x[j]) 
          then plotlist[j]:=[op(plotlist[j]),solve(equat[i],x[j])];
               break;
        end if;
      end do;
    end do;
    if Scov_[1][1]=0 and Scov_[1][2]=0 
      then extend_:=extend;
      else extend_:=0;
    end if;

    planes[1]:=plot3d({seq([plotlist[1][j],x[2],x[3]],j=1..nops(plotlist[1]))},
               x[2]=Scov_[2][1]-extend_..Scov_[2][2]+extend_,
               x[3]=Scov_[3][1]-extend_..Scov_[3][2]+extend_,style=wireframe,color=grey,grid=[2,2]):
    planes[2]:=plot3d({seq([x[1],plotlist[2][j],x[3]],j=1..nops(plotlist[2]))},
               x[1]=Scov_[1][1]-extend_..Scov_[1][2]+extend_,
               x[3]=Scov_[3][1]-extend_..Scov_[3][2]+extend_,style=wireframe,color=grey,grid=[2,2]):
    planes[3]:=plot3d({seq([x[1],x[2],plotlist[3][j]],j=1..nops(plotlist[3]))},
               x[1]=Scov_[1][1]-extend_..Scov_[1][2]+extend_,
               x[2]=Scov_[2][1]-extend_..Scov_[2][2]+extend_,style=wireframe,color=grey,grid=[2,2]):
    if 4<nargs then if args[5]<0 then planes[1]:=NULL;planes[2]:=NULL;planes[3]:=NULL; end if end if;

    Smax_:=seq(plottools[point](Smax_pts[i],
                                    symbol=circle,color=grey,symbolsize=10),
              i=1..nops(Smax_pts));
    Sstd_:=seq(seq(plottools[point](Sstd_pts[i],
                                    symbol=circle,color=black,symbolsize=10-j),
              i=1..nops(Sstd_pts)),j=0..9);
    bdpts_:=seq(plottools[point](bdpts[i],
                                    symbol=box,color=black,symbolsize=15),
              i=1..nops(bdpts));

    return(plots[display3d](planes[1],planes[2],planes[3],Sstd_,Smax_,bdpts_,
                 axes=normal,scaling=constrained,
                 view=[Scov_[1,1]-extend..Scov_[1,2]+extend,
                       Scov_[2,1]-extend..Scov_[2,2]+extend,
                       Scov_[3,1]-extend..Scov_[3,2]+extend]));
  end if;      
end proc:


########################################################### 
#                                                         #
#  'aniSmax' veranschaulicht die Vergrößerungs-           #
#  Strategien der Standard-Strukturmengen.                #
#                                                         #
###########################################################

aniSmax:=proc(t,k,n,K)
local S,i,j,m,min_IJ,r;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  r:=nops(k);
  S:=selectSmax(t,k,n,K);
  for i from 1 to nops(S) do
    min_IJ[i]:=[seq(infinity,m=1..r+1)];
    for j from 1 to nops(S[i][2]) do
      for m from 1 to r+1 do 
        if S[i][2][j][m]<min_IJ[i][m] then min_IJ[i][m]:=S[i][2][j][m] end if;
      end do;
    end do;
  end do;
 
  S:=seq(POINTS(op(map(x->x-min_IJ[i],S[i][2]))),i=1..nops(S));
  return(plots[display](S,symbol=circle,symbolsize=15,axes=normal,
         color=black,scaling=constrained,insequence=true));
end proc:


########################################################### 
#                                                         #
#  'm_selectSmax' sortiert die P-maximalen Struktur-      #
#  Mengen aufsteigend nach der Grösse.                    #
#  Als obere Schranke dient K.                            #                                    
#                                                         #
###########################################################

m_selectSmax:=proc(t,k,n,K)
local i,r,x,y,IJ,Smax_,Smax_sorted,num_IJ;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  r:=nops(k);
  IJ:=Sstd([seq(K,i=1..r+1)]) minus {[seq(0,i=1..r+1)]};
  num_IJ:=nops(IJ);
  Smax_:=[seq([IJ[i],m_Smax(t,k,n,IJ[i])],i=1..num_IJ)];
  Smax_:=sort(Smax_,(x,y)->evalb(nops(x[2])<=nops(y[2])));
  
  Smax_sorted:=[Smax_[1]];
  for i from 2 to num_IJ do
    if evalb(nops(Smax_[i,2])<>nops(Smax_[i-1,2]))
      then Smax_sorted:=[op(Smax_sorted),Smax_[i]]
    end if; 
  end do;
      
  return(Smax_sorted);
end proc:


########################################################### 
#                                                         #
#  'replacevar' ersetzt die auftretenden Variablen in     #
#  einer Rekursionsgleichung durch möglichst einfache     #
#  Zahlen (0 oder 1).                                     #
#                                                         #
###########################################################

replacevar:=proc(recurrence,F,n,a)
local i,j,var,num_var,subs_S,rec_wvar,rules,recurrences;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  var:=select(has,indets(recurrence),a);
  num_var:=nops(var);
  recurrences:={};
  
  if num_var>0 
    then
      for i from 1 to num_var do
        rec_wvar:=subs(var[i]=1,recurrence);
        rec_wvar:=subs({seq(var[j]=0,j=1..num_var)},rec_wvar);
        if rec_wvar<>0 
          then rec_wvar:=numer(rec_wvar);
               rec_wvar:=contentrec_(rec_wvar,F,n,a);
               rec_wvar:=collect(rec_wvar,F,factor);
               recurrences:=recurrences union {rec_wvar=0};
        end if; 
      end do;
  end if;
  
  return(recurrences);
end proc:


########################################################### 
#                                                         #
#  'replacevar_sum' ersetzt die auftretenden Variablen    # 
#  in einer Rekursionsgleichung so, daß möglichst viele   #
#  Summanden wegfallen.                                   #
#                                                         #
###########################################################

replacevar_sum:=proc(recurrence,S,n,a)
local i,rec,max_r,min_r,eq,temp_rec;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  rec:=recurrence;
  min_r:=minrec(rec)[1];
  max_r:=maxrec(rec)[1];
  
  for i from min_r to max_r-2 do
    try eq:=coeff(collect(lhs(rec)-rhs(rec),S),S(n+i));
        temp_rec:=subs(solve(eq,select(has,indets(rec),a)),rec);
        temp_rec:=numer(lhs(temp_rec)-rhs(temp_rec))=0;
        if lhs(temp_rec)<>0 or rhs(temp_rec)<>0 
          then rec:=temp_rec;
        end if; 
    catch:
    end try;
  end do;

  rec:=contentrec_(lhs(rec),S,n,a);
  return(collect(rec,S,factor)=0);      
end proc:


########################################################### 
#                                                         #
#  'contentrec' bestimmt den Inhalt einer Rekursion       #
#  und teilt die Rekursionsgleichung durch diesen.        #
#                                                         #
###########################################################

content_r:=proc(recurrence)
local i,content,rec,F_list,F,recurrence_;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  if type(recurrence,equation) 
    then recurrence_:=lhs(recurrence)-rhs(recurrence);
    else recurrence_:=recurrence;
  end if;
  F_list:=select(type,indets(recurrence_),function);
  F:=op(0,F_list[1]);

  rec:=factor(recurrence_);
  content:=1;
  if type(rec,`*`) 
    then 
      for i from 1 to nops(rec) do
        if type(op(i,rec),freeof(F))
          then content:=content*op(i,rec)
        end if;
      end do;
  end if;

  return([content,collect(normal(rec/content),F,factor)=0]);  
end proc:

contentrec:=proc(recurrence)
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  if type(recurrence,list) or type(recurrence,set)
    then return(map(x->content_r(x),[op(recurrence)]));
    else return(content_r(recurrence));
  end if;
end proc:


########################################################### 
#                                                         #
#  'contentrec_' bestimmt den Inhalt einer Rekursion      #
#  und teilt die Rekursionsgleichung durch Faktoren von   #
#  diesen, dessen Nullstellen nichtnegativ sind.          #
#                                                         #
###########################################################

contentrec_:=proc(recurrence,F,n,a)
local i,content,rec,order;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  rec:=recurrence=0;
  order:=orderrec(rec)-maxrec(rec)[1];
  rec:=factor(lhs(rec));
  content:=1;
  if type(rec,`*`) 
    then 
      for i from 1 to nops(rec) do
        if type(op(i,rec),freeof(F))
          then if has(op(i,rec),a) 
                 then content:=content*op(i,rec);
                 else if nops(select(type,[solve(subs(n=n+order,op(i,rec)),n)],nonnegint))=0
                        then content:=content*op(i,rec);
                      end if;
               end if;
        end if;
      end do;
  end if;

  return(normal(rec/content));  
end proc:


########################################################### 
#                                                         #
#  'rat_linear' gibt die Zahlen [m,l1,...,lr] bei Eingabe # 
#  eines (m,l1,...lr)-fachen hypergeometr. Terms aus.     #
#                                                         #
###########################################################

rat_linear:=proc(t,k,n)
local term,i,m,p,q,pp,qq,r,denom_list;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  r:=nops(k);
  term:=convertfacex(op(convertgam(select(has,simpcomb(t),GAMMA))),k,n);
  pp:=nops(term[1]);
  qq:=nops(term[4]);
  for i from 1 to r+1 do 
    denom_list[i]:=[];
  end do;
  for p from 1 to pp do
    denom_list[1]:=[op(denom_list[1]),denom(term[1][p])];
    for i from 1 to r do 
      denom_list[i+1]:=[op(denom_list[i+1]),denom(term[2][i][p])];
    end do;
  end do;
  for q from 1 to qq do
    denom_list[1]:=[op(denom_list[1]),denom(term[4][q])];
    for i from 1 to r do 
      denom_list[i+1]:=[op(denom_list[i+1]),denom(term[5][i][q])];
    end do;
  end do; 
  
  return([seq(lcm(op(denom_list[m])),m=1..r+1)]);
end proc:


########################################################### 
#                                                         #
#  'updateF' wandelt einen (m,l1,...lr)-fachen hyper-     #
#  geometr. Terms in einen hypergeometr. Term um.         #
#                                                         #
###########################################################

updateF:=proc(t,k,n)
local rat,i,var,r;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  r:=nops(k);
  var:=[n,op(k)];
  rat:=rat_linear(t,k,n);
  return(subs(seq(var[i]=rat[i]*var[i],i=1..r+1),t));
end proc:


########################################################### 
#                                                         #
#  'updateF_back' wandelt eine Rekursion in eine          # 
#  Rekursion bzgl. des Ausgangsterms um.                  #
#                                                         #
###########################################################

updateF_back:=proc(recurrence,t,k,n,F)
local i,m,new_rec,rat_,tt,pol,var,r;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  r:=nops(k);
  var:=[n,op(k)];

  rat_:=rat_linear(t,k,n);
  new_rec:=0;
  for i from 1 to nops(recurrence) do
    tt:=select(has,op(i,recurrence),F);
    pol:=subs({seq(var[m]=var[m]/rat_[m],m=1..r+1)},normal(op(i,recurrence)/tt));
    new_rec:=new_rec+pol*subs({seq(op(m,tt)=var[m]+expand(rat_[m]*(op(m,tt)-var[m])),m=1..r+1)},tt);
  end do;  
    
  return(new_rec);
end proc:


########################################################### 
#                                                         #
#  'krec' gibt eine k-freie Rekursionsgleichung aus,      #
#  sofern eine existiert.                                 #
#                                                         #
###########################################################

krec:=proc(t_,k,n,IJ)
local t,tt,i,j,m,num_IJ,num_Smax,F,a,Smax_,subs_S,var,variables,eq,IJ_,
      sol,recurrence,print,rat_term,arg;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  if 4<nargs then F:=args[5] end if;
  if 5<nargs then a:=args[6] end if;

  tt:=updateF(t_,k,n);

  if _EnvFilter=1
    then t:=filter(tt,k,n); 
    else t:=tt;
  end if;  
  
  print:=1;
  var:=[n,op(k)]; 
  if type(IJ[1],`list`)
    then num_IJ:=nops(IJ[1]);
         IJ_:=IJ[1];
         if nops(IJ)>1 
           then
             if not(type(IJ[2],`list`))
               then Smax_:=IJ[2];
               else Smax_:=IJ;
                    print:=0;
             end if;
           else
             Smax_:=IJ;
             print:=0;
         end if;
    else Smax_:=Smax(t,k,n,IJ);
         IJ_:=IJ;
         num_IJ:=nops(IJ);
  end if;
  num_Smax:=nops(Smax_);

  for i from 1 to num_Smax do
    rat_term[i]:=simpcomb(subs({seq(var[m]=var[m]-Smax_[i,m],m=1..num_IJ)},t)/t);
  end do;
  eq:=sum('a[op(Smax_[i])]*rat_term[i]','i'=1..num_Smax);
  variables:=select(has,indets(eq),a); 
  eq:={coeffs(expand(numer(eq)),k)};
  if _EnvP=1 then  
    if print=1 
      then printf("structureset: %a \t number of equations: %d \t number of variables: %d \n",
                  IJ_,nops(eq),nops(variables));
      else printf("number of equations: %d \t number of variables: %d \n",
                  nops(eq),nops(variables));
    end if;
  end if;
  
  sol:=_EnvSolve(eq,variables);
  if {seq(op(2,sol[i]),i=1..nops(sol))}={0}
    then ERROR(`No k-free recurrence equation exists`)
  end if;
  
  for i from 1 to num_Smax do
    arg[i]:=seq(var[m]-Smax_[i,m],m=1..num_IJ); 
  end do;

  recurrence:=sum('a[op(Smax_[i])]*F(arg[i])','i'=1..num_Smax);
  if _EnvFilter=1 then recurrence:=filter_back(recurrence,tt,k,n,F) end if;
  recurrence:=subs(sol,recurrence);
  recurrence:=updateF_back(recurrence,t_,k,n,F); 
  if nargs<=5 
    then recurrence:=replacevar(recurrence,F,n,a); 
    else recurrence:=numer(normal(recurrence));
         recurrence:=contentrec_(recurrence,F,n,a);
         recurrence:=collect(recurrence,F,factor);
         recurrence:={recurrence=0};
  end if;
   
  return(recurrence);
end proc:


########################################################### 
#                                                         #
#  'findkfreerec' versucht eine k-freie Rekursion         #
#  zu finden. Es werden keine P-maximalen Strukturmengen  #
#  verwendet. Die Strukturmengen werden dabei solange     #
#  vergrössert, bis eine Rekursionsgleichung gefunden     #
#  wird.                                                  #
#                                                         #
###########################################################

findkfreerec:=proc(t,k,n)
local i,j,r,IJ,recurrence,F,a,S;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  F:=NULL; a:=NULL; S:=NULL;
  if 4<nargs then F:=args[5] end if;
  if 5<nargs then a:=args[6] end if;   
  if 3<nargs 
    then if not(type(args[4],set)) and not(type(args[4],list))
            then F:=args[4];
                 if 4<nargs then a:=args[5] end if;
            else S:=args[4];
         end if;
  end if;

  if type(S,set) then return(krec(t,k,n,S,F,a)) end if;
  if type(S,list) then return(krec(t,k,n,Sstd(S),F,a)) end if; 

  r:=nops(k);
  recurrence:={};
  IJ:=[1,seq(0,j=1..r)];
  i:=1; 
  while recurrence={} do
    try recurrence:=krec(t,k,n,[IJ,Sstd(IJ)],F,a);
    catch "No k-free recurrence equation exists":
    end try;
    if i<=r 
      then i:=i+1;
      else i:=1;
    end if;
    IJ[i]:=IJ[i]+1;
  end do; 
  
  return(recurrence);
end proc:


########################################################### 
#                                                         #
#  'maxorder' gibt eine obere Schranke für die Ordnung    #
#  der k-freien Rekursion an.                             #
#                                                         #
###########################################################

maxorder:=proc(t,k,n)
local i,p,q,r,max_ord,propterm,fac_diff;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  propterm:=properterm(t,k,n);
  fac_diff:=facdiff(t,k,n);
  r:=nops(k);
  max_ord:=[seq(0,i=1..r)];
  for i from 1 to r do 
    max_ord[i]:=max_ord[i]+max(-fac_diff[i+1],0);
    for p from 1 to nops(propterm[2,2,1]) do
      max_ord[i]:=max_ord[i]+max(propterm[2][2,i,p],0);
    end do;
    for q from 1 to nops(propterm[2,5,1]) do
      max_ord[i]:=max_ord[i]+max(-propterm[2][5,i,q],0);
    end do;
  end do;
   
  return(ceil((1/(r!))*(sum('max_ord[i]','i'=1..r))^r));
end proc:


########################################################### 
#                                                         #
#  'orderrec' gibt die Ordnungen von Rekursions-          #
#  Gleichungen in der ersten Variablen an (optional:      #
#  in einer bestimmten Variablen).                        #
#                                                         #
###########################################################

order_r:=proc(recurrence)
local i,j,nvar,F_list,v_list;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  F_list:=[op(select(type,indets(lhs(recurrence)-rhs(recurrence)),function))];
  if nargs>1 
    then 
      i:=args[2];
      for j from 1 to nops(F_list[1]) do
        if op(indets(op(j,F_list[1])))=i
          then nvar:=j;
        end if;
      end do;
      v_list:=map(x->(op(nvar,x)-i),F_list);
    else 
      v_list:=map(x->(op(1,x)-op(indets(op(1,F_list[1])))),F_list);
  end if; 
  return(max(op(v_list))-min(op(v_list)));
end proc:

orderrec:=proc(recurrence)
local var;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  if nargs>1 then var:=args[2] else var:=NULL end if;
  if type(recurrence,list) or type(recurrence,set)
    then return(map(x->order_r(x,var),[op(recurrence)]));
    else return(order_r(recurrence,var));
  end if;
end proc:


########################################################### 
#                                                         #
#  'ordercer' gibt die Ordnungen von Zertifikats-         #
#  Rekursionen in der ersten Variablen an (optional:      #
#  in einer bestimmten Variablen).                        #
#                                                         #
###########################################################

ordercer:=proc(certificate)
local var;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  if nargs>1 then var:=args[2] else var:=NULL end if;
  if type(certificate,list)=true or type(certificate,set)=true
    then return(orderrec(map(x->torec(x),certificate),var));
    else return(orderrec(torec(certificate),var));
  end if;
end proc:


########################################################### 
#                                                         #
#  'checkrec' überprüft, ob t die eingegebene             #
#  Rekursionsgleichung erfüllt.                           #
#                                                         #
###########################################################

check_r:=proc(recurrence,t)
local i,j,var,l_rec,F_list,fulfilled,F,rules,subs_S,rec1,rec2,err;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  l_rec:=lhs(recurrence)-rhs(recurrence);

  if type(t,equation)
    then if orderrec(recurrence)>orderrec(t)
           then rec1:=t;
                rec2:=recurrence;
           else rec1:=recurrence;
                rec2:=t;
         end if;
         err:=0;
         while orderrec(rec2)>orderrec(rec1) do
           rec2:=addrec(rec1,rec2);
           err:=err+1;
           if err>200 then ERROR(`Recurrences must contain the same variables`) end if;
         end do;
         rec2:=addrec(rec1,rec2);
         
         if lhs(rec2)=0 and rhs(rec2)=0
           then return(true) 
           else return(false);
         end if;  
  end if;  

  if type(l_rec,`+`)
    then 
      F_list:=[op(select(type,indets(l_rec),function))];
      var:=[];
      for i from 1 to nops(F_list[1]) do
        var:=[op(var),op(indets(op(i,F_list[1])))];
      end do; 

      rules:={};
      for i from 1 to nops(F_list) do
        subs_S:={};
        for j from 1 to nops(F_list[i]) do
          subs_S:=subs_S union {var[j]=op(j,F_list[i])};
        end do;
        rules:=rules union {F_list[i]=subs(subs_S,t)};
      end do;
        
      fulfilled:=simpcomb(subs(rules,l_rec))=0;
      if fulfilled 
        then return(true) 
        else return(false);
      end if;   
  end if;

end proc:

checkrec:=proc(recurrence,t)
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  if type(recurrence,list) or type(recurrence,set)
    then return(map(x->check_r(x,t),[op(recurrence)]));
    else return(check_r(recurrence,t));
  end if;
end proc:


########################################################### 
#                                                         #
#  'checkcer' überprüft, ob t die eingegebene             #
#  Rekursionsgleichung erfüllt.                           #
#                                                         #
###########################################################

check_c:=proc(recurrence,t)
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  return(checkrec(torec(recurrence),t));
end proc:

checkcer:=proc(recurrence,t)
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  if type(recurrence,list) or type(recurrence,set)
    then return(map(x->check_c(x,t),[op(recurrence)]));
    else return(check_c(recurrence,t));
  end if;
end proc:


########################################################### 
#                                                         #
#  'prinpart' gibt den Hauptteil eines Zertifikats aus.   #
#                                                         #
###########################################################

prinpart:=proc(certificate)
local F_list,delta,l_cer,i,k;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  l_cer:=lhs(certificate)-rhs(certificate);
  F_list:=select(type,indets(l_cer),function);
  delta:={};
  for i from 1 to nops(F_list) do
    if type(op(0,F_list[i]),indexed)
      then delta:=delta union {F_list[i]}
    end if;
  end do; 
  return(l_cer-sum(delta[k],k=1..nops(delta)));
end proc:


########################################################### 
#                                                         #
#  'minrec' bestimmt die Minima der Indizes der           #
#  zulässigen hypergeometrischen Terme.                   # 
#                                                         #
###########################################################

minrec:=proc(recurrence)
local i,j,min,r,l_rec,var,F_list;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  l_rec:=lhs(recurrence)-rhs(recurrence);

  F_list:=select(type,indets(prinpart(l_rec=0)),function);
  
  var:=[];
  for i from 1 to nops(F_list[1]) do
    var:=[op(var),op(indets(op(i,F_list[1])))];
  end do; 
  r:=nops(var)-1; 
           
  for j from 1 to r+1 do 
    min[j]:=infinity
  end do;
     
  for i from 1 to nops(F_list) do
    for j from 1 to nops(op(i,F_list)) do
      if evalb(op(j,op(i,F_list))-var[j]<min[j]) 
        then min[j]:=op(j,op(i,F_list))-var[j]
      end if;
    end do;
  end do;   
             
  return([seq(min[i],i=1..r+1)]);
end proc:


########################################################### 
#                                                         #
#  'maxrec' bestimmt die Maxima der Indizes der           #
#  zulässigen hypergeometrischen Terme.                   # 
#                                                         #
###########################################################

maxrec:=proc(recurrence)
local i,j,max,r,l_rec,var,F_list;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  l_rec:=lhs(recurrence)-rhs(recurrence);

  F_list:=select(type,indets(prinpart(l_rec=0)),function);
  
  var:=[];
  for i from 1 to nops(F_list[1]) do
    var:=[op(var),op(indets(op(i,F_list[1])))];
  end do; 
  r:=nops(var)-1; 
           
  for j from 1 to r+1 do 
    max[j]:=-infinity
  end do;
     
  for i from 1 to nops(F_list) do
    for j from 1 to nops(op(i,F_list)) do
      if evalb(op(j,op(i,F_list))-var[j]>max[j]) 
        then max[j]:=op(j,op(i,F_list))-var[j]
      end if;
    end do;
  end do;   
             
  return([seq(max[i],i=1..r+1)]);
end proc:


########################################################### 
#                                                         #
#  'shiftrec' führt eine Indexverschiebung aus.           #
#                                                         #
###########################################################

shift_r:=proc(recurrence)
local i,j,max,min,r,l_rec,var,F_list,a;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  l_rec:=lhs(recurrence)-rhs(recurrence);
  
  F_list:=select(type,indets(l_rec),function);
  var:=[];
  for i from 1 to nops(F_list[1]) do
    var:=[op(var),op(indets(op(i,F_list[1])))];
  end do; 
  
  r:=nops(var)-1;
 
  if nargs>1
    then min:=minrec(recurrence);
         l_rec:=subs(seq(var[j]=var[j]-min[j],j=1..r+1),l_rec);
    else max:=maxrec(recurrence);
         l_rec:=subs(seq(var[j]=var[j]-max[j],j=1..r+1),l_rec);
  end if;
  l_rec:=contentrec_(l_rec,op(0,F_list[1]),op(1,var),a); 
  l_rec:=collect(l_rec,op(0,F_list[1]),factor);  
              
  return(l_rec=0);
end proc:

shiftrec:=proc(recurrence)
local args_;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  if nargs>1 
    then args_:=args[2];
    else args_:=NULL;
  end if;
  if type(recurrence,list) or type(recurrence,set)
    then return(map(x->shift_r(x,args_),{op(recurrence)}));
    else return(shift_r(recurrence,args_));
  end if;
end proc:


########################################################### 
#                                                         #
#  'shiftcer' führt eine Indexverschiebung aus.           #
#                                                         #
###########################################################

shift_c:=proc(recurrence)
local i,j,min,r,l_rec,var,F_list;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  l_rec:=lhs(recurrence)-rhs(recurrence);
  if nargs>1 
    then min:=maxrec(recurrence);
         min[1]:=minrec(recurrence)[1];
    else min:=minrec(recurrence);
         min[1]:=maxrec(recurrence)[1];
  end if;

  F_list:=select(type,indets(l_rec),function);
  
  var:=[];
  j:=1;

  while j<nops(F_list) do 
    if not(type(op(0,F_list[j]),indexed))
      then 
        for i from 1 to nops(F_list[j]) do
          var:=[op(var),op(indets(op(i,F_list[j])))];
        end do;
        break;
      else
        j:=j+1;
    end if;
  end do;  
  r:=nops(var)-1; 
 
  for i from 1 to nops(F_list) do
    if type(op(0,F_list[i]),indexed)
      then 
        for j from 2 to r+1 do
          if evalb(op(1,op(0,F_list[i]))=var[j])
            then 
              l_rec:=subs(op(0,op(0,F_list[i]))[var[j]]=op(0,op(0,F_list[i]))[var[j]+min[j]],l_rec)
          end if;
        end do;
    end if;
  end do;
       
  l_rec:=subs(seq(var[j]=var[j]-min[j],j=1..r+1),l_rec);
                
  return(l_rec=0);
end proc:

shiftcer:=proc(recurrence)
local args_;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  if nargs>1 
    then args_:=args[2];
    else args_:=NULL;
  end if;
  if type(recurrence,list) or type(recurrence,set)
    then return(map(x->shift_c(x,args_),{op(recurrence)}));
    else return(shift_c(recurrence,args_));
  end if;
end proc:


########################################################### 
#                                                         #
#  'tocer' wandelt eine Rekursionsgleichung               #
#  (nicht notwendig k-frei) in eine Rekursionsgleichung   #
#  mit Delta-Ausdrücken (Zertifikats-Rekursion) um.       #
#                                                         #
###########################################################

to_c:=proc(recurrence)
local i,j,r,kk,l_rec,max,rec_list,F,F_list,cer_recurrence,
      delta,Delta,exp,p_part,p_part_test,coeff_,m;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  l_rec:=lhs(recurrence)-rhs(recurrence);
  
  if not(type(l_rec,`+`)) 
    then ERROR(`No non-trivial certificate recurrence exists`)
  end if;
  
  F_list:=[];
  for i from 1 to nops(l_rec) do
    F_list:=[op(F_list),op(select(type,indets(op(i,l_rec)),function))];
  end do;
  F:=op(0,F_list[1]);
  kk:=[];
  for i from 2 to nops(F_list[1]) do
    kk:=[op(kk),op(indets(op(i,F_list[1])))];
  end do;
  r:=nops(kk);

  rec_list:=[];
  for i from 1 to nops(l_rec) do
    rec_list:=[op(rec_list),[op(i,l_rec)/F_list[i],
                            [op(F_list[i])]]];
  end do;

  max:=-minrec(recurrence);
   
  cer_recurrence:=0;  
  for i from 1 to r do 
    delta:=0;
    for j from 1 to nops(rec_list) do
      exp:=[seq(-(rec_list[j,2,l+1]-kk[l]),l=1..r)];
      delta:=delta+sum(subs({seq(kk[m]=kk[m]-(max[m+1]-exp[m]),m=1..i-1),kk[i]=kk[i]-l},
             rec_list[j,1]*F(op(rec_list[j,2]))),l=0..max[i+1]-exp[i]-1);          
    end do;
    delta:=collect(delta,F,factor);
    if delta<>0 
      then cer_recurrence:=cer_recurrence+Delta[kk[i]](delta);
    end if;
  end do;
  
  p_part:=0; 
  p_part_test:=0; 
  for j from 1 to nops(rec_list) do
      exp:=[seq(-(rec_list[j,2,l+1]-kk[l]),l=1..r)];
      coeff_:=subs({seq(kk[m]=kk[m]-(max[m+1]-exp[m]),m=1..r)},rec_list[j,1]);
      p_part:=p_part+coeff_*F(rec_list[j,2,1],seq(kk[m]-max[m+1],m=1..r));
      p_part_test:=p_part_test+coeff_*F(rec_list[j,2,1]);          
  end do;  
  p_part:=collect(p_part,F,factor);
  p_part_test:=collect(p_part_test,F,factor);
  
  if p_part=0 or not(type(p_part_test,freeof(kk)))
    then ERROR(`No non-trivial certificate recurrence exists`)
  end if;
 
  cer_recurrence:=p_part+cer_recurrence;
                  
  if nargs>1 
    then return(p_part)
    else return(cer_recurrence=0)
  end if;  
end proc:

tocer:=proc(recurrence)
local opt;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  if nargs>1 
    then opt:=1 
    else opt:=NULL
  end if;
  if type(recurrence,list) or type(recurrence,set)
    then return(map(x->to_c(x,opt),{op(recurrence)}));
    else return(to_c(recurrence,opt));
  end if;
end proc:


########################################################### 
#                                                         #
#  'torec' wandelt eine Rekursionsgleichung mit Delta-    #
#  Ausdrücken (Zertifikats-Rekursion) in eine             #
#  Rekursionsgleichung (nicht notwendig k-frei) um.       #
#                                                         #
###########################################################

to_r:=proc(certificate)
local i,m,l_cer,F_list,delta_list,F,recurrence,index;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  l_cer:=lhs(certificate)-rhs(certificate);

  F_list:=select(type,indets(l_cer),function);
  delta_list:={}; 
  for i from 1 to nops(F_list) do
    if type(op(0,F_list[i]),indexed) 
      then delta_list:=delta_list union {F_list[i]}
    end if;
  end do;
  F_list:=F_list minus delta_list;  

  F:=op(0,F_list[1]);
  
  recurrence:=subs({seq(delta_list[m]=0,m=1..nops(delta_list))},l_cer);
  for i from 1 to nops(delta_list) do
    index:=op(1,op(0,delta_list[i]));
    recurrence:=recurrence+op(1,delta_list[i]);
    recurrence:=recurrence-subs(index=index-1,op(1,delta_list[i]));
  end do;
  
  return(collect(recurrence,F,factor)=0);  
end proc:

torec:=proc(recurrence)
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  if type(recurrence,list) or type(recurrence,set)
    then return(map(x->to_r(x),{op(recurrence)}));
    else return(to_r(recurrence));
  end if;
end proc:


########################################################### 
#                                                         #
#  'Jmax' bestimmt die Maxima der Komponenten der Punkte  #
#  einer beliebigen Strukturmenge S, wobei die erste      #
#  Komponente nicht betrachtet wird.                      #
#                                                         #
###########################################################

Jmax:=proc(S)
local i,j,m,r,Jmax_;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  r:=nops(S[1])-1;
  Jmax_:=[seq(-infinity,m=1..r)];
  for i from 1 to nops(S) do
    for j from 1 to r do
      if S[i][j+1]>Jmax_[j] 
        then Jmax_[j]:=S[i][j+1]
      end if;
    end do;
  end do;
  return(Jmax_);
end proc:


########################################################### 
#                                                         #
#  '_rec' versucht eine Rekursionsgleichung zu finden     #
#  (bzgl.der P-maximalen Strukturmenge der Standard-      #
#  Strukturmenge Sstd), bei der der Hauptteil k-frei ist. #
#               ( 1. Verallgemeinerung )                  #
#                                                         #
###########################################################

_rec:=proc(t_,k,n,IJ,M)
local t,tt,i,j,kk,l,m,variables,eq,var,F,a,Smax_,num_Smax,num_IJ,pol,pol_set,Jmax_,
      ansatz,rec_eq,reduced_eq,reduced_variables,recurrence,reduced_sol,sol,IJ_;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  if 5<nargs then F:=args[6] end if;
  if 6<nargs then a:=args[7] end if;  
  
  tt:=updateF(t_,k,n);

  if _EnvFilter=1
    then t:=filter(tt,k,n); 
    else t:=tt;
  end if;   

  kk:=[op(k)];
  var:=[n,op(k)];   
  if type(IJ[1],`list`)
    then Smax_:=IJ[2];
         IJ_:=IJ[1];
         num_IJ:=nops(IJ[1]);
    else Smax_:=Smax(t,k,n,IJ);
         IJ_:=IJ;
         num_IJ:=nops(IJ);
  end if;  
  num_Smax:=nops(Smax_);
  Jmax_:=Jmax(Smax_);
  
  ########################################################### 
  
  for i from 1 to num_Smax do
    ansatz[i]:=a[op(Smax_[i]),seq(l[m],m=1..num_IJ-1)]*
               product(kk[m]^l[m],m=1..num_IJ-1)*
               simpcomb(subs({seq(var[m]=var[m]-Smax_[i,m],m=1..num_IJ)},t)/t);
    recurrence[i]:=a[op(Smax_[i]),seq(l[m],m=1..num_IJ-1)]*
                   product(kk[m]^l[m],m=1..num_IJ-1)*
                   F(op(var-Smax_[i]));
    for j from 1 to num_IJ-1 do
      ansatz[i]:=sum(ansatz[i],l[j]=0..M);
      recurrence[i]:=sum(recurrence[i],l[j]=0..M);
    end do;
  end do;
  
  ansatz:=sum('ansatz[i]','i'=1..num_Smax);
  recurrence:=sum('recurrence[i]','i'=1..num_Smax);
  
  ########################################################### 

  for i from 1 to num_Smax do
    pol[op(Smax_[i])]:=a[op(Smax_[i]),seq(l[m],m=1..num_IJ-1)]*
                       product((kk[m]+Smax_[i,m+1]-Jmax[m])^l[m],m=1..num_IJ-1);
    pol[Smax_[i,1]]:=0;
    for j from 1 to num_IJ-1 do
      pol[op(Smax_[i])]:=sum(pol[op(Smax_[i])],l[j]=0..M);
    end do;
  end do;
  
  for i from 1 to num_Smax do
    pol[Smax_[i,1]]:=expand(pol[Smax_[i,1]]+pol[op(Smax_[i])])
  end do;

  pol_set:=[op({seq(pol[Smax_[i,1]],i=1..num_Smax)})];
  for i from 1 to nops(pol_set) do
    if ldegree(pol_set[i],{op(kk)})=0
      then pol_set[i]:=pol_set[i]-tcoeff(pol_set[i],kk)
    end if;
    reduced_variables[i]:=select(has,indets(pol_set[i]),a); 
    reduced_eq[i]:={coeffs(pol_set[i],k)};
    reduced_sol[i]:=solve(reduced_eq[i],reduced_variables[i]);
  end do;
 
  ########################################################### 

  eq:=subs({seq(op(reduced_sol[i]),i=1..nops(pol_set))},ansatz);
  variables:=select(has,indets(eq),a);
  eq:={coeffs(numer(eq),k)};
  
  if _EnvP=1 then    
  printf("structureset: %a \t number of equations: %d \t number of variables: %d \n",
          IJ_,nops(eq),nops(variables));
  end if;
  
  sol:=_EnvSolve(eq,variables); 
  if {seq(op(2,sol[i]),i=1..nops(sol))}={0}
    then ERROR(`No non-trivial certificate recurrence exists`)
  end if;

  recurrence:=subs({seq(op(reduced_sol[i]),i=1..nops(pol_set))},recurrence);
  if _EnvFilter=1 then recurrence:=filter_back(recurrence,tt,k,n,F) end if;  
  recurrence:=subs(sol,recurrence);
  recurrence:=updateF_back(recurrence,t_,k,n,F); 
  if nargs<=6 
    then recurrence:=replacevar(recurrence,F,n,a) 
    else recurrence:=numer(normal(recurrence));
         recurrence:=contentrec_(recurrence,F,n,a);
         recurrence:=collect(recurrence,F,factor);
         recurrence:={recurrence=0};
  end if;
   
  return(recurrence); 
end proc:


########################################################### 
#                                                         #
#  'findrec' versucht eine Rekursionsgleichung zu         #
#  finden (M=0: k-frei und sonst: nicht k-frei).          #
#  Die Strukturmengen werden dabei solange vergrössert,   #
#  bis eine Rekursionsgleichung gefunden wird.            #
#                                                         #
###########################################################

findrec:=proc(t,k,n,K,M)
local i,j,r,IJ,recurrence,cer,F,a,Smax_,temp_rec;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  cer:=1; F:=NULL; a:=NULL;
  if 5<nargs 
    then if type(args[6],name) 
           then F:=args[6];
                if 6<nargs 
                  then a:=args[7]; 
                  else a:=NULL; 
                end if; 
           else cer:=args[6]; 
                if 6<nargs 
                  then F:=args[7];
                       if 7<nargs 
                         then a:=args[8]; 
                         else a:=NULL;
                       end if; 
                  else F:=NULL;
                       a:=NULL; 
                end if;
         end if;
  end if;
  
  recurrence:={};
  if type(K,`list`)
    then if M=0 
           then recurrence:=krec(t,k,n,K,F,a);
           else recurrence:=_rec(t,k,n,K,M,F,a);
         end if;
         if cer=1 
           then
             temp_rec:={};
             for j from 1 to nops(recurrence) do
               try temp_rec:=temp_rec union {tocer(recurrence[j])};
               catch "No non-trivial certificate recurrence exists":
               end try;
             end do;
             recurrence:=temp_rec;                                  
         end if;
    else if K<>0 
           then Smax_:=selectSmax(updateF(t,k,n),k,n,K);
                i:=1;
                if M=0
                  then while i<=nops(Smax_) do
                         try recurrence:=krec(t,k,n,Smax_[i],F,a);
                             if cer=1 
                               then recurrence:=map(x->tocer(x),recurrence);
                             end if;
                             break;
                         catch "No k-free recurrence equation exists":  
                             i:=i+1;
                         end try;     
                       end do;
                  else while i<=nops(Smax_) do
                         try recurrence:=_rec(t,k,n,Smax_[i],M,F,a);
                             if cer=1 
                               then
                                 temp_rec:={};
                                 for j from 1 to nops(recurrence) do
                                   try temp_rec:=temp_rec union {tocer(recurrence[j])};
                                   catch "No non-trivial certificate recurrence exists":
                                   end try;
                                 end do;
                                 recurrence:=temp_rec;                                  
                             end if;
                             if recurrence={} 
                               then i:=i+1;
                               else break;
                             end if;
                         catch "No non-trivial certificate recurrence exists":  
                             i:=i+1;
                         end try;     
                       end do;
                end if;
           else r:=nops(k);
                IJ:=[1,seq(0,i=1..r)];
                i:=1;
                if M=0
                  then while i<=r+1 do
                         try recurrence:=krec(t,k,n,IJ,F,a);
                             if cer=1 
                               then recurrence:=map(x->tocer(x),recurrence);
                             end if;
                             break;
                         catch "No k-free recurrence equation exists":  
                             if i=r+1 then i:=1 else i:=i+1 end if;
                             IJ[i]:=IJ[i]+1;
                         end try;     
                       end do;
                  else while i<=r+1 do
                         try recurrence:=_rec(t,k,n,IJ,M,F,a);
                             if cer=1 
                               then
                                 temp_rec:={};
                                 for j from 1 to nops(recurrence) do
                                   try temp_rec:=temp_rec union {tocer(recurrence[j])};
                                   catch "No non-trivial certificate recurrence exists":
                                   end try;
                                 end do;
                                 recurrence:=temp_rec;                                  
                             end if;
                             if recurrence={} 
                               then if i=r+1 then i:=1 else i:=i+1 end if;
                                    IJ[i]:=IJ[i]+1;
                               else break;
                             end if;
                         catch "No non-trivial certificate recurrence exists":  
                             if i=r+1 then i:=1 else i:=i+1 end if;
                             IJ[i]:=IJ[i]+1;
                         end try;     
                       end do;
                end if;
         end if;
  end if;
  
  if recurrence={}
    then if M=0 
           then ERROR(`No k-free recurrence equation exists`)
           else ERROR(`No non-trivial certificate recurrence exists`)
         end if;
  end if;  

  return(recurrence);
end proc:


########################################################### 
#                                                         #
#  'sumrec' summiert über alle Summationsvariablen        #
#  der Rekursionsgleichungen in der Menge recurrence.     #
#                                                         #
###########################################################

sumrec:=proc(recurrence)
local i,j,m,S,rules,sum_recurrence,recurrences,p_part,F_list,
      minimal,sum_rec,ord,min_ord,new_rec,a;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  minimal:=true;
  if 1<nargs 
    then if args[2]=1 
           then minimal:=false;
           else S:=args[2] 
         end if;
  end if;
  if not(type(recurrence,set)) 
    then recurrences:={recurrence} 
    else recurrences:=recurrence
  end if; 
  
  sum_recurrence:={};
  for i from 1 to nops(recurrences) do
    try p_part:=tocer(recurrences[i],0);
        F_list:=select(type,indets(p_part),function);
        rules:=[];
        for j from 1 to nops(F_list) do
          rules:=[op(rules),F_list[j]=S(op(1,F_list[j]))];
        end do;
        p_part:=subs(rules,p_part);
        p_part:=contentrec_(p_part,S,op(indets(op(1,F_list[1]))),a);    
        sum_recurrence:=sum_recurrence union {collect(p_part,S,factor)=0};
    catch "No non-trivial certificate recurrence exists": 
    end try;
  end do;

  if minimal and nargs<3 then
    sum_rec:=[];
    ord:=orderrec(sum_recurrence);
    min_ord:=min(op(ord));
    for j from 1 to nops(ord) do
      if ord[j]=min_ord
        then sum_rec:=[op(sum_rec),shiftrec(sum_recurrence[j])];
      end if;
    end do; 
    
    if nops(sum_rec)>1 and min_ord>1 
      then for i from 1 to nops(sum_rec) do
             for j from i+1 to nops(sum_rec) do
               new_rec:=addrec(sum_rec[1],sum_rec[2]);
               if new_rec<>(0=0)
                 then return(new_rec); 
               end if;
             end do;
           end do;              
    end if;

    return(sum_rec[1]);
  end if;
    
  return(sum_recurrence);
end proc:


########################################################### 
#                                                         #
#  'sumcer' summiert über alle Summationsvariablen        #
#  der Zertifikats-Rekursionen in der Menge certificate.  #
#                                                         #
###########################################################

sumcer:=proc(certificate)
local i,j,m,S,rules,subs_S,sum_certificate,certificates,p_part,F_list,l_cer,delta_list,
      minimal,sum_rec,ord,min_ord,new_rec,a;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  minimal:=true;
  if 1<nargs 
    then if args[2]=1 
           then minimal:=false;
           else S:=args[2] 
         end if;
  end if;
  if not(type(certificate,set)) 
    then certificates:={certificate} 
    else certificates:=certificate
  end if;

  sum_certificate:={};
  for i from 1 to nops(certificates) do
    l_cer:=lhs(certificates[i])-rhs(certificates[i]);
    F_list:=select(type,indets(l_cer),function);
    delta_list:={}; 
    for j from 1 to nops(F_list) do
      if type(op(0,F_list[j]),indexed) 
        then delta_list:=delta_list union {F_list[j]}
      end if;
    end do;
    F_list:=F_list minus delta_list;
        
    p_part:=subs({seq(delta_list[m]=0,m=1..nops(delta_list))},l_cer);

    rules:=[];
    for j from 1 to nops(F_list) do
      rules:=[op(rules),F_list[j]=S(op(1,F_list[j]))];
    end do;
    p_part:=subs(rules,p_part);
    p_part:=contentrec_(p_part,S,op(indets(op(1,F_list[1]))),a);           
    sum_certificate:=sum_certificate union {collect(p_part,S,factor)=0};
  end do;
  
  if minimal and nargs<3 then
    sum_rec:=[];
    ord:=orderrec(sum_certificate);
    min_ord:=min(op(ord));
    for j from 1 to nops(ord) do
      if ord[j]=min_ord
        then sum_rec:=[op(sum_rec),shiftrec(sum_certificate[j])];
      end if;
    end do; 
     
    if nops(sum_rec)>1 and min_ord>1 
      then for i from 1 to nops(sum_rec) do
             for j from i+1 to nops(sum_rec) do
               new_rec:=addrec(sum_rec[1],sum_rec[2]);
               if new_rec<>(0=0) 
                 then return(new_rec); 
               end if;
             end do;
           end do;            
    end if;

    return(sum_rec[1]);
  end if;
    
  return(sum_certificate);
end proc:


########################################################### 
#                                                         #
#  'decomposeS' zerlegt die Strukturmenge S in Mengen     #
#  S_0,...,S_r, wobei x ein Element von S_i ist, wenn es  #
#  reduzibel bzgl. S und i ist (bei i=0 irreduzibel       #
#  bzgl. S).                                              #
#                                                         #
###########################################################

decomposeS:=proc(S)
local i,j,k,r,de_S,S_;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  r:=nops(S[1])-1;
  S_:=S;

  for i from 1 to r do
    de_S[i]:={};
    for j from 1 to nops(S_) do
      for k from 1 to nops(S) do
        if S_[j]=subsop(i+1=op(i+1,S[k])-1,S[k])
          then de_S[i]:=de_S[i] union {S_[j]};
          break;
        end if;
      end do;
    end do;
    S_:=S_ minus de_S[i];
  end do;
  
  return([S_,seq(de_S[i],i=1..r)]);
end proc:


########################################################### 
#                                                         #
#  '_cer' versucht eine Rekursionsgleichung zu finden     #
#  (bzgl.der P-maximalen Strukturmenge der Standard-      #
#  Strukturmenge Sstd), bei der der Hauptteil k-frei ist. #
#               ( 2. Verallgemeinerung )                  #
#                                                         #
###########################################################

_cer:=proc(tt,k,n,IJ,M)
local t,i,j,m,p,kk,var,Smax_,de_Smax,num_IJ,IJ_,F,a,ansatz,recurrence,Jmax_,
      ansatz_,recurrence_,reduced_variables,reduced_eq,reduced_sol,pol,
      eq,variables,sol,recurrences,pol_set,rules;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  if 5<nargs then F:=args[6] end if;
  if 6<nargs then a:=args[7] end if;  
  
  if _EnvFilter=1
    then t:=filter(tt,k,n); 
    else t:=tt;
  end if;
  
  kk:=[op(k)];
  var:=[n,op(k)];   
  if type(IJ[1],`list`)
    then Smax_:=IJ[2];
         IJ_:=IJ[1];
         num_IJ:=nops(IJ[1]);
    else Smax_:=Smax(t,k,n,IJ);
         IJ_:=IJ;
         num_IJ:=nops(IJ);
  end if; 
  de_Smax:=decomposeS(Smax_);
  Jmax_:=Jmax(Smax_);
  
  ########################################################### 
  
  for i from 1 to num_IJ do
    for j from 1 to nops(de_Smax[i]) do
      rules[i,j]:={seq(var[m]=var[m]-de_Smax[i,j,m],m=1..num_IJ)};
      if i<>1
        then 
          ansatz_[i,j]:=a[op(de_Smax[i,j]),seq(l[m],m=1..num_IJ-1)]*
                        product(kk[m]^l[m],m=1..num_IJ-1)*
                        simpcomb(subs(rules[i,j],t)/t)-
                        subs(var[i]=var[i]-1,a[op(de_Smax[i,j]),seq(l[m],m=1..num_IJ-1)]*
                        product(kk[m]^l[m],m=1..num_IJ-1))*
                        simpcomb(subs(rules[i,j],subs(var[i]=var[i]-1,t))/t);                   
        else     
          ansatz_[i,j]:=a[op(de_Smax[i,j]),seq(l[m],m=1..num_IJ-1)]*
                        product(kk[m]^l[m],m=1..num_IJ-1)*
                        simpcomb(subs(rules[i,j],t)/t);          
      end if;

      recurrence_[i,j]:=a[op(de_Smax[i,j]),seq(l[m],m=1..num_IJ-1)]*
                        product(kk[m]^l[m],m=1..num_IJ-1);
      for p from 1 to num_IJ-1 do
        ansatz_[i,j]:=sum(ansatz_[i,j],l[p]=0..M);
        recurrence_[i,j]:=sum(recurrence_[i,j],l[p]=0..M);
      end do;

      recurrence_[i,j]:=recurrence_[i,j]*F(op(var-de_Smax[i,j]));
    end do;
    
    ansatz_[i]:=sum(ansatz_[i,m],m=1..nops(de_Smax[i]));
    recurrence_[i]:=sum(recurrence_[i,m],m=1..nops(de_Smax[i]));
  end do;
  
  ansatz:=sum('ansatz_[i]','i'=1..num_IJ);
    
  ########################################################### 
  
  for i from 1 to nops(de_Smax[1]) do
    pol[op(de_Smax[1,i])]:=a[op(de_Smax[1,i]),seq(l[m],m=1..num_IJ-1)]*
                           product((kk[m]+de_Smax[1,i,m+1]-Jmax_[m])^l[m],m=1..num_IJ-1);
    pol[de_Smax[1,i,1]]:=0;
    for j from 1 to num_IJ-1 do
      pol[op(de_Smax[1,i])]:=sum(pol[op(de_Smax[1,i])],l[j]=0..M);
    end do;
  end do;
  
  for i from 1 to nops(de_Smax[1]) do
    pol[de_Smax[1,i,1]]:=expand(pol[de_Smax[1,i,1]]+pol[op(de_Smax[1,i])])
  end do;

  pol_set:=[op({seq(pol[de_Smax[1,i,1]],i=1..nops(de_Smax[1]))})];
  for i from 1 to nops(pol_set) do
    if ldegree(pol_set[i],{op(kk)})=0
      then pol_set[i]:=pol_set[i]-tcoeff(pol_set[i],kk)
    end if;
    reduced_variables[i]:=select(has,indets(pol_set[i]),a); 
    reduced_eq[i]:={coeffs(pol_set[i],k)};
    reduced_sol[i]:=solve(reduced_eq[i],reduced_variables[i]);
  end do;
    
  ########################################################### 
  
  eq:=subs({seq(op(reduced_sol[i]),i=1..nops(pol_set))},ansatz);
  variables:=select(has,indets(eq),a);
  eq:={coeffs(numer(eq),k)};
  
  if _EnvP=1 then  
  printf("structureset: %a \t number of equations: %d \t number of variables: %d \n",
          IJ_,nops(eq),nops(variables));
  end if; 
 
  sol:=_EnvSolve(eq,variables);
   
  if {seq(op(2,sol[i]),i=1..nops(sol))}={0}
    then ERROR(`No non-trivial certificate recurrence exists`)
  end if;
  
  ########################################################### 
 
  recurrences:=subs({seq(op(reduced_sol[i]),i=1..nops(pol_set))},
                     [seq(recurrence_[i],i=1..num_IJ)]);
  recurrences:=subs(sol,recurrences);

  recurrence:=recurrences[1];
  for i from 2 to num_IJ do
    recurrence:=recurrence+recurrences[i]-subs(var[i]=var[i]-1,recurrences[i]);
  end do;

  if _EnvFilter=1 then recurrence:=filter_back(recurrence,tt,k,n,F) end if;

  if nargs<=6 
    then recurrence:=replacevar(recurrence,F,n,a) 
    else recurrence:=numer(normal(recurrence));
         recurrence:=contentrec_(recurrence,F,n,a);
         recurrence:=collect(recurrence,F,factor);
         recurrence:={recurrence=0};
  end if;

  return(recurrence);
end proc:


########################################################### 
#                                                         #
#  'findcer' versucht eine Rekursionsgleichung mit        # 
#  Hilfe der zweiten Verallgemeinerung zu finden.         #
#  Die Strukturmengen werden dabei solange vergrössert,   #
#  bis eine Rekursionsgleichung gefunden wird.            #
#                                                         #
###########################################################

findcer:=proc(t,k,n,K,M)
local i,j,r,IJ,recurrence,cer_,F,a,Smax_,temp_rec;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  cer_:=1; F:=NULL; a:=NULL;
  if 5<nargs 
    then if type(args[6],name) 
           then F:=args[6];
                if 6<nargs 
                  then a:=args[7]; 
                  else a:=NULL; 
                end if; 
           else cer_:=args[6]; 
                if 6<nargs 
                  then F:=args[7];
                       if 7<nargs 
                         then a:=args[8]; 
                         else a:=NULL;
                       end if; 
                  else F:=NULL;
                       a:=NULL; 
                end if;
         end if;
  end if;
  
  recurrence:={};
  if type(K,`list`)
    then recurrence:=_cer(t,k,n,K,M,F,a);
         if cer_=1 
           then
             temp_rec:={};
             for j from 1 to nops(recurrence) do
               try temp_rec:=temp_rec union {tocer(recurrence[j])};
               catch "No non-trivial certificate recurrence exists":
               end try;
             end do;
             recurrence:=temp_rec;                                  
         end if;     
    else if K<>0 
           then Smax_:=selectSmax(updateF(t,k,n),k,n,K);
                i:=1;
                while i<=nops(Smax_) do
                  try recurrence:=_cer(t,k,n,Smax_[i],M,F,a);
                      if cer_=1 
                        then 
                          temp_rec:={};
                          for j from 1 to nops(recurrence) do
                            try temp_rec:=temp_rec union {tocer(recurrence[j])};
                            catch "No non-trivial certificate recurrence exists":
                            end try;
                          end do;
                          recurrence:=temp_rec;                                                                
                      end if;
                      if recurrence={} 
                        then i:=i+1;
                        else break;
                      end if;
                  catch "No non-trivial certificate recurrence exists":  
                         i:=i+1;
                  end try;     
                end do;
           else r:=nops(k);
                IJ:=[1,seq(0,i=1..r)];
                i:=1;
                while i<=r+1 do
                  try recurrence:=_cer(t,k,n,IJ,M,F,a);
                      if cer_=1 
                        then
                          temp_rec:={};
                          for j from 1 to nops(recurrence) do
                            try temp_rec:=temp_rec union {tocer(recurrence[j])};
                            catch "No non-trivial certificate recurrence exists":
                            end try;
                          end do;
                          recurrence:=temp_rec;                                  
                      end if;
                      if recurrence={} 
                        then if i=r+1 then i:=1 else i:=i+1 end if;
                             IJ[i]:=IJ[i]+1;
                        else break;
                      end if;
                  catch "No non-trivial certificate recurrence exists":  
                        if i=r+1 then i:=1 else i:=i+1 end if;
                        IJ[i]:=IJ[i]+1;
                  end try;     
                end do;
         end if;
  end if;
  
  if recurrence={}
    then ERROR(`No non-trivial certificate recurrence exists`)
  end if; 

  return(recurrence);
end proc:


########################################################### 
#                                                         #
#  'solverec' bestimmt die Lösung einer Rekursion, deren  #
#  Strukturmenge genau zwei Elemente besitzt.             #
#                                                         #
###########################################################

solverec:=proc(recurrence)
local i,j,n,S,S_list,s0,rat,cont,num,den,numlist,denlist,lc,shift,sh,order,sol,m,arg_pw;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  order:=orderrec(recurrence);
  S_list:=select(type,indets(recurrence),function);
  if nops(S_list)>2 
    then ERROR(`solverec expects only two summands in recurrence equation`)  
  end if;
  n:=op(indets(op(1,S_list[1])));
  S:=op(0,S_list[1]);
  rat:=normal(solve(shiftrec(recurrence,1),S(n+order))/S(n));
  cont:=contentrec(shiftrec(recurrence,1))[1];
  num:=numer(rat)*cont;
  den:=denom(rat)*cont; 
  lc:=lcoeff(num,n)/lcoeff(den,n);
  numlist:=normal([solve(num,n)]);
  numlist:=[seq(-j,j=numlist)];
  denlist:=normal([solve(den,n)]);
  denlist:=[seq(-j,j=denlist)];

  if order=1 
    then shift:=0;
         if nargs=2 
           then shift:=op(1,lhs(args[2]));
                numlist:=map(x->x+shift,numlist);
                denlist:=map(x->x+shift,denlist);
         end if;
         sh:=-min(1,op(select(type,{op(numlist),op(denlist)},nonposint)))+1; 
         numlist:=map(x->x+sh,numlist);
         denlist:=map(x->x+sh,denlist);
         shift:=shift+sh;
         s0:=S(shift);
           
         if member(1,denlist,'i') 
           then denlist:=subsop(i=NULL,denlist);
           else numlist:=[op(numlist),1];
         end if;

         sol:=eval(subs({n=n-shift,pochhammer=`simpcomb/pochhammer`},s0*hyperterm(numlist,denlist,lc,n)));
         if nargs=2 
           then sol:=subs(args[2],sol);
         end if;
         if sol<>0 then
           sol:=applyrule([(-1)^(k::anything)*pochhammer(n::anything,k::anything)/(k::anything)!=binomial(-n,k),
                           (-1)^(k::anything)*pochhammer(n::anything,k::anything)/pochhammer(1,k::anything)=binomial(-n,k)],
                           sol);
         end if; 
         #sol:=convert(sol,binomial);
         return(sol);
    else for i from 1 to nops(numlist) do
           numlist[i]:=(numlist[i]+l)/order;
         end do;
         for i from 1 to nops(denlist) do
           denlist[i]:=(denlist[i]+l)/order;
         end do;
         numlist:=[op(numlist),1];
                 
         sol:=eval(subs(pochhammer=`simpcomb/pochhammer`,S(l)*order^((nops(numlist)-1-nops(denlist))*((n-l)/order))*hyperterm(numlist,denlist,lc,(n-l)/order)));
         arg_pw:=[];
         for m from 1 to order do
           arg_pw:=[op(arg_pw),irem(n,order)=m-1];
           arg_pw:=[op(arg_pw),subs(l=m-1,sol)]; 
         end do;
         if nargs=2 
           then arg_pw:=subs(args[2],arg_pw); 
         end if;
         arg_pw:=applyrule([(-1)^(k::anything)*pochhammer(n::anything,k::anything)/(k::anything)!=binomial(-n,k),
                            (-1)^(k::anything)*pochhammer(n::anything,k::anything)/pochhammer(1,k::anything)=binomial(-n,k)],
                            arg_pw);
         #arg_pw:=convert(arg_pw,binomial); 
         return(piecewise(op(arg_pw)));
  end if;

end proc:


########################################################### 
#                                                         #
#  'multclosedform' bestimmt zu einem zulässigen Term t   #
#  eine geschlossene Form, falls eine Rekursion           #
#  erster Ordnung für die Summe existiert.                #
#                                                         #
###########################################################

multclosedform:=proc(t,k,sn,init)
local i,K,j,r,M,sumrec,s0,rat,num,den,numlist,denlist,lc,S,n,F,a,hypsol;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  r:=nops(k);
  if type(sn,function) 
    then S:=op(0,sn); 
         n:=op(1,sn);
    else n:=sn;
  end if;
  if type(init,set)
    then S:=init[1];
    else S:=init;
  end if;
  S:=op(0,lhs(S));
 
  K:=1;
  M:=1;
  if 4<nargs 
    then M:=args[5];
  end if; 
  if 5<nargs
    then K:=args[5];
         M:=args[6];
  end if;
             
  sumrec:={};
  while sumrec={} do
    try #sumrec:=sumcer(findrec(t,k,n,K,M),S);
        sumrec:=sumcer(findrec(t,k,n,K,M,F,a),S);
        sumrec:=replacevar_sum(sumrec,S,n,a);
    catch "No non-trivial certificate recurrence exists":
    end try;
    K:=K+1;
  end do;
    
  try return(solverec(sumrec,init))
  catch: 
  end try; 

  hypsol:=petkovsek(sumrec,init);
  if hypsol={} 
    then ERROR(`There exists no closed form`); 
    else return(hypsol);
  end if;    
end proc:


########################################################### 
#                                                         #
#  'multsumrecursion' bestimmt zu einem zulässigen Term t #
#  eine Rekursion minimaler Ordnung für die Summe.        #
#                                                         #
###########################################################

multsumrecursion:=proc(t,k,sn)
local i,K,j,r,M,sumrec,s0,rat,num,den,numlist,denlist,lc,S,n,init,cl_form,F,a;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  r:=nops(k);
  cl_form:=false;
  if type(sn,function) 
    then S:=op(0,sn); 
         n:=op(1,sn);
    else n:=sn;
  end if;

  K:=1;
  if 3<nargs
    then 
      if type(args[4],integer)
        then M:=args[4];
             if 4<nargs
               then K:=args[4];
                    M:=args[5];
             end if;
        else if type(args[4],set) or type(args[4],equation)
               then cl_form:=true;
                    init:=args[4];
                    if type(init,set)
                      then S:=init[1];
                      else S:=init;
                    end if;
                    S:=op(0,lhs(S));           
             end if;
             if 4<nargs 
               then M:=args[5];
               else M:=1;
             end if; 
             if 5<nargs
               then K:=args[5];
                    M:=args[6];
             end if;
      end if;
    else M:=1;
  end if;
             
  sumrec:={};
  while sumrec={} do
    try #sumrec:=sumcer(findrec(t,k,n,K,M),S);
        sumrec:=sumcer(findrec(t,k,n,K,M,F,a),S);
        sumrec:=replacevar_sum(sumrec,S,n,a);
    catch "No non-trivial certificate recurrence exists":
    end try;
    K:=K+1;
  end do;
    
  if cl_form
    then try return(solverec(sumrec,init))
         catch:
         end try; 
  end if;
 
  return(sumrec);
end proc:


########################################################### 
#                                                         #
#  'addrec' addiert zwei Rekursionsgleichungen. Dabei     #
#  wird ein Summand eliminiert.                           #
#                                                         #
###########################################################

addrec:=proc(rec1,rec2)
local i,r1,r2,cf1,cf2,F_list,var,num_IJ,F,new_rec,a;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  r1:=shiftrec(rec1);
  r2:=shiftrec(rec2);
  F_list:=select(type,indets(r1),function);
  num_IJ:=nops(F_list[1]);
  var:=seq(op(indets(op(i,F_list[1]))),i=1..num_IJ);
  F:=op(0,F_list[1]);
  cf1:=coeff(lhs(r1),F(var));  
  cf2:=coeff(lhs(r2),F(var));
  new_rec:=lhs(expand(cf2*r1-cf1*r2));
  if new_rec=0 
    then return(0=0);
    else return(shiftrec(contentrec_(new_rec,F,var[1],a)=0));
  end if;
end proc:


########################################################### 
#                                                         #
#  'addcer' addiert zwei Zertifikats-Rekursionsgleich-    #
#  ungen. Dabei wird ein Summand eliminiert, falls man    #
#  kein zusätzliches Argument eingibt.                    #
#                                                         #
###########################################################

addcer:=proc(cer1,cer2)
local new_cer,i,j,F_list,delta,num_IJ,var,c1,c2,F,cf1,cf2,delta_c1,delta_c2,all_j,all_i;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  F_list:=select(type,indets(prinpart(shiftcer(cer1))),function);
  num_IJ:=nops(F_list[1]);
  var:=seq(op(indets(op(i,F_list[1]))),i=1..num_IJ);
  F:=op(0,F_list[1]);
  if nargs<=2 
    then c1:=shiftcer(cer1);
         c2:=shiftcer(cer2);
         cf1:=coeff(prinpart(c1),F(var));  
         cf2:=coeff(prinpart(c2),F(var));
    else c1:=cer1;
         c2:=cer2;
         if type(args[3],negative) 
           then cf1:=1;
           else cf1:=-1;
         end if;  
         cf2:=1;
  end if;
  new_cer:=collect(expand(cf2*prinpart(c1)-cf1*prinpart(c2)),F,factor);
  delta_c1:=indets(lhs(c1)-prinpart(c1),'typefunc(anything, indexed)');
  delta_c2:=indets(lhs(c2)-prinpart(c2),'typefunc(anything, indexed)');

  all_i:=[seq(false,k=1..nops(delta_c2))];
  for i from 1 to nops(delta_c1) do
    all_j:=false;
    for j from 1 to nops(delta_c2) do
      if op(1,op(0,delta_c1[i]))=op(1,op(0,delta_c2[j])) 
        then new_cer:=new_cer+Delta[op(1,op(0,delta_c1[i]))](collect(
             expand(cf2*op(delta_c1[i])-cf1*op(delta_c2[j])),F,factor));
             all_j:=true;
             all_i[j]:=true;
      end if;
    end do;
    if all_j=false 
      then new_cer:=new_cer+Delta[op(1,op(0,delta_c1[i]))](collect(
           expand(cf2*op(delta_c1[i])),F,factor));
    end if;
  end do;
  for j from 1 to nops(delta_c2) do
    if all_i[j]=false 
      then new_cer:=new_cer+Delta[op(1,op(0,delta_c2[j]))](collect(
           expand(-cf1*op(delta_c2[j])),F,factor));
    end if;
  end do;

  return(new_cer=0);
end proc:


########################################################### 
#                                                         #
#  'filter' filtert alle Faktoren von einem zulässigen    #
#   Term aus, die nicht von den Summationsvariablen k     #
#   abhängen.                                             #
#                                                         #
###########################################################

filter:=proc(t,k,n)
local kfree,tt;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  tt:=simpcomb(t);
  kfree:=select(type,tt,freeof(k));
  
  return(normal(tt/kfree));
end proc:


########################################################### 
#                                                         #
#  'filter_back' macht das Filtern in der eingegebenen    #
#   Rekursion wieder rückgängig.                          #
#                                                         #
###########################################################

filter_back:=proc(recurrence,t,k,n,F)
local i,m,kfree,new_rec;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  kfree:=select(type,simpcomb(t),freeof(k));
  new_rec:=0;
  for i from 1 to nops(recurrence) do
    m:=op(1,select(has,op(i,recurrence),F));
    new_rec:=new_rec+op(i,recurrence)/normal(simpcomb(subs(n=m,kfree)/kfree));
  end do;  
    
  return(new_rec);
end proc:


########################################################### 
#                                                         #
#  'petkovsek' erleichtert die Eingabe und verbessert die #
#   Ausgabe der Prozedur 'rechyper' aus dem Package hsum. #
#                                                         #
###########################################################

petkovsek:=proc(rec)
local i,L,hyper_sol,order,opt_arg,S,n,S_list;
  order:=orderrec(rec);
  S_list:=select(type,indets(rec),function);
  n:=op(indets(op(1,S_list[1])));
  S:=op(0,S_list[1]);

  L:=rechyper(lhs(subs(n=n+order,shiftrec(rec))),S(n));
  if nargs>1 
    then opt_arg:=args[2];
         if not(type(opt_arg,equation)) 
           then return(L);
         end if; 
    else opt_arg:=NULL;
  end if;
  hyper_sol:={};
  for i from 1 to nops(L) do
    hyper_sol:={solverec(denom(L[i])*S(n+1)-numer(L[i])*S(n)=0,opt_arg)} union
               hyper_sol;
  end do;
 
  return(hyper_sol); 
end proc:


########################################################### 
#                                                         #
#  'duration' misst die Zeit einer Prozedur proc.         #
#  Als Argument erwartet die Funktion 'proc'.             #
#                                                         #
###########################################################

duration:=proc(f)
local start,sol;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  start:=time():
  sol:=eval(f);
  return([time()-start,sol]);
end proc:


########################################################### 


if nops([libname])>1 
  then multsum_LIB:=cat(libname[-1],"/multsum.dll") 
  else multsum_LIB:=cat(libname,"/multsum.dll")
end if:


c_properterm:=proc(t,k,n)
local pol_gam,abc_uvw,var;
option remember, `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  var:=[n,op(k)];
  check_hyper(t,k,n);
  pol_gam:=check_proper(t,k,n);
  abc_uvw:=convertgam(pol_gam[2]);
  
  return([nops(abc_uvw[1]),
          nops(abc_uvw[2]),
          ListTools[Flatten](convertfacex(abc_uvw[1],abc_uvw[2],k,n,1))]);  
end proc:

try
c_setfacex:=
define_external(
  'set_fac_ex',
    'n'           :: integer[4],
    'pp'          :: integer[4],
    'qq'          :: integer[4],
    'ab_uv'       :: ARRAY(datatype=integer[4]),
  LIB=multsum_LIB
):


c_facdiff:=
define_external(
  'fac_diff',
  LIB=multsum_LIB
):


c_structurefct:=
define_external(
  'structure_fct',
  LIB=multsum_LIB
):


c_inequat:=
define_external(
  'inequat',
    'IJ'          :: ARRAY(datatype=integer[4]),
  LIB=multsum_LIB
):


c_setminmaxIJ:=
define_external(
  'set_min_max_IJ',
    'min_max_IJ'  :: ARRAY(datatype=integer[4]),
  LIB=multsum_LIB
):


c_getSmax:=
define_external(
  'Smax',
    'Smax_'       :: ARRAY(datatype=integer[4]),
  RETURN          :: boolean[1],  
  LIB=multsum_LIB
):


c_getnumSmax:=
define_external(
  'get_num_Smax',
  RETURN          :: integer[4],
  LIB=multsum_LIB
):


c_getS:=
define_external(
  'select_Smax',
    'max_IJ'      :: integer[4],
    'S'           :: ARRAY(datatype=integer[4]),
  RETURN          :: boolean[1],  
  LIB=multsum_LIB
):


c_getnumS:=
define_external(
  'get_num_S',
  RETURN          :: integer[4],
  LIB=multsum_LIB
):
catch: _EnvC:=0;
end try:

c_Smax:=proc(t,k,n,IJ)
local num_IJ,propterm,IJ_,pp,qq,ab_uv,Smax_,i,j,min_max_IJ,extend,all_Smax,max_num_Smax;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  if 4<nargs then extend:=args[5] else extend:=_EnvSmax end if;
  num_IJ:=nops(IJ);
  max_num_Smax:=100000;

  propterm:=c_properterm(t,k,n);
  pp:=propterm[1];
  qq:=propterm[2];
  ab_uv:=Array(propterm[3],datatype=integer[4],order=C_order);  
  IJ_:=Array(IJ,datatype=integer[4],order=C_order);
  Smax_:=Matrix(max_num_Smax,num_IJ,datatype=integer[4],order=C_order);
  min_max_IJ:=Array(ListTools[Flatten]([seq([-extend,extend],i=1..num_IJ)]),datatype=integer[4],order=C_order);

  c_setfacex(num_IJ,pp,qq,ab_uv); 
  c_setminmaxIJ(min_max_IJ);
  c_facdiff();
  c_structurefct();
  c_inequat(IJ_);
  
  all_Smax:=c_getSmax(Smax_);
  if not(all_Smax)
    then WARNING(`Smax contains probably more elements`) 
  end if; 
 
  return({seq([seq(Smax_[i,j],j=1..num_IJ)],i=1..c_getnumSmax())});
end proc:  


c_selectSmax:=proc(t,k,n,max_IJ)
local num_IJ,propterm,IJ_,pp,qq,ab_uv,S,i,j,m,min_max_IJ,extend,all_S,max_num_S,
      num_Smax,S_,Smax_,Smax_sorted;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  if 4<nargs then extend:=args[5] else extend:=_EnvSmax end if;
  num_IJ:=nops(k)+1;
  max_num_S:=1000000;

  propterm:=c_properterm(t,k,n);
  pp:=propterm[1];
  qq:=propterm[2];
  ab_uv:=Array(propterm[3],datatype=integer[4],order=C_order);  
  S:=Array(1..max_num_S,datatype=integer[4],order=C_order);
  min_max_IJ:=Array(ListTools[Flatten]([seq([-extend,extend],i=1..num_IJ)]),datatype=integer[4],order=C_order);
  
  c_setfacex(num_IJ,pp,qq,ab_uv); 
  c_setminmaxIJ(min_max_IJ);
  c_facdiff();
  c_structurefct();
    
  all_S:=c_getS(max_IJ,S);
  if not(all_S)
    then WARNING(`Some Smax contain probably more elements`) 
  end if; 
  
  S_:=[];
  i:=1;
  while i<c_getnumS() do
    IJ_:=[seq(S[j],j=i..i+num_IJ-1)];
    num_Smax:=S[i+num_IJ];
    Smax_:={};
    for j from i+num_IJ+1 by num_IJ to i+num_IJ+1+num_Smax*num_IJ-1 do
      Smax_:=Smax_ union {[seq(S[m],m=j..j+num_IJ-1)]};
    end do;
    S_:=[op(S_),[IJ_,Smax_]];
    i:=i+num_IJ+1+num_Smax*num_IJ;   
  end do;
 
  S_:=sort(S_,(x,y)->evalb(nops(x[2])<=nops(y[2])));
  Smax_sorted:=[S_[1]];
  for i from 2 to nops(S_) do
    if evalb(nops(S_[i,2])<>nops(S_[i-1,2]))
      then Smax_sorted:=[op(Smax_sorted),S_[i]]
    end if; 
  end do;
   
  return(Smax_sorted);
end proc: 


selectSstd:=proc(t,k,n,max_IJ)
local i,IJ_Smax,extend;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  if 4<nargs then extend:=args[5] else extend:=NULL end if;
  IJ_Smax:=selectSmax(t,k,n,max_IJ,extend);
  
  return([seq(IJ_Smax[i][1],i=1..nops(IJ_Smax))]);
end proc: 


selectSmax:=proc(t,k,n,max_IJ)
local extend;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  if 4<nargs then extend:=args[5] else extend:=NULL end if;
  if _EnvC=1
    then return(c_selectSmax(t,k,n,max_IJ,extend))
    else return(m_selectSmax(t,k,n,max_IJ))
  end if;
end proc: 


Smax:=proc(t,k,n,max_IJ)
local extend;
option `Copyright 2003-2005 Torsten Sprenger, University of Kassel`;
  if 4<nargs then extend:=args[5] else extend:=NULL end if;
  if _EnvC=1
    then return(c_Smax(t,k,n,max_IJ,extend))
    else return(m_Smax(t,k,n,max_IJ,extend))
  end if;
end proc: 


protect(check_hyper,check_proper,convertgam,convertfacex,properterm,checkcer,findkfreerec,
        facdiff,structurefct,boundpts_Sstd,boundpts_S,Sstd,inequat,isfinite,Scov,multclosedform,
        Smax,drawSmax,selectSmax,krec,checkrec,shiftrec,_rec,findcer,torec,maxrec,
        tocer,findrec,sumrec,sumcer,decomposeS,_cer,duration,replacevar,order_r,multsumrecursion,
        c_properterm,c_setfacex,c_structurefct,c_inequat,c_setminmaxIJ,aniSmax,
        c_getSmax,c_Smax,c_getnumSmax,c_getS,c_getnumS,c_facdiff,orderrec,ordercer,
        c_selectSmax,selectSstd,m_selectSmax,m_Smax,minrec,shiftcer,Jmax,multsum_LIB,
        to_c,to_r,shift_c,shift_r,check_c,check_r,solverec,drawsupport,addrec,
        addcer,prinpart,binomial_,maxorder,filter,updateF,rat_linear,petkovsek,replacevar_sum,
        checkhyper,checkproper,boundptsSstd,contentrec_,contentrec,content_r):


print(`Package "Multiple Summation", Maple 8 - Maple 9`);
print(`Copyright 2003-2005, Torsten Sprenger, University of Kassel`);