File r37/packages/specfn/linrec.red from the latest check-in

module linrec;

% solves (simple in) homogenous linear recursion relation REC which
% has to have the form var(n)... with initial conditions
% IC in the form of list {var(1)=?,var(2)=...}
% the method used - substitution s(n)=x^n
% like e^(lambda*x) in linear homogenous ODE.

% The inhomogenous term has to be a constant.

% The following code has been developed by Richard Liska from Prague
% for the "Barry Simon Tests for PC Magazine".
% Some checks and generalizations included by W. Neun , ZIB Berlin

% the switch trlinrec turns the verbose mode on.

switch trlinrec;

algebraic procedure rsolve (rec,ic);

begin scalar mvar1,inde,constpart,lowestind,modic,modrec,indis;

mvar1 := part(mainvar rec ,0);
inde  := mainvar part(mainvar rec,1);
indis := the_indices(foreach kk in ic collect lhs kk,mvar1);
highestind := if indis neq {} then max(indis) else NIL;
nindis := the_indices(rec,mvar1);

           % applying rule : rec where r(~n) => 0

constpart := lisp aeval list('whereexp,
                list('list,list('replaceby,list(mvar1,list('!~,'n)),0)),
                rec);

if not freeof(constpart,inde) then rederr
   ("Cant solve recurrence equations with non-constant coefficients");

if eq(constpart , 0) then return
        solve_lin_rec(rec,ic)
  else
        <<  modrec := sub (inde = inde +1,rec);
            modrec := modrec -rec;
	    if (highestind) then <<
% Propagate the recursion to get additional start value
            modic := sub(inde=inde -max(nindis)+ highestind+1,rec);
            for each aa in ic do modic := sub(aa,modic);
            modic := (first solve(modic,mainvar modic)) . ic >>
            else modic := ic;
            return solve_lin_rec(modrec,modic);
         >>;
end;


fluid '(linrecx!*,linrecvar!*);
 
algebraic procedure solve_lin_rec(rec,ic);
   % solves homogenous linear recursion relation REC which
   % has to have the form var(n)... with initial conditions
   % IC in the form of list {var(1)=?,var(2)=...}
   % the method used - substitution s(n)=x^n
   % like e^(lambda*x) in linear homogenous ODE.
   % Of course some checking should be added. (Done WN)
   begin
     scalar lrec,sol,msol,gsol,j,flagg,c,linrecvar!*,errflag,nsave;
     clear n;
     linrecvar!* := part (mainvar rec,0);
     linrecx!* := lisp gensym();
     c:= lisp mkquote gensym(); %this is the dirty part. WN
     operator c;
     if(part(rec,0) eq linrecvar!*) and arglength(ic)=1 and
		part(mainvar lhs first ic,0) = linrecvar!*
		then return rhs (first ic);
      if(part(rec,0) neq plus) then return  rederr
       "Cant solve recurrence equations with non-constant coefficients";
      lrec := arglength rec;
      lrec := part(rec,lrec);
     for all n let linrecvar!*(n) = linrecx!*^n;
     lrec := lrec;
     rec:=rec /lrec;
     for all n clear linrecvar!* (n);
     rec:=num rec;
     for each j in coeff(rec,linrecx!*) do
                if (not freeof(j,part(part (part rec,1),1)))
			then errflag := 17; %???
     if (errflag = 17) then return rederr
       "Cant solve recurrence equations with non-constant coefficients";
     j:=1;
     for each a in solve(rec,linrecx!*) do
       <<a:=rhs a^n;
         gsol:=gsol+c(j)*a;
         j:=j+1;
         msol:=first multiplicities!*;
         multiplicities!*:=rest multiplicities!*;
         while msol>1 do
           <<a:=n*a;
             gsol:=gsol+c(j)*a;
             j:=j+1;
             msol:=msol-1 >> >>;
     if lisp !*trlinrec then write "General solution: ",linrecvar!*
		,"(N) := ",gsol;
     if ic = {} then sol := {} else
     sol:=solve(for each a in ic collect 
			sub(n=part(lhs a,1),gsol)=rhs a,
                for i:=1:arglength ic collect c(i));
   % If some c(i) remains it can be arbitrary  complex;
     sol := lisp subla('((equal . replaceby)),sol);
     sol := lisp subla('((equal . replaceby)),sol);
     let sol;
     gsol:=gsol;
     clearrules sol;
     let moivre_expt;
     gsol:=gsol;
     clearrules moivre_expt;

     for i:=1:j do  if coeff(gsol,lisp list(c, i)) = list(gsol) then nil
	else gsol:=
	   sub(lisp list(c, i) = lisp caaar makearbcomplex(),gsol);
     return gsol
   end;

% (1 + i)**n =>     %applying Moivre's formula for complex numbers

%		 N/2       N*PI          N*PI
%		2   *(COS(------) + SIN(------)*I)
%		            4             4

algebraic (moivre_expt := { (~z)^(~k) =>
	      Moivre(z,k) when not  freeof(z,i)});

algebraic procedure Moivre(z,k);

  begin scalar rho,phi; % what ( will happen
  rho := sqrt( (repart z)^2 + (impart z)^2);
  if repart  z = 0 then phi := pi/2 else
  	phi := atan((impart z)/(repart  z));
  return rho^k *(cos(k*phi) + i * sin (k*phi));
  end;

algebraic procedure the_indices(ex,mvar);

if part(ex,0) = list then for each kk in ex join the_indices(kk,mvar)
else
begin scalar eqq,L1,L2,kern;
  eqq := ex;
  lisp (kern := union (kernels !*q2f  (numr simp eqq ./ 1),
                        kernels !*q2f (denr simp eqq ./ 1)));

  L1 := 'list . lisp  foreach k in kern join
	 if atom k then nil else if eqcar(k,mvar) then
                                list cadr k else nil;

  return L1;
end;
endmodule;
end;