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;