Artifact 06d799fc0d9e0ac11fafd7b17d45000b9abdbe7f493c6c70e02ac98f6b6de97a:
- Executable file
r38/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: 17401) [annotate] [blame] [check-ins using] [more...]
module diff; % Differentiation package. % Author: Anthony C. Hearn. % Modifications by: Francis J. Wright. % Copyright (c) 2000 Anthony C. Hearn. All rights reserved. fluid '(!*depend frlis!* powlis!* subfg!* wtl!*); fluid '(!*allowdfint !*dfint); global '(mcond!*); % Contains a reference to RPLACD (a table update), commented out. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Turning on the switch `allowdfint' allows "differentiation under the % integral sign", i.e. df(int(y, x), v) -> int(df(y, v), x), if this % results in a simplification. If the switch `dfint' is also turned % on then this happens regardless of whether the result simplifies. % Both switches are off by default. switch allowdfint, dfint; 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. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Consider df(u,x,y,z). % If none of x,y,z are equal to u then the order of differentiation is % commuted into a canonical form, unless the switch `nocommutedf' is % turned on, in which case the order of differentiation is not % commuted at all. The switch `nocommutedf' is off by default. % If one (or more) of x,y,z is equal to u then the order of % differentiation is NOT commuted and the derivative is NOT simplified % to zero, unless the switch `commutedf' is turned on. It is off by % default. (CRACK needs to turn it on!) % The new default behaviour should match the behaviour of REDUCE 3.6. % Turning on the switch `commutedf' should reproduce the default % behaviour of REDUCE 3.7. % If `commutedf' is off and the switch `simpnoncomdf' is on then % simplify df(u,x,u) -> df(u,x,2)/df(u,x), df(u,x,n,u) -> % df(u,x,n+1)/df(u,x), as suggested by Alain Moussiaux, PROVIDED u % depends only on the one variable x. This simplification removes the % non-commutative aspect of the derivative. switch commutedf, nocommutedf, simpnoncomdf; % Turning either `commutedf' or `nocommutedf' on turns the other off. % Turning commutation on or noncommutation off, or turning % simplification of noncommutative derivatives on, causes % resimplification. deflist('((commutedf ((t (off1 'nocommutedf) (rmsubs)))) (nocommutedf ((t (off1 'commutedf)) (nil (rmsubs)))) (simpnoncomdf ((t (rmsubs))))), 'simpfg); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % If the switch `expanddf' is turned on then REDUCE uses the chain % rule to expand symbolic derivatives of indirectly dependent % variables provided the result is unambiguous, i.e. provided there is % no direct dependence. It is off by default. Thus, for example, % given % depend f, u, v; depend {u, v}, x; % then, if `expanddf' is on, % df(f,x) -> df(f,u)*df(u,x) + df(f,v)*df(v,x) % whereas after % depend f, x; % df(f,x) does not expand at all (since the result would be ambiguous % and the algorithm would loop). % For similar handling in the case of explicit dependence, % e.g. df(f(u(x),v(x)),x), please use the standard package `DFPART' by % Herbert Melenk. switch expanddf; deflist('((expanddf ((t (rmsubs))))), 'simpfg); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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. % Take care with noncommuting expressions. if n>1 and noncomp u then return addsq(multsq(simpdf {u,v},simpexpt {u,n - 1}), multpq(u .** 1,diffp(u . (n - 1),v))) else 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; % Derivative of a dependent identifier; maybe apply chain % rule. Suppose u(v) = u(a(u),b(u),...), i.e. given % depend {u}, a, b, {a, b}, v; % then (essentially) depl!* = ((b v) (a v) (u b a)) if !*expanddf and not(v memq (x:=cdr atsoc(u, depl!*))) then << w := nil ./ 1; for each a in x do w := addsq(w, multsq(simp{'df,u,a},simp{'df,a,v})); go to e >>; 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 then << % multiple derivative if cadr u eq v then % (df (df v x y z ...) v) ==> 0 if commutedf if !*commutedf and null !*depend then return nil ./ 1 else if !*simpnoncomdf and (w:=atsoc(v, depl!*)) and null cddr w % and (cadr w eq (x:=caddr u)) then % (df (df v x) v) ==> (df v x 2)/(df v x) etc. % if single independent variable << x := caddr u; % w := simp {'quotient, {'df,u,x}, {'df,v,x}}; w := quotsq(simp{'df,u,x},simp{'df,v,x}); go to e >> else if 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!) 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 . merge!-ind!-vars(u,v), 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 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; 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=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 merge!-ind!-vars(u,v); % Consider (df u v) where u = (df a b c d ...) % It is non-commuting if a = v or if a in (b c d ...) % i.e. if a in (v b c d ...) if !*nocommutedf or (not !*commutedf and (cadr u memq (v . cddr u))) then derad!*(v,cddr u) else derad(v,cddr u); symbolic procedure derad!*(u,v); % Non-commuting derad %% Return the canonical list of differentiation variables %% equivalent to v,u, where v is a LIST of previus differentiation %% variables, when df(df(f(v,u), v), u) is simplified to %% df(f(v,u), v, u). Essentially just cons u onto v. reverse if u eq car(v:=reverse v) then % x,y, y 2 . v else if numberp car v and u eq cadr v then % x,y,n, y (car v + 1) . cdr v else u . v; % x,y, z 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;