reduce-historical/r38/packages/crack/crintfix.red

%********************************************************************
module intfix$  % Further fixes to the integration package.
%********************************************************************
%  Routines to extend the REDUCE integrator or to fix problems
%  Author: Francis Wright
%
%  $Id$
%

if lisp !*comp then apply1('load!-package, 'int)$   

fluid '(!*depend !*nolnr !*failhard)$ 

% die folgende Aenderung verhindert das Erzeugen von int* ...

remd('simpint!*)$

symbolic procedure simpint!* u$
   begin scalar x$
      return if (x := opmtch('int . u)) then simp x
              else simpiden('int . u)
% statt       else simpiden('int!* . u)
   end$

% ein Patch fuer das REDUCE 3.5 EZGCD

%symbolic procedure simpexpt u$
%   % We suppress reordering during exponent evaluation, otherwise
%   % internal parts (as in e^(a*b)) can have wrong order.
%   begin scalar expon;
%      expon := simpexpon carx(cdr u,'expt) where kord!*=nil;
%      expon := resimp expon;  % We still need right order. <--- change.
%      return simpexpt1(car u,expon,nil)
%   end$

% Zum Integrieren

% put('int, 'simpfn, 'SimpIntPatch)$

%algebraic <<
% % fuer reelle Rechnungen:
% let {abs(~r)**(~n) => r**n when (fixp(n) and evenp(n))}$
% let {
% int(1/~x^(~n),~x) => -x/(x^n*(n-1)) when numberp n,
%     ~x^(~m/~n)*~x => x**((m+n)/n) when (numberp n and numberp m),
%     int(~z/~y,~x) => log(y) when z = df(y,x)}$
%
% if sin(!%x)**2+cos(!%x)**2 neq 1 then
% let {sin(~x)**2 => 1-cos(x)**2}$
%
% if cosh(!%x)**2 neq (sinh(!%x)**2 + 1) then
% let {cosh(~x)**2 => (sinh(x)**2 + 1)}$
%
% if sin(!%x)*tan(!%x/2)+cos(!%x) neq 1 then
% let {tan(~x/2) => (1-cos(x))/sin(x)}$
%
% if sin(!%x)*cot(!%x/2)-cos(!%x) neq 1 then
% let {cot(~x/2) => (1+cos(x))/sin(x)}$
%
% if sqrt(!%x**2-!%y**2)-sqrt(!%x-!%y)*sqrt(!%x+!%y) neq 0 then
% let {sqrt(~x)*sqrt(~y) => sqrt(x*y)}
%>>$

endmodule$


module dfint$

% Patch to improve differentiation, mainly of integrals.
% This version specifically for use by the crack package.

% Francis J. Wright <F.J.Wright@QMW.ac.uk>, 27 December 1997

fluid '(!*fjwflag)$  !*fjwflag := t$

switch allowdfint, dfint$               % dfint OFF by default
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.

!*allowdfint := t$                      % allowdfint ON by default
put('int, 'dfform, 'dfform_int)$

% The switch allowdfint ALLOWS differentiation under the integral sign
% provided the result simplies, and should normally be on.

% The switch dfint FORCES differentiation under the integral sign,
% PROVIDED ALLOWDFINT IS ALSO ON, and should normally be turned on
% only when required.


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$


% Author: Francis J. Wright <F.J.Wright@QMW.ac.uk>
% Last revised: 27 December 1997

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))$

endmodule$


module intdf$

% Patch to simpint1 in src/int/trans/driver.red to provide better
% simplification of integrals of derivatives.  (I think -- hope --
% this is the right place to hook this patch into the integrator!)
% This patch was motivated by the needs of crack.

% F.J.Wright@Maths.QMW.ac.uk, 31 December 1997

%% load_package int$
%apply1('load!-package, 'int)$           % not at compile time!

switch PartialIntDf$                    % off by default
deflist('((PartialIntDf ((t (rmsubs))))), 'simpfg)$

% If the switch PartialIntDf is turned on then integration by parts is
% performed if the result simplifies in the sense that it integrates a
% symbolic derivative and does not introduce new symbolic derivatives.
% However, because the initial integral contains an unevaluated
% derivative then the result must still contain an unevaluated
% integral.

symbolic procedure simpint1 u;
   % Varstack* rebound, since FORMLNR use can create recursive
   % evaluations.  (E.g., with int(cos(x)/x**2,x)).
   begin scalar !*keepsqrts,v,varstack!*;
      u := 'int . prepsq car u . cdr u;
      if (v := formlnr u) neq u
        then if !*nolnr
	       then <<v := simp subst('int!*,'int,v);
                      return remakesf numr v ./ remakesf denr v>>
	      else <<!*nolnr := nil . !*nolnr;
		     v:=errorset!*(list('simp,mkquote v),!*backtrace);
		     if pairp v then v := car v else v := simp u;
		     !*nolnr := cdr !*nolnr; 
                     return v>>;
      % FJW: At this point linearity has been applied.
      return if (v := opmtch u) then simp v
         % FJW: Check for a directly integrable derivative:
      else if (v := NestedIntDf(cadr u, caddr u)) then mksq(v,1)
      else if !*failhard then rerror(int,4,"FAILHARD switch set")
         % FJW: Integrate by parts if the result simplifies:
      else if !*PartialIntDf and
         (v := PartialIntDf(cadr u, caddr u)) then mksq(v,1)
      else mksq(u,1)
   end$

symbolic procedure NestedIntDf(y, x);
   %% int( ... df(f,A,x,B) ..., x) -> ... df(f,A,B) ...
   %% Find a df(f,A,x,B) among possibly nested int's and df's within
   %% the integrand y in int(y,x), and return the whole structure y
   %% but with the derivative integrated; otherwise return nil.
   %% [A,B are arbitrary sequences of kernels.]
   not atom y and
   begin scalar car_y, nested;
      return
         if (car_y := car y) eq 'df and memq(x, cddr y) then
            %% int( df(f, A, x, B), x ) -> df(f, A, B)
            'df . cadr y . delete(x, cddr y)
               %% use delete for portability!
               %% deleq is defined in CSL, delq in PSL -- oops!
         else if memq(car_y, '(df int)) and
            (nested := NestedIntDf(cadr y, x)) then
            %% int( df(int(df(f, A, x, B), c), C), x ) ->
            %%      df(int(df(f, A, B), c), C)
            %% int( int(df(f, A, x, B), c), x ) ->
            %%      int(df(f, A, B), c)
            car_y . nested . cddr y
   end$

symbolic procedure PartialIntDf(y, x);
   %% int(u(x)*df(v(x),x), x) -> u(x)*v(x) - int(df(u(x),x)*v(x), x)
   %% Integrate by parts if the resulting integral simplifies [to
   %% avoid infinite loops], which means that df(u(x),x) may not
   %% contain any unevaluated derivatives; otherwise return nil.
   not atom y and
   begin scalar denlist, facs, df, u, v;
      if car y eq 'quotient then <<
         denlist := cddr y;
         % y := numerator:
         if atom(y := cadr y) then return % no derivative
      >>;
      % y := list of factors:
      if car y eq 'times then y := cdr y
      else if denlist then y := y . nil
      else return;
      % Find an integrable derivative among the factors:
      facs := y;
      while facs and not
         (eqcar(df := car facs, 'df) and memq(x, cddr df)) do
            facs := cdr facs;
      if null facs then return;         % no integrable derivative
      % Construct u(x) and v(x) [v(x) may still be a derivative]:
      u := delete(df, y);               % list of factors
      u := if null u then 1 else if cdr u then 'times . u else car u;
      if denlist then u := 'quotient . u . denlist;
      v := cadr df;                     % kernel being differentiated
      if (df := delete(x, cddr df)) then v := 'df . v . df;
      % Check that df(u(x),x) simplifies:
      if smemq('df, df := reval {'df,u,x}) then return;
      return reval {'difference,
         {'times,u,v}, {'int, {'times, df, v}, x}}
   end$

endmodule$

end$