Artifact b37eb2b77ebfff79b8cfcd26d085092c560d1af3eba3beede89ef8b2291cb691:
- Executable file
r37/packages/specfn/linrec.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 5267) [annotate] [blame] [check-ins using] [more...]
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;