Artifact 29bcc207f50531ff89f0de8e2f53dd18e56042e7ec8030cf6d68815ad455aba4:
- Executable file
r37/packages/scope/codhrn.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: 31053) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/scope/codhrn.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: 31053) [annotate] [blame] [check-ins using]
module ghorner; % Generalized Horner support. % ------------------------------------------------------------------- ; % Copyright : J.A. van Hulzen, Twente University, Dept. of Computer ; % Science, P.O.Box 217, 7500 AE Enschede, the Netherlands.; % Author : M.C. van Heerwaarden. ; % ------------------------------------------------------------------- ; % ------------------------------------------------------------------- ; % This module contains procedures which implement a generalized Horner; % scheme. There are two generalizations: ; % 1. It is possible to offer a set of assignment statements. Each RHS ; % will be transformed into a Horner scheme.; ; % 2. A list of identifiers is accepted as input.The polynomial will be; % hornered w.r.t. the first identifier in the list, then the ; % coefficients are hornered w.r.t. the second identifier, etc. ; % ; % The following steps are taken to achieve this result. ; % ; % The polynomial P is expanded by turning on the switch EXP and by ; % applying Aeval on P. Each term of the expanded polynomial is brought; % in a normal form. The terms are sorted using a binary tree represen-; % tation. From this tree a list of terms is extracted with the powers ; % in descending order.This list is rewritten into a Horner scheme. ; % ; % The 'normal form' of a term is: ; % (TIMES COEF (EXPT X N)) ; % It may be degenerated to: ; % (EXPT X N) for COEF = 1 ; % (TIMES COEF X) for N = 1 ; % (COEF) for N = 0 ; % When a term is a minus term, the minus is handled as a part of the ; % coefficient. ; % ------------------------------------------------------------------- ; global '(!*algpri autohorn); switch algpri; !*algpri := t; % ------------------------------------------------------------------- ; % ALGEBRAIC MODE COMMAND PARSER ; % ------------------------------------------------------------------- ; % The -STAT and FORM- procedures provide an interface with the ; % algebraic mode. To horner a set of expressions, one can use the ; % HORNER command, which has the following syntax: ; % <HORNER command> ::= GHORNER <ass. list> [VORDER <ID-list>] ; % <ass. list> ::= <assignment statement> | ; % << <assignment statement> ; % {; <assignment statement>} >>; % <ID-list> ::= <ID> | <ID> {, <ID>} ; % When the switch ALGPRI is ON, results are printed using Assgnpri, ; % When used inside a SCOPE-command the switch ALGPRI is turned OFF ; % automatically. However the current ALGPRI-setting is automatically ; % restored by SCOPE. ; % ------------------------------------------------------------------- ; put('ghorner, 'stat, 'ghornerstat); symbolic procedure ghornerstat; begin scalar x,y; % --------------------------------------------------------------- ; % GHORNER has already been read. ; % --------------------------------------------------------------- ; flag('(vorder), 'delim); flag('(!*rsqb), 'delim); if car(x := xread t) = 'progn % Read expressions; then x := cdr x % Remove keyword PROGN; else x := list x; % An assignment is also an asslist; if not(cursym!* eq 'vorder) then if cursym!* eq '!*semicol!* then autohorn := t else symerr('ghorner, t) else << autohorn := nil; y := remcomma xread nil % Read variable ordering list; >>; remflag('(vorder), 'delim); remflag('(!*rsqb), 'delim); return list('ghorner, x, y) end; put('ghorner, 'formfn, 'formghorner); symbolic procedure formghorner(u, vars, mode); list('ghorner, mkquote cadr u, mkquote caddr u); symbolic procedure ghorner(assset, varlist); % ------------------------------------------------------------------- ; % arg: assset = set of assignment statements ; % varlist = a list of variables ; % eff: For each assignment statement in assset, the RHS is turned into; % a Horner scheme. When varlist is not empty, the first variable ; % from varlist is used to form the scheme. The cdr of varlist is ; % used to transform the coefficients into a Horner scheme. ; % Implicitly, the assignment is executed by putting the SQ-form ; % of the Horner scheme on the property-list of the LHS-variable. ; % This means that the Horner scheme is available in the algebraic; % mode. When the switch ALGPRI is ON, the list of assignment ; % statements is printed. ; % res: If ALGPRI is OFF the list with hornered assignment statements ; % is returned. Nothing is returned when ALGPRI is ON. ; % ------------------------------------------------------------------- ; begin scalar h, hexp, res; hexp := !*exp; !*exp := nil; res := for each ass in assset collect if not eqcar(ass, 'setq) then rederr("Assignment statement expected") else << h:=inithorner(caddr ass, varlist); if !*algpri then << if eqcar(h, 'quotient) then put(cadr ass,'avalue, list('scalar, mk!*sq(numr !*f2q !*a2f cadr h ./ numr !*f2q !*a2f caddr h))) else put(cadr ass,'avalue, list('scalar, mk!*sq !*f2q !*a2f h)); assgnpri(h, list cadr ass, t); terpri() >> else list(car ass,cadr ass,h) >>; autohorn := nil; !*exp := hexp; if not !*algpri then return res end; symbolic procedure inithorner(p, varlist); % ------------------------------------------------------------------- ; % arg: p = polynomial ; % varlist = list of variables ; % eff: p is expanded and hornered to the various variables ; % res: the hornered version of p ; % ------------------------------------------------------------------- ; begin scalar n, d, hmcd, res; hmcd := !*mcd; !*mcd := t; p := reval p; res := hornersums(p, varlist); !*mcd := hmcd; return res end; symbolic procedure hornersums(p, varlist); if (atom(p) or domprop(p)) % JB 9/3/94 then p else if eqcar(p, 'plus) then horner(p, varlist) else append(list car p, for each elt in cdr p collect hornersums(elt, varlist)); symbolic procedure horner(p, varlist); % ------------------------------------------------------------------- ; % arg: p = polynomial ; % varlist = a list of variables for which the scheme must be made; % res: A Horner scheme of p with respect to first variable in varlist ; % ------------------------------------------------------------------- ; begin scalar hexp, tree, var; hexp := !*exp; !*exp := t; p := reval p; tree := '(nil nil nil); var := if varlist then car varlist else if autohorn then mainvar2 p else nil; if var then << for each kterm in cdr p do tree := puttree(tree, orderterm(kterm, var), var); p := gathertree(tree, var . cdr varlist); p := schema(p, var, kpow(car p, var)) >>; !*exp := hexp; return p end; symbolic procedure hornercoef(term, varlist); % ------------------------------------------------------------------- ; % arg: term = term of a polynomial in 'normal form' ; % varlist = the list of variables, including the one which just ; % has been used to create the scheme. ; % res: The same term is returned, but the coefficient has been turned ; % into a Horner scheme, using the second variable of varlist as ; % main variable. ; % ------------------------------------------------------------------- ; begin scalar n, cof; return if null(cof := kcof(term, (n := kpow(term, car varlist)))) then nil else if atom cof then term else if n = 0 then hornersums(cof, cdr varlist) else list(car term, hornersums(cof, cdr varlist), caddr term) end; symbolic procedure puttree(tree, term, var); % ------------------------------------------------------------------- ; % arg: tree = tree structure ( node, left edge, right edge), in which ; % the ordered terms are present. ; % term = the term which has to be put in ; % var = the variable for which the Horner scheme must be made ; % res: When the power of term is higher than the power of the node of ; % the root, puttree is called to place term in the right subtree ; % If the power is lower, term is placed in the left subtree. If ; % the powers are equal the coefficients are added. ; % ------------------------------------------------------------------- ; begin scalar c, n, m; return if null tree or null car tree then list (term, nil, nil) else if (n := kpow(term, var)) < (m := kpow(car tree, var)) then list(car tree, puttree(cadr tree, term, var), caddr tree) else if n > m then list(car tree, cadr tree, puttree(caddr tree, term, var)) else << % n = m so at least one term has been ; % inserted. Terms are added using only ; % one plus. ; c := kcof(car tree, n); if pairp c and car c = 'plus then c := cdr c else c := list c; if n = 0 then list(append('(plus), append(list(kcof(term,n)),c)), cadr tree, caddr tree) else list(list('times, append('(plus), append(list(kcof(term,n)),c)), if car c = 1 then car tree else caddar tree ), cadr tree, caddr tree)>> end; symbolic procedure gathertree(tree, varlist); % ------------------------------------------------------------------- ; % arg: tree = a tree as made by puttree ; % varlist = list of variables ; % res: a list of the terms which are stored in the tree. The term with; % the highest power is first in the list. For every term found, a; % Horner-scheme is made for the coefficients of this term.At this; % point the current variable remains on varlist. ; % ------------------------------------------------------------------- ; begin % This is a reversed depth-first search; return if null tree then nil else append(gathertree(caddr tree, varlist), append(list hornercoef(car tree, varlist), gathertree(cadr tree, varlist))) end; symbolic procedure orderterm(tt, var); % ------------------------------------------------------------------- ; % arg: tt = one term from the expanded polynomial ; % var = the variable for which the Horner scheme must be made ; % res: the term tt is returned in the 'normal form' which is described; % in the opening section. ; % ------------------------------------------------------------------- ; begin scalar h, res, factr, min; min := nil; if tt = var then res := tt else << if eqcar(tt, 'minus) then << min := t; tt := cadr tt >>; if not eqcar(tt,'times) then if min then if tt=var or (eqcar(tt,'expt) and cadr tt=var) then res := list('times, '(minus 1), tt) else res := list('minus, tt) else res := tt else << while not null (tt := cdr tt) do << if pairp(h := car tt) and eqcar(h, 'minus) then << min := not min; h := cadr h >>; if h = var then factr := h else << if eqcar(h, 'expt) and cadr h = var then factr := h else res := append(res, list h) >> >>; if min then << h := list('minus, car res); if null cdr res then res := list h else res := append(list h, cdr res) >>; res := if null factr then cons('times, res) else if null cdr res then list('times, car res, factr) else list('times, append('(times), res), factr) >> >>; return res end; symbolic procedure schema(pn, var, n); % ------------------------------------------------------------------- ; % arg: pn = the polynomial pn given as a list of terms in 'normal ; % form' in decsending order w.r.t. the powers of these ; % terms. ; % var = the Horner-scheme variable. ; % n = degree of the polynomial. ; % eff: Some effort is made to change "(TIMES var 1)" to "var" and to ; % turn "...(TIMES var (TIMES var..." into ; % "...(TIMES (EXPT var n) ..." ; % res: Horner scheme for the polynomial pn. ; % ------------------------------------------------------------------- ; begin scalar hn, k, k!+1mis; hn := kcof(car pn, n); % The n-th term always exists; if null (pn := cdr pn) then pn:=list(nil); % Else car(NIL) could be evaluated. for k := (n - 1) step -1 until 0 do << % --------------------------------------------------------- ; % hn contains the coefficients for the terms var^n upto ; % var^(k+1). The var for term var^(k+1) is still missing. ; % This is correct when for k=0 the last var will be added. ; % --------------------------------------------------------- ; if kpow(car pn, var) = k then << % k-th term exists; hn := list('plus, kcof(car pn, k), if hn = 1 then var else if not (k = n-1) and k!+1mis then if pairp hn and car hn = 'times then list('times,list('expt,var, kpow(cadr hn, var) + 1), caddr hn) else list('expt,var, kpow(hn, var) + 1) else list('times, var, hn) ); k!+1mis := nil; if null (pn := cdr pn) then pn:=list(nil) >> else << % k-th term misses; hn := if hn = 1 then var else if not (k = n-1) and k!+1mis then if pairp hn and car hn = 'times then list('times,list('expt,var, kpow(cadr hn, var) + 1), caddr hn) else list('expt, var, kpow(hn, var) + 1) else list('times, var, hn); k!+1mis := t >> >>; return hn end; symbolic procedure kpow(term, var); % ------------------------------------------------------------------- ; % arg: term = term of a polynomial in 'normal form' ; % var = the variable for which the Horner scheme must be made ; % res: the power of var in term ; % ------------------------------------------------------------------- ; begin scalar h; return if null term then nil else if (h := term) = var then 1 else if eqcar(h, 'expt) and eqcar(cdr h, var) then caddr h else if eqcar(h, 'times) then if (h := caddr h) = var then 1 else if not atom h and eqcar(cdr h, var) then caddr h else 0 else 0 end; symbolic procedure kcof(term, n); % ------------------------------------------------------------------- ; % arg: term = term of a polynomial in 'normal form' ; % n = the power of term ; % res: the coefficient of var in term ; % ------------------------------------------------------------------- ; if null n then nil else if n = 0 then term else if n = 1 then if not eqcar(term, 'times) then 1 else cadr term else if eqcar(term, 'expt) then 1 else cadr term; symbolic procedure mainvar2 u; % ------------------------------------------------------------------- ; % Same procedure as mainvar from ALG2.RED, but returns NIL instead of ; % 0 and does not allow a mainvar of the form (EXPT E X) to be returned; % ------------------------------------------------------------------- ; begin scalar res; res := if domainp(u := numr simp!* u) then nil else if sfp(u := mvar u) then prepf u else u; if eqcar(res, 'expt) then res := nil; return res end; %----------------------------------------------------------------------- % Algebraic mode psop function definition. %----------------------------------------------------------------------- symbolic procedure alghornereval u; % -------------------------------------------------------------------- ; % Variant of ghorner-command. Accepts 1 or 2 arguments. The first has % to be a list of equations. Their rhs's have to be hornered. The % second argument is optional. It defines the list of identifiers to % be used for the ordering. % -------------------------------------------------------------------- ; begin scalar algpri,assset,res,varlist; integer nargs; nargs:=length u; if nargs<3 then << assset:=foreach el in (if atom car u then cdr reval car u else cdar u )collect list('setq,cadr el,caddr el); if nargs=2 then varlist:=cdadr u >> else assset:='!*!*error!*!*; if eq(assset,'!*!*error!*!*) then rederr("WRONG NUMBER OF ARGUMENTS ALGHORNER") else << algpri:=!*algpri; !*algpri:=nil; res:=apply('ghorner,list(assset,varlist)); if (!*algpri:=algpri) then return algresults1(foreach el in res collect cons(cadr el,caddr el)) else return res >> end; put('alghorner,'psopfn,'alghornereval)$ %------------------------------------------------------------------ % Construction of Krukyov Horner's form of polynomial % JB 9/3/94 %------------------------------------------------------------------ algebraic procedure horner0(p,x)$ %---------------------------------------------------------- % p is a polynomial, % x is a Horner's variable$ % return p transformed to Horner's form %---------------------------------------------------------- begin scalar c,h$ on exp$ p:=p$ c:=reverse coeff(p,x)$ off exp$ h:=0$ while c neq {} do << h:=h*x+first c$ c:=rest c$ >>$ return h$ end$ algebraic procedure horner1(p)$ %---------------------------------------------------------- % p is a polynomial, % return p transformed to Horner's form % the MAINVAR of p use as a Horner's variable %---------------------------------------------------------- if numberp p then p else begin scalar c,h,x$ on exp$ p:=p$ x:=mainvar p$ c:=reverse coeff(p,x)$ off exp$ h:=0$ while c neq {} do << h:=h*x+horner1 first c$ c:=rest c$ >>$ return h$ end$ lisp global '(hvlst)$ % use for debug purposes. algebraic procedure horner2(p)$ %---------------------------------------------------------- % p is a polynomial, % return p transformed to Horner's form % Horner's variable is defined by HVAR1 procedure. % Outside effect: clear HVLST variable. % HVLST variable content the result of work of HVAR1 % (use for debug purposes). %---------------------------------------------------------- << clhvlist()$ horner20 p >>$ algebraic procedure horner20(p)$ %---------------------------------------------------------- % p is a polynomial, % return p transformed to Horner's form % Horner's variable is defined by HVAR1 procedure. %---------------------------------------------------------- if numberp p then p else begin scalar q,x,c$ on exp$ q:=p$ x:=hvar1 q$ c:=sub(x=0,q)$ q:=(q-c)/x$ off exp$ q:=horner20(q)*x+horner20(c)$ return q$ end$ symbolic procedure hvar1 q$ %---------------------------------------------------------- % q is a polynomial, % return Horner's variable. % Outside effect: set HVLST variable (use for debug only). % HVLST::=((expr . alst)...) % expr::=polynomial % alst::=((var.number)...) % Here the Horner's variable define as a variable % that entry to q in more tems then others. % For example: X+X**2+Y+1. The Horner's variable is X. %---------------------------------------------------------- if numberp q then rederr "HVAR1: impossible!" else begin scalar x,y,v$ q:=reval q$ % usually it is not needed. if null atom q and car q eq 'plus then q:=cdr q else q:=list q$ for each z in q do << if null atom z and car z eq 'minus then z:=cadr z$ if null atom z and car z eq 'times then z:=cdr z else z:=list z$ for each w in z do << if null atom w and car w eq 'expt then w:=cadr w else if numberp w then w:=nil$ if w and (y:=assoc(w,v)) then rplacD(y,cdr y + 1) else if w then v:=(w . 1).v$ >>$ >>$ x:=car v$ for each z in cdr v do if cdr z > cdr x then x:=z$ hvlst:=(q.v).hvlst$ return car x$ end$ algebraic procedure khorner20(p,vlst)$ %---------------------------------------------------------- % p is a polynomial, vlst is a list of horner-variables. % return p transformed to Horner's form % Horner's variable is defined by the khvar1-procedure. %---------------------------------------------------------- if numberp p then p else begin scalar q,x,c; on exp; q:=p; if (x:=khvar1(q,vlst)) then << c:=sub(x=0,q); q:=(q-c)/x; off exp; return(khorner20(q,vlst)*x+khorner20(c,vlst)) >> else << off exp; return(nestedfac q) >> end$ symbolic procedure khvar1(q,vlst); %---------------------------------------------------------- % q is a polynomial, vlst is a list of horner-variables. % return Horner's variable. % Here the Horner's variable is defined as a variable % that occurs in more q-terms than the others. % For example: X in q = X+X**2+Y+1. %---------------------------------------------------------- if numberp q then rederr "HVAR1: impossible!" else begin scalar x,y,v; vlst:=cdr vlst; q:=reval q;% redefinition q usually not needed. if null atom q and car q eq 'plus then q:=cdr q else q:=list q; foreach z in q do << if null atom z and car z eq 'minus then z:=cadr z; if null atom z and car z eq 'times then z:=cdr z else z:=list z; for each w in z do << if null atom w and car w eq 'expt then w:=cadr w else if numberp w then w:=nil; if w and memq(w,vlst) then if (y:=assoc(w,v)) then rplacd(y,cdr y + 1) else v:=(w . 1).v >> >>; if v then << x:=car v; foreach z in cdr v do if cdr z > cdr x then x:=z$ return car x >> else return nil end$ symbolic procedure hvlist()$ %---------------------------------------------------------- % Procedure for printing HVLST variable. % Debug utility. %---------------------------------------------------------- for each x in hvlst do print x$ symbolic procedure clhvlist()$ %---------------------------------------------------------- % Procedure for clearing HVLST variable. % Debug utility. %---------------------------------------------------------- hvlst:=nil$ symbolic operator khvar1,hvar1,hvlist,clhvlist$ % Interface with REDUCE % ------------------------------------------------------------------- % Interface for generalised facilities, based on the use of the % procedure gkhorner. This procedure can be used with one argument % only, being a list of equations of the form lhsi=rhsi, where the % i-th lhs is a name and the i-th rhs a (multivariate) polynomial, % to be hornered either exhaustively using horner20(rhsi) or restric- % tively using the second argument vlst, being a list of horner- % variables, and procedure khorner20. When further vlst variables are % absent the remaining parts of q are further polished using nestedfac. % ------------------------------------------------------------------- symbolic procedure khornereval u; begin scalar poly,varlst; integer nargs; nargs:=length u; if nargs<3 then << poly:=aeval car u; if nargs=2 then varlst:=aeval cadr u>> else poly:='!*!*error!*!*; if eq(poly,'!*!*error!*!*) then rederr("WRONG NUMBER OF ARGUMENTS KHORNER") else return if nargs=1 then reval horner2 poly else reval khorner20(poly,varlst) end; put('khorner,'psopfn,'khornereval)$ symbolic procedure gkhornereval u; begin scalar poly_s,varlst; integer nargs; nargs:=length u; if nargs<3 then << poly_s:=cdar u; if nargs=2 then varlst:=cadr u>> else poly_s:='!*!*error!*!*; if eq(poly_s,'!*!*error!*!*) then rederr("WRONG NUMBER OF ARGUMENTS GKHORNER") else return if pairp(car poly_s) and eq(caar poly_s,'equal) then append(list('list), foreach poly in poly_s collect list('equal, cadr poly, khornereval if nargs=1 then cddr poly else list(caddr poly,varlst))) else append(list('list), foreach poly in poly_s collect khornereval if nargs=1 then list(poly) else list(poly,varlst)) end$ put('gkhorner,'psopfn,'gkhornereval)$ symbolic procedure alggkhornereval u; begin scalar poly_s,varlst; integer nargs; nargs:=length u; if nargs<3 then << poly_s:=cdar u; if nargs=2 then varlst:=cadr u >> else poly_s:='!*!*error!*!*; if eq(poly_s,'!*!*error!*!*) then rederr("WRONG NUMBER OF ARGUMENTS GKHORNER") else return algresults1(foreach poly in poly_s collect cons(cadr poly, khornereval if nargs=1 then cddr poly else list(caddr poly,varlst))) end; put('alggkhorner,'psopfn,'alggkhornereval)$ endmodule; end;