Artifact ea2094fdee5b4806b2665c64b0bd62996b44642c3f70fe1a4e44e808de3ca20a:
- Executable file
r37/packages/poly/diff.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: 12786) [annotate] [blame] [check-ins using] [more...]
module diff; % Differentiation package. % Author: Anthony C. Hearn. % Modifications by: Francis J. Wright. % Copyright (c) 1991 RAND. All rights reserved. fluid '(!*allowdfint !*depend !*fjwflag frlis!* powlis!* subfg!* wtl!*); switch allowdfint, dfint, fjwflag; !*fjwflag := t; % Let's try it on for a while. global '(mcond!*); % Contains a reference to RPLACD (a table update), commented out. symbolic procedure simpdf u; % U is a list of forms, the first an expression and the remainder % kernels and numbers. % Value is derivative of first form wrt rest of list. begin scalar v,x,y,z; if null subfg!* then return mksq('df . u,1); v := cdr u; u := simp!* car u; a: if null v or null numr u then return u; x := if null y or y=0 then simp!* car v else y; if denr x neq 1 or atom numr x then typerr(prepsq x,"kernel or integer") else (if domainp z then if get(car z,'domain!-diff!-fn) then begin scalar dmode!*,alglist!*; x := prepf z; if null prekernp x then typerr(x,"kernel") end else typerr(prepf z,"kernel") else if null red z and lc z = 1 and ldeg z = 1 then x := mvar z else typerr(prepf z,"kernel")) where z = numr x; v := cdr v; if null v then <<u := diffsq(u,x); go to a>>; y := simp!* car v; % At this point, y must be a kernel or equivalent to an integer. % Any other value is an error. if null numr y then <<v := cdr v; y := nil; go to a>> else if not(z := d2int y) then <<u := diffsq(u,x); go to a>>; v := cdr v; for i:=1:z do u := diffsq(u,x); y := nil; go to a end; symbolic procedure d2int u; if denr u neq 1 then nil else if numberp(u := numr u) then u else if not domainp u or not(car u eq '!:rd!:) then nil else (if abs(float x - u)<!!fleps1 then x else nil) where x=fix u where u=rd2fl u; put('df,'simpfn,'simpdf); symbolic procedure prekernp u; if atom u then idp u else idp car u and null((car u memq '(plus minus times quotient recip)) or ((car u eq 'expt) and fixp caddr u)); symbolic procedure diffsq(u,v); % U is a standard quotient, V a kernel. % Value is the standard quotient derivative of U wrt V. % Algorithm: df(x/y,z)= (x'-(x/y)*y')/y. multsq(addsq(difff(numr u,v),negsq multsq(u,difff(denr u,v))), 1 ./ denr u); symbolic procedure difff(u,v); %U is a standard form, V a kernel. %Value is the standard quotient derivative of U wrt V. % Allow for differentiable domains. if atom u then nil ./ 1 else if atom car u then (if diff!-fn then apply2(diff!-fn,u,v) else nil ./ 1) where (diff!-fn = get(car u,'domain!-diff!-fn)) else addsq(addsq(multpq(lpow u,difff(lc u,v)), multsq(diffp(lpow u,v),lc u ./ 1)), difff(red u,v)); symbolic procedure diffp(u,v); % U is a standard power, V a kernel. % Value is the standard quotient derivative of U wrt V. begin scalar n,w,x,y,z; integer m; n := cdr u; % integer power. u := car u; % main variable. if u eq v and (w := 1 ./ 1) then go to e else if atom u then go to f %else if (x := assoc(u,dsubl!*)) and (x := atsoc(v,cdr x)) % and (w := cdr x) then go to e % deriv known. % DSUBL!* not used for now. else if (not atom car u and (w:= difff(u,v))) or (car u eq '!*sq and (w:= diffsq(cadr u,v))) then go to c % extended kernel found. else if x := get(car u,'dfform) then return apply3(x,u,v,n) else if x:= get(car u,dfn_prop u) then nil else if car u eq 'plus and (w := diffsq(simp u,v)) then go to c else go to h; % unknown derivative. y := x; z := cdr u; a: w := diffsq(simp car z,v) . w; if caar w and null car y then go to h; % unknown deriv. y := cdr y; z := cdr z; if z and y then go to a else if z or y then go to h; % arguments do not match. y := reverse w; z := cdr u; w := nil ./ 1; b: % computation of kernel derivative. if caar y then w := addsq(multsq(car y,simp subla(pair(caar x,z), cdar x)), w); x := cdr x; y := cdr y; if y then go to b; c: % save calculated deriv in case it is used again. % if x := atsoc(u,dsubl!*) then go to d % else x := u . nil; % dsubl!* := x . dsubl!*; % d: rplacd(x,xadd(v . w,cdr x,t)); e: % allowance for power. % first check to see if kernel has weight. if (x := atsoc(u,wtl!*)) then w := multpq('k!* .** (-cdr x),w); m := n-1; % Evaluation is far more efficient if results are rationalized. return rationalizesq if n=1 then w else if flagp(dmode!*,'convert) and null(n := int!-equiv!-chk apply1(get(dmode!*,'i2d),n)) then nil ./ 1 else multsq(!*t2q((u .** m) .* n),w); f: % Check for possible unused substitution rule. if not depends(u,v) and (not (x:= atsoc(u,powlis!*)) or not depends(cadddr x,v)) and null !*depend then return nil ./ 1; w := list('df,u,v); w := if x := opmtch w then simp x else mksq(w,1); go to e; h: % Final check for possible kernel deriv. if car u eq 'df % multiple derivative then if depends(cadr u,v) % FJW - my version of above test was simply as follows. Surely, inner % derivative will already have simplied to 0 unless v depends on A! and not(cadr u eq v) % (df (df v A) v) ==> 0 %% and not(cadr u eq v and not depends(v,caddr u)) %% % (df (df v A) v) ==> 0 unless v depends on A. then <<if !*fjwflag and eqcar(cadr u, 'int) then % (df (df (int F x) A) v) ==> (df (df (int F x) v) A) ? % Commute the derivatives to differentiate the integral? if caddr cadr u eq v then % Evaluating (df u v) where u = (df (int F v) A) % Just return (df F A) - derivative absorbed << w := 'df . cadr cadr u . cddr u; go to j >> else if !*allowdfint and % Evaluating (df u v) where u = (df (int F x) A) % (If dfint is also on then this will not arise!) % Commute only if the result simplifies: not_df_p(w := diffsq(simp!* cadr cadr u, v)) then << % Generally must re-evaluate the integral (carefully!) % FJW. Bug fix! % w := aeval{'int, mk!*sq w, caddr cadr u} . cddr u; w := 'df . reval{'int, mk!*sq w, caddr cadr u} . cddr u; go to j >>; % derivative absorbed if (x := find_sub_df(w:= cadr u . derad(v,cddr u), get('df,'kvalue))) then <<w := simp car x; for each el in cdr x do for i := 1:cdr el do w := diffsq(w,car el); go to e>> else w := 'df . w >> else if null !*depend then return nil ./ 1 else w := {'df,u,v} else w := {'df,u,v}; j: if (x := opmtch w) then w := simp x else if not depends(u,v) and null !*depend then return nil ./ 1 else w := mksq(w,1); go to e end; deflist('((dfint ((t (rmsubs)))) (allowdfint ((t (progn (put 'int 'dfform 'dfform_int) (rmsubs))) (nil (remprop 'int 'dfform))))), 'simpfg); % There is no code to reverse the df-int commutation, % so no reason to call rmsubs when the switch is turned off. symbolic procedure dfform_int(u, v, n); % Simplify a SINGLE derivative of an integral. % u = '(int y x) [as main variable of SQ form] % v = kernel % n = integer power % Return SQ form of df(u**n, v) = n*u**(n-1)*df(u, v) % This routine is called by diffp via the hook % "if x := get(car u,'dfform) then return apply3(x,u,v,n)". % It does not necessarily need to use this hook, but it needs to be % called as an alternative to diffp so that the linearity of % differentiation has already been applied. begin scalar result, x, y; y := simp!* cadr u; % SQ form integrand x := caddr u; % kernel result := if v eq x then y % df(int(y,x), x) -> y replacing the let rule in INT.RED else if not !*intflag!* and % not in the integrator % If used in the integrator it can cause infinite loops, % e.g. in df(int(int(f,x),y),x) and df(int(int(f,x),y),y) !*allowdfint and % must be on for dfint to work << y := diffsq(y, v); !*dfint or not_df_p y >> % it has simplified then simp{'int, mk!*sq y, x} % MUST re-simplify it!!! % i.e. differentiate under the integral sign % df(int(y, x), v) -> int(df(y, v), x). % (Perhaps I should use prepsq - kernels are normally true prefix?) else !*kk2q{'df, u, v}; % remain unchanged if not(n eq 1) then result := multsq( (((u .** (n-1)) .* n) .+ nil) ./ 1, result); return result end; symbolic procedure not_df_p y; % True if the SQ form y is not a df kernel. not(denr y eq 1 and not domainp (y := numr y) and eqcar(mvar y, 'df)); % Compute a dfn-property name corresponding to the argument number % of an operator expression. Here we assume that most functions % will have not more than 3 arguments. symbolic procedure dfn_prop(w); (if n=1 then 'dfn else if n=2 then 'dfn2 else if n=3 then 'dfn3 else mkid('dfn,n)) where n=length cdr w; % The following three functions, and the hooks to this code above, were % suggested by Gerhard Post and Marcel Roelofs. symbolic procedure find_sub_df(df_args,df_values); df_values and (is_sub_df(df_args,car df_values) or find_sub_df(df_args,cdr df_values)); symbolic procedure is_sub_df(df_args,df_value); begin scalar df_set,kernel,n,entry; if car(df_args) neq cadar(df_value) then return nil; % check fns. df_args := dot_df_args cdr df_args; df_set := cddar df_value; while df_set and df_args do % Check differentiations. <<kernel := car df_set; if cdr df_set and fixp(n := cadr df_set) then df_set := cdr df_set else n := 1; if (entry := atsoc(kernel,df_args)) and (n := cdr entry-n) geq 0 then rplacd(entry,n) else df_args:=nil; df_set := cdr df_set>>; return if df_args then (cadr(df_value) . df_args); end; symbolic procedure dot_df_args l; begin scalar kernel,n,df_args; while l do <<kernel := car l; if cdr l and fixp(n := cadr l) then l := cdr l else n := 1; df_args := (kernel . n) . df_args; l := cdr l>>; return df_args; end; symbolic procedure derad(u,v); if null v then list u else if numberp car v then car v . derad(u,cdr v) else if u=car v then if cdr v and numberp cadr v then u . (cadr v + 1) . cddr v else u . 2 . cdr v else if ordp(u,car v) then u . v else car v . derad(u,cdr v); symbolic procedure letdf(u,v,w,x,b); begin scalar y,z,dfn; if atom cadr x then go to b else if not idp caadr x then typerr(caadr x,"operator") else if not get(caadr x,'simpfn) then <<redmsg(caadr x,"operator"); mkop caadr x>>; rmsubs(); dfn := dfn_prop cadr x; if not(mcond!* eq 't) or not frlp cdadr x or null cddr x or cdddr x or not frlp cddr x or not idlistp cdadr x or repeats cdadr x or not(caddr x member cdadr x) then go to b; z := lpos(caddr x,cdadr x); if not get(caadr x,dfn) then put(caadr x, dfn, nlist(nil,length cdadr x)); w := get(caadr x,dfn); if length w neq length cdadr x then rerror(poly,17, list("Incompatible DF rule argument length for", caadr x)); a: if null w or z=0 then return errpri1 u else if z neq 1 then <<y := car w . y; w := cdr w; z := z-1; go to a>> else if null b then y := append(reverse y,nil . cdr w) else y := append(reverse y,(cdadr x . v) . cdr w); return put(caadr x,dfn,y); b: %check for dependency; if caddr x memq frlis!* then return nil else if idp cadr x and not(cadr x memq frlis!*) then depend1(cadr x,caddr x,t) else if not atom cadr x and idp caadr x and frlp cdadr x then depend1(caadr x,caddr x,t); return nil end; symbolic procedure frlp u; null u or (car u memq frlis!* and frlp cdr u); symbolic procedure lpos(u,v); if u eq car v then 1 else lpos(u,cdr v)+1; endmodule; end;