% 29.11.2003

% This procedure calculates a number of unknowns coefficients
% of the first order ODE.
% m is maximal degree if derivative.
% Solutions tend to infinity as 1/t**p.

 procedure quvar(m,p);
 begin
 scalar numb, k, j;
 numb:=0;
 for k:=0:m-1 do
     << j:=0;
        while p*j <= (p+1)*(m-k) do
        << numb := numb+1;
           j:=j+1;
        >>;
     >>;
 return numb;
 end;


% This procedure constructs the first order ODE.
% a(i,j) is an operator of coefficients.
% m is maximal degree if derivative.
% Solutions tend to infinity as 1/t**p.

 procedure equa(a,m,p,y,dy);
 begin
 scalar equ, k, j;
 equ:=0;
 for k:=0:m do
     << j:=0;
        while p*j <= (p+1)*(m-k) do
        << equ := equ+a(j,k)*y^j*dy^k;
           j:=j+1;
        >>;
     >>;
 return equ;
 end;

% This procedure constructs  the Laurent series for the first order ODE.
% The straight algorithm  is used.
% a(i,j) is an operator of coefficients.
% m is maximal degree if derivative.
% Solutions tend to infinity as 1/t**p.
% c(j) are coefficients of the Laurent series of the function y.

 procedure equalaur(a,m,p,ove,c,tt);
 begin
 scalar equ, equlist, Laurentmax,y,dy;
 Laurentmax:=quvar(m,p)+ove;
 y:=for k:= -p: Laurentmax-p sum c(k)*tt^k;
 dy:=df(y,tt);
 equ:=equa(a,m,p,y,dy);
 equ:=equ*tt^((p+1)*m);
 tt**(Laurentmax+2):=0;
 equ:=equ;
 equlist:={};
 for k:=1:Laurentmax+1 do
    <<  equlist:=append(equlist,{sub(tt=0,equ)});
        equ:=df(equ,tt)/k;
    >>;
 return equlist;
 end;



% This procedure constructs the first order ODE as the following list:
% {{a(0,0),0,0},{a(1,0),1,0},.....,{a(0,m),0,m}}
% a(i,j) is an operator of coefficients.
% m is maximal degree if derivative.
% Solutions tend to infinity as 1/t**p.

 procedure equalist(a,m,p);
 begin
 scalar fequlist, k, j;
 fequlist:={};
 for k:=0:m do
     << j:=0;
        while p*j <= (p+1)*(m-k) do
        << fequlist := append(fequlist,{{a(j,k),j,k}});
           j:=j+1;
        >>;
     >>;
 return fequlist;
 end;



% The following procedures monomlaur(c,mon,j,p), dydegree(c,dydeg,j,p)
% and ydegree(c,ydeg,j,p) construct j-th term of the Laurent series of
% monom coef*y^ydeg*dy^dydeg, where y=for k=-p:INFINITY sum c(k)*t^k,
% dy:=df(y,t) (derivative of y).
% mon={coef,ydeg,dydeg}.


 procedure ydegree(c,n,j,p);
 begin
 scalar stepp;
 if p<1 then stepp:=p else stepp:=1;
 return if n=1 then c(j)
           else
           for k:=-p STEP stepp UNTIL (j+p*n) sum c(k)*ydegree(c,n-1,j-k,p);
 end;


 procedure dydegree(c,n,j,p);
 begin
 scalar stepp;
 if p<1 then stepp:=p else stepp:=1;
 return if n=1 then (j+1)*c(j+1)
           else
               for k:=-p-1 STEP stepp UNTIL (j+(p+1)*n)
                   sum (k+1)*c(k+1)*dydegree(c,n-1,j-k,p);
 end;


 procedure monomlaur(c,mon,j,p);
 begin
 scalar stepp,coef,ydeg,dydeg,k,res;
 coef:=part(mon,1);
 ydeg:=part(mon,2);
 dydeg:=part(mon,3);
 if p<1 then stepp:=p else stepp:=1;
 if ydeg=0
        then << if dydeg=0
             then << if j=0
                       then res:=coef
                       else res:=0
                  >>
             else res:=coef*dydegree(c,dydeg,j,p)
        >>
   else << if dydeg=0
             then res:=coef*ydegree(c,ydeg,j,p)
             else 
	       << if p<1 then stepp:=p else stepp:=1;
	          res:= coef*(for k:=-p*ydeg STEP stepp UNTIL j+(p+1)*dydeg
                          sum ydegree(c,ydeg,k,p)*dydegree(c,dydeg,j-k,p));
	       >>; 		  
        >>;
 return res;
 end;



% This procedure constructs j-th term of the Laurent series of the first order
% ODE.
% fequlist is the first order ODE presented as the list results by procedure
% equalist.
% y=for k=-p:INFINITY sum c(k)*t^k,

 procedure oneequlaur(c,fequlist,j,p);
 begin
 scalar k, equj;
 equj:=0;
 for k:=1:length(fequlist) do
     equj:=equj+monomlaur(c,part(fequlist,k),j,p);
 return equj;
 end;


% This procedure constructs the Laurent series of the first order ODE.
% fequlist is the first order ODE presented as the list results by procedure
% equalist.
% y=for k=-p:Laurentmax-p sum c(k)*t^k,

 procedure equlaurlist(a,m,p,ove,c,pr);
 begin
 scalar stepp, k, laurlist, fequlist, Laurentmax, equk;
 if p<1 then stepp:=p else stepp:=1;
 Laurentmax:=(ove-1)*stepp;
 fequlist:=equalist(a,m,p);
 for k:=-m*(p+1) STEP stepp UNTIL -p-stepp do c(k):=0;
 laurlist:={};
 for k:=-m*(p+1) STEP stepp UNTIL Laurentmax-m*(p+1) do
    << equk:=oneequlaur(c,fequlist,k,p);
       if pr=1 then write "equ(",k,")=",equk;
       laurlist:=append(laurlist,{equk});
    >>;
 return laurlist;
 end;



% This procedure assists to symplify
% the obtaining system of equations.

 procedure simplequ(m,p);
 begin
 scalar simplist, list2, len, pow, k,j,rr;
 len:=1;
 simplist:={{0,0,0}};
 for k:=0:m-1 do
     << j:=0;
        while p*j <= (p+1)*(m-k) do
             << rr:=1;
                list2:={};
                pow:=-(p*j+(p+1)*k);
                while (rr<=len) and (part(simplist,rr,1)<pow) do
                    << list2:=append(list2,{part(simplist,rr)});
                       rr:=rr+1;
                    >>;
                if (rr>len) or (part(simplist,rr,1)>pow) then
                   << list2:=append(list2,{{-(p*j+(p+1)*k),j,k}});
                      while (rr<=len) do
                          << list2:=append(list2,{part(simplist,rr)});
                             rr:=rr+1;
                          >>;
                      simplist:=list2;
                      len:=len+1;
                   >>;
                j:=j+1;
             >>;
     >>;
 return simplist;
 end;

% This procedure checks that term {i,j} belong  to equlist.
% The structure of equlist is equlist:={{something1,n1,m1}{something2,n2,m2}}.

 procedure belonglist(equlist,term);
 begin
 return if equlist={} then 0
           else if part(equlist,1,2)=part(term,1) and
               part(equlist,1,3)=part(term,2) then 1
        else belonglist(rest(equlist),term);
 end;



% This procedure assists to symplify
% the obtaining system of equations.

 procedure simplequ2(equlist,m,p);
 begin
 scalar simplist, list2, len, pow, k,j,rr;
 simplist:={{0,0,0}};
 len:=1;
 for k:=0:m-1 do
     << j:=0;
        while p*j <= (p+1)*(m-k) do
         << if belonglist(equlist,{j,k})=1 then
            << rr:=1;
               list2:={};
               pow:=-(p*j+(p+1)*k);
               while (rr<=len) and (part(simplist,rr,1)<pow) do
                   << list2:=append(list2,{part(simplist,rr)});
                      rr:=rr+1;
                   >>;
               if (rr>len) or (part(simplist,rr,1)>pow) then
                  <<list2:=append(list2,{{-(p*j+(p+1)*k),j,k}});
                    while (rr<=len) do
                         << list2:=append(list2,{part(simplist,rr)});
                            rr:=rr+1;
                         >>;
                    simplist:=list2;
                    len:=len+1;
                  >>;
            >>;
            j:=j+1;
         >>;
     >>;
 return simplist;
 end;



% This procedure gives the list of remainining unknowns a(i,j).

 procedure varlist(a,m,p,simlist);
 begin
 scalar k, j, js, len, vlist, var;
 len:=length(simlist);
 var:={};
 for k:=0:m-1 do
     << j:=0;
        while p*j <= (p+1)*(m-k) do
        << vlist:={a(j,k)};
           for js:=1:len do
              if k=part(simlist,js,3) and j=part(simlist,js,2) then vlist:={};
           var:=append(var,vlist);
           j:=j+1;
        >>;
     >>;
 return var;
 end;



% This procedure gives the list of remainining unknowns a(i,j) in
% the form {...,{a(i,j),i,j},...}.

procedure minusimp(a,simlist);
 begin
 scalar k, j, list3, len, term, ad;
 len:=length(simlist);
 list3:={};
 for j:=1:length(a) do
 << term:=part(a,j);
    ad:=1;
    for k:=1:len do
       if(part(term,2)=part(simlist,k,2)
          and part(term,3)=part(simlist,k,3))
         then ad:=0;
    if ad=1 then list3:=append(list3,{term});
 >>;
 return list3;
 end;


procedure avoida(res,resgroeb,var);
begin
 scalar k,j,res2,dft,dfrg,term;
 res2:={};
 for j:=1:length(res) do
   << term:=part(res,j);
   for k:=1:length(var) do
   << dft:=df(term,part(var,k));
      dfrg:=df(part(resgroeb,k),part(var,k));
      if(dft*dfrg) neq 0 then term:=dft*part(resgroeb,k)-dfrg*term;
   >>;
   res2:=append(res2,{term});
 >>;
 return res2;
end;

end;