Artifact 96288f155566b971f6613829e216b074438ab2190d3db43cd15685b0a24e63f2:
- Executable file
r37/packages/crack/crintfix.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: 13614) [annotate] [blame] [check-ins using] [more...]
module intfix$ % Further fixes to the integration package. 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$