reduce-historical/r38/packages/poly/diff.red

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;