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

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;