Artifact 36a6c0a08a98064d7c6e9123e28f064afadd05960b96798785e6dd43006e74b8:
- Executable file
r37/packages/int/driver.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: 30166) [annotate] [blame] [check-ins using] [more...]
module driver; % Driving routines for integration program. % Author: Mary Ann Moore and Arthur C. Norman. % Modifications by: John P. Fitch, David Hartley, Francis J. Wright. fluid '(!*algint !*backtrace !*exp % !*failhard !*gcd !*intflag!* !*keepsqrts !*limitedfactors !*mcd !*noncomp !*nolnr !*partialintdf !*precise !*purerisch !*rationalize !*structure !*trdint !*trint !*uncached basic!-listofnewsqrts basic!-listofallsqrts gaussiani intvar kord!* listofnewsqrts listofallsqrts loglist powlis!* sqrt!-intvar sqrt!-places!-alist subfg!* varlist varstack!* xlogs zlist); exports integratesq,simpint,simpint1; imports algebraiccase,algfnpl,findzvars,getvariables,interr,printsq, transcendentalcase,varsinlist,kernp,simpcar,prepsq,mksq,simp, opmtch,formlnr; switch algint,nolnr,trdint,trint; switch hyperbolic; % Form is int(expr,var,x1,x2,...); % meaning is integrate expr wrt var, given that the result may % contain logs of x1,x2,... % x1, etc are intended for use when the system has to be helped % in the case that expr is algebraic. % Extended arguments x1, x2, etc., are not currently supported. symbolic procedure simpint u; % Simplifies an integral. First two components of U are the integrand % and integration variable respectively. Optional succeeding % components are log forms for the final integral. if atom u or null cdr u or cddr u and (null cdddr u or cddddr u) then rerror(int,1,"Improper number of arguments to INT") else if cddr u then simpdint u % then if getd 'simpdint then simpdint u % else rerror(int,2,"Improper number of arguments to INT") else begin scalar ans,dmod,expression,variable,loglist,oldvarstack, !*intflag!*,!*purerisch,cflag,intvar,listofnewsqrts, listofallsqrts,sqrtfn,sqrt!-intvar,sqrt!-places!-alist, basic!-listofallsqrts,basic!-listofnewsqrts,coefft, varchange,w; !*intflag!* := t; % Shows we are in integrator. variable := !*a2k cadr u; if not(idp variable or pairp variable and numlistp cdr variable) % then typerr(variable,"integration variable"); then <<varchange := variable . intern gensym(); if !*trint then printc {"Integration kernel", variable, "replaced by simple variable", cdr varchange}; variable := cdr varchange>>; intvar := variable; % Used in SIMPSQRT and algebraic integrator. w := cddr u; if w then rerror(int,3,"Too many arguments to INT"); listofnewsqrts:= list mvar gaussiani; % Initialize for SIMPSQRT. listofallsqrts:= list (argof mvar gaussiani . gaussiani); sqrtfn := get('sqrt,'simpfn); put('sqrt,'simpfn,'proper!-simpsqrt); % We need explicit settings of several switches during integral % evaluation. In addition, the current code cannot handle domains % like floating point, so we suppress it while the integral is % calculated. UNCACHED is turned on since integrator does its own % caching. % Any changes made to these settings must also be made in wstrass. if dmode!* then << % added by Alan Barnes if (cflag:=get(dmode!*, 'cmpxfn)) then onoff('complex, nil); if (dmod := get(dmode!*,'dname)) then onoff(dmod,nil)>> where !*msg := nil; begin scalar dmode!*,!*exp,!*gcd,!*keepsqrts,!*limitedfactors,!*mcd, !*rationalize,!*structure,!*uncached,kord!*, ans1,denexp,badbit,nexp,oneterm,!*precise; !*keepsqrts := !*limitedfactors := t; % !*sqrt := t; !*exp := !*gcd := !*mcd := !*structure := !*uncached := t; dmode!* := nil; if !*algint then << % The algint code now needs precise off. % !*precise := t; % Start a clean slate (in terms of SQRTSAVE) for this % integral. sqrt!-intvar:=!*q2f simpsqrti variable; if (red sqrt!-intvar) or (lc sqrt!-intvar neq 1) or (ldeg sqrt!-intvar neq 1) then interr "Sqrt(x) not properly formed" else sqrt!-intvar:=mvar sqrt!-intvar; basic!-listofallsqrts:=listofallsqrts; basic!-listofnewsqrts:=listofnewsqrts; sqrtsave(basic!-listofallsqrts,basic!-listofnewsqrts, list(variable . variable))>>; coefft := (1 ./ 1); % Collect simple coefficients. expression := int!-simp car u; if varchange then <<depend1(car varchange,cdr varchange,t); expression := int!-subsq(expression,{varchange})>>; denexp := 1 ./ denr expression; % Get into two bits nexp := numr expression; while not atom nexp and null cdr nexp and not depends(mvar nexp,variable) do <<coefft := multsq(coefft,(((caar nexp) . 1) . nil) ./ 1); nexp := lc nexp>>; ans1 := nil; while nexp do begin % Collect by zvariables scalar x,zv,tmp; if atom nexp then << x:=!*f2q nexp; nexp:=nil >> else << x:=!*t2q car nexp; nexp:=cdr nexp >>; x := multsq(x,denexp); zv := findzvars(getvariables x,list variable,variable,nil); begin scalar oldzlist; while oldzlist neq zv do << oldzlist := zv; foreach zz in oldzlist do zv:=findzvars(distexp(pseudodiff(zz,variable)), zv,variable,t)>>; % The following line was added to make, for example, % int(df(sin(x)/x),x) return the expected result. zv := sort(zv, function ordp) end; tmp := ans1; while tmp do <<if zv=caar tmp then <<rplacd(car tmp,addsq(cdar tmp,x)); tmp := nil; zv := nil>> else tmp := cdr tmp>>; if zv then ans1 := (zv . x) . ans1 end; if length ans1 = 1 then oneterm := t; % Efficiency nexp := ans1; ans := nil ./ 1; badbit:=nil ./ 1; % SQ zero while nexp do % Run down the terms <<u := cdar nexp; if !*trdint then <<princ "Integrate"; printsq u; princ "with Zvars "; print caar nexp>>; ans1 := errorset!*(list('integratesq,mkquote u, mkquote variable,mkquote loglist, mkquote caar nexp), !*backtrace); nexp := cdr nexp; if errorp ans1 then badbit := addsq(badbit,u) else <<ans := addsq(caar ans1, ans); badbit:=addsq(cdar ans1,badbit)>>>>; if !*trdint then <<prin2 "Partial answer="; printsq ans; prin2 "To do="; printsq badbit>>; % We have run down the terms. If there are any bad bits, redo % them. However, since a non-zero badbit implies that % integratesq aborted, the internal variable order may be % confused. So we reset kord!* and reorder expressions in this % case. if badbit neq '(nil . 1) then <<setkorder nil; badbit := reordsq badbit; ans := reordsq ans; coefft := reordsq coefft; if !*trdint then <<princ "Retrying..."; printsq badbit>>; if oneterm and ans = '(nil . 1) then ans1 := nil else ans1 := errorset!*(list('integratesq,mkquote badbit, mkquote variable,mkquote loglist,nil), !*backtrace); if null ans1 or errorp ans1 then ans := addsq(ans,simpint1(badbit . variable . w)) else <<ans := addsq(ans,caar ans1); %% FJW: It is possible for ans here to be just a %% spurious constant term, in which case we discard it. if not smemq(variable, ans) then ans := nil ./ 1; %% This may not be the best place for this fix, but I %% don't see how it can ever do any harm. [I don't %% think we need a full depend test here.] if cdar ans1 neq '(nil . 1) then ans := addsq(ans, simpint1(cdar ans1 . variable . w)) >>>>; end; ans := multsq(coefft,ans); %Put back coefficient, preserving order. % if errorp ans % then return <<put('sqrt,'simpfn,sqrtfn); % if !*failhard then error1(); % simpint1(expression . variable . w)>> % else ans := car ans; % expression := sqrtchk numr ans ./ sqrtchk denr ans; if !*trdint then << printc "Resimp and all that"; printsq ans >>; % We now need to check that all simplifications have been done % but we have to make sure INT is not resimplified, and that SIMP % does not complain at getting the same argument again. put('int,'simpfn,'simpiden); put('sqrt,'simpfn,sqrtfn); << if dmod then onoff(dmod,t); % added by Alan Barnes if cflag then onoff('complex,t)>> where !*msg := nil; oldvarstack := varstack!*; varstack!* := nil; % ans := errorset!*(list('resimp,mkquote ans),t); ans := errorset!*(list('int!-resub,mkquote ans,mkquote varchange),t); put('int,'simpfn,'simpint); varstack!* := oldvarstack; return if errorp ans then error1() else car ans end; symbolic procedure int!-resub(x,v); % {sq,alist} -> sq % Undo any variable change and resimplify. if v then <<x := int!-subsq(x,{revpr v}); depend1(car v,cdr v,nil); resimp x>> else resimp x; symbolic procedure int!-subsq(x,v); % {sq,alist} -> sq % A version of subsq with the int and df operators unprotected. % Intended for straightforward change of variable names only. begin scalar subfuncs,subfg!*; subfuncs := {remprop('df,'subfunc),remprop('int,'subfunc)}; x := subsq(x,v); put('df,'subfunc,car subfuncs); put('int,'subfunc,cadr subfuncs); return x end; symbolic procedure numlistp u; % True if u is a list of numbers. null u or numberp car u and numlistp cdr u; % symbolic procedure sqrtchk u; % % U is a standard form. Result is another standard form with square % % roots replaced by half powers. % if domainp u then u % else if not eqcar(mvar u,'sqrt) % then addf(multpf(lpow u,sqrtchk lc u),sqrtchk red u) % % else if mvar u = '(sqrt -1) % % then addf(multpf(mksp('i,ldeg u),sqrtchk lc u),sqrtchk red u) % else addf(multpf(mksp(list('expt,cadr mvar u,'(quotient 1 2)), % ldeg u), % sqrtchk lc u), % sqrtchk red u); symbolic procedure int!-simp u; % Converts U to canonical form, including the resimplification of % *sq forms. subs2 resimp simp!* u; put('int,'simpfn,'simpint); symbolic procedure integratesq(integrand,var,xlogs,zv); begin scalar varlist,x,zlist,!*noncomp; if !*trint then << printc "Start of Integration; integrand is "; printsq integrand >>; !*noncomp := noncomfp numr integrand or noncomfp denr integrand; varlist:=getvariables integrand; varlist:=varsinlist(xlogs,varlist); %in case more exist in xlogs if zv then zlist := zv else <<zlist := findzvars(varlist,list var,var,nil); % Important kernels. % The next section causes problems with nested exponentials or logs. begin scalar oldzlist; while oldzlist neq zlist do <<oldzlist := zlist; foreach zz in oldzlist do zlist := findzvars(distexp(pseudodiff(zz,var)), zlist,var,t)>> end>>; if !*trint then << printc "Determination of the differential field descriptor"; printc "gives the functions:"; print zlist >>; %% Look for rational powers in the descriptor %% If there is make a suitable transformation and do the sub integral %% and return the revised integral x := look_for_substitute(integrand, var, zlist); if x then return x; %% End of rational patch if !*purerisch and not allowedfns zlist then return (nil ./ 1) . integrand; % If it is not suitable for Risch. varlist := setdiff(varlist,zlist); % varlist := purge(zlist,varlist); % Now zlist is list of things that depend on x, and varlist is list % of constant kernels in integrand. if !*algint and cdr zlist and algfnpl(zlist,var) then return algebraiccase(integrand,zlist,varlist) else return transcendentalcase(integrand,var,xlogs,zlist,varlist) end; symbolic procedure distexp(l); if null l then nil else if atom car l then car l . distexp cdr l else if (caar l = 'expt) and (cadar l = 'e) then begin scalar ll; ll:=caddr car l; if eqcar(ll,'plus) then << ll:=foreach x in cdr ll collect list('expt,'e,x); return ('times . ll) . distexp cdr l >> else return car l . distexp cdr l end else distexp car l . distexp cdr l; symbolic procedure pseudodiff(a,var); if atom a then % **** Treat diffs correctly?? if depends(a,var) then list prepsq simpdf(list(a,var)) else nil else if car a memq '(atan equal log plus quotient sqrt times minus) then begin scalar aa,bb; foreach zz in cdr a do << bb:=pseudodiff(zz,var); aa:= union(bb,aa) >>; return aa end else if car a eq 'expt then if depends(cadr a,var) then if depends(caddr a,var) then prepsq simp list('log,cadr a) . %% a(x)^b(x) cadr a . caddr a . union(pseudodiff(cadr a,var),pseudodiff(caddr a,var)) else cadr a . pseudodiff(cadr a,var) %% a(x)^b else caddr a . pseudodiff(caddr a,var) %% a^b(x) else list prepsq simpdf(list(a,var)); symbolic procedure look_for_substitute(integrand, var, zz); % Search for rational power transformations begin scalar res; if atom zz then return nil else if (res := look_for_rational(integrand, var, zz)) then return res else if (res := look_for_quad(integrand, var, zz)) then return res else if (res := look_for_substitute(integrand, var, car zz)) then return res else return look_for_substitute(integrand, var, cdr zz) end; symbolic procedure look_for_rational(integrand, var, zz); % Look for a form x^(n/m) in the field descriptor, and transform % the integral if it is found. Note that the sqrt form may be used % as well as exponentials. Return nil if no transformation if (car zz = 'sqrt and cadr zz = var) then look_for_rational1(integrand, var, 2) else if (car zz = 'expt) and (cadr zz = var) and (listp caddr zz) and (caaddr zz = 'quotient) and (numberp cadr caddr zz) and (numberp caddr caddr zz) then look_for_rational1(integrand, var, caddr caddr zz) else nil; symbolic procedure look_for_rational1(integrand, var, m); % Actually do the transformation and integral begin scalar newvar, res, ss, mn2m!-1; newvar := gensym(); mn2m!-1 := !*f2q(((newvar .** (m-1)) .* m) .+ nil); %% print ("Integrand was " . integrand); % x => y^m, and dx => m y^(m-1) integrand := multsq(subsq(integrand, list(var . list('expt,newvar,m))), mn2m!-1); if !*trint then << prin2 "Integrand is transformed to "; printsq integrand >>; begin scalar intvar; intvar := newvar; % To circumvent an algint bug. res := integratesq(integrand, newvar, nil, nil); end; ss := list(newvar . list('expt,var, list('quotient, 1, m))); res := subsq(car res, ss) . subsq(quotsq(cdr res, mn2m!-1), ss); if !*trint then << printc "Transforming back..."; printsq car res; prin2 " plus a bad part of "; printsq cdr res >>; return res end; symbolic procedure look_for_quad(integrand, var, zz); % Look for a form sqrt(a+bx+cx^2) in the field descriptor % and transform to the appropriate asin, acosh or asinh. % Return nil if no transformation found if !*algint then nil % as Algint does it better?? else begin if (car zz = 'sqrt and listp cadr zz and caadr zz = 'plus) or (car zz = 'expt and listp cadr zz and caadr zz = 'plus and listp caddr zz and car caddr zz = 'quotient and fixp caddr caddr zz) then << zz := simp cadr zz; if (cdr zz = 1) then << zz := cdr coeff1(prepsq zz, var, nil); if length zz = 2 then return begin % Linear scalar a, b; scalar nvar, res, ss; a := car zz; b := cadr zz; if (depends(a,var) or depends(b,var)) then return nil; nvar := gensym(); if !*trint then << prin2 "Linear shift suggested "; prin2 a; prin2 " "; prin2 b; terpri(); >>; integrand := subsq(integrand, % Make the substitution list(var . list('quotient, list('difference, list('expt,nvar,2),a), b))); integrand := multsq(integrand, % and the dx component simp list('quotient,list('times,nvar,2), b)); % integrand := subsq(integrand, % list(var . list('difference, nvar, a))); % integrand := multsq(integrand, simp b); if !*trint then << prin2 "Integrand is transformed by substitution to "; printsq integrand; prin2 "using substitution "; prin2 var; prin2 " -> "; printsq simp list('quotient, list('difference,list('expt,nvar,2),a), b); >>; res := integratesq(integrand, nvar, nil, nil); ss := list(nvar . list('sqrt,list('plus,list('times,var,b), a))); res := subsq(car res, ss) . subsq(multsq(cdr res, simp list('quotient,b, list('times,nvar,2))), ss); %% Should one reject if there is a bad bit?? return res; end else if length zz = 3 then return begin % A quadratic scalar a, b, c; a := car zz; b := cadr zz; c:= caddr zz; if (depends(a,var) or depends(b,var) or depends(c,var)) then return nil; a := simp list('difference, a, % Re-centre list('times,b,b, list('quotient,1,list('times,4,c)))); if null numr a then return nil; % Power occurred. b := simp list('quotient, b, list('times, 2, c)); c := simp c; return if minusf numr c then << if minusf numr a then begin scalar !*hyperbolic; !*hyperbolic := t; return look_for_invhyp(integrand,nil,var,a,b,c) end else look_for_asin(integrand,var,a,b,c)>> else << if minusf numr a then look_for_invhyp(integrand,t,var,a,b,c) else look_for_invhyp(integrand,nil,var,a,b,c) >> end else if length zz = 5 then return begin % A quartic scalar a, b, c, d, e, nn, dd, mm; a := car zz; b := cadr zz; c:= caddr zz; d := cadddr zz; e := car cddddr zz; if not(b = 0) or not(d = 0) then return nil; if (depends(a,var) or depends(c,var)) or depends(e,var) then return nil; nn := numr integrand; dd := denr integrand; if denr(mm :=quotsq(nn ./ 1, !*kk2q var)) = 1 and even_power(numr mm, var) and even_power(dd, var) then << % substitute x -> sqrt(y) return sqrt_substitute(numr mm, dd, var); >>; if denr(mm :=quotsq(dd ./ 1, !*kk2q var)) = 1 and even_power(nn, var) and even_power(numr mm, var) then << % substitute x -> sqrt(y) return sqrt_substitute(nn, multf(dd,!*kk2f var), var); >>; return nil; end; >>>>; return nil; end; symbolic procedure look_for_asin(integrand, var, a, b, c); % Actually do the transformation and integral begin scalar newvar, res, ss, sqmn, onemth, m, n; m := prepsq a; n := prepsq c; b := prepsq b; newvar := gensym(); sqmn := prepsq apply1(get('sqrt, 'simpfn), list list('quotient, list('minus,n), m)); onemth := list('cos, newvar); ss := list('sin, newvar); powlis!* := list(ss, 2, '(nil . t), list('difference,1,list('expt,onemth,2)), nil) . powlis!*; integrand := subs2q multsq(subsq(integrand, list(var . list('difference, list('quotient,ss,sqmn), b))), quotsq(onemth := simp onemth, simp sqmn)); if !*trint then << prin2 "Integrand is transformed by substitution to "; printsq integrand; prin2 "using substitution "; prin2 var; prin2 " -> "; printsq simp list('difference, list('quotient, ss, sqmn), b); >>; res := integratesq(integrand, newvar, nil, nil); powlis!* := cdr powlis!*; ss:= list(newvar . list('asin,list('times,list('plus,var,b),sqmn))); res := subsq(car res, ss) . subsq(quotsq(cdr res, onemth), ss); if !*trint then << printc "Transforming back..."; printsq car res; prin2 " plus a bad part of "; printsq cdr res >>; if (car res = '(nil . 1)) then return nil; return res; end; symbolic procedure look_for_invhyp(integrand, do_acosh, var, a, b, c); % Actually do the transformation and integral; uses acosh/asinh form % depending on second argument begin scalar newvar, res, ss, sqmn, onemth, m, n, realdom; m := prepsq a; n := prepsq c; b := prepsq b; newvar := gensym(); if do_acosh then << sqmn := prepsq apply1(get('sqrt, 'simpfn), list list('quotient, n, list('minus, m))); onemth := list('sinh, newvar); ss := list('cosh, newvar) >> else << sqmn:= prepsq apply1(get('sqrt,'simpfn),list list('quotient,n,m)); onemth := list('cosh, newvar); ss := list('sinh, newvar) >>; powlis!* := list(ss, 2, '(nil . t), list((if do_acosh then 'plus else 'difference), list('expt, onemth, 2),1), nil) . powlis!*; % print ("sqmn" . sqmn); print("onemth" . onemth); print ("ss" . ss); % print cdddar powlis!*; integrand := subs2q multsq(subsq(integrand, list(var . list('difference,list('quotient,ss,sqmn),b))), quotsq(onemth := simp onemth, simp sqmn)); if !*trint then << prin2 "Integrand is transformed by substitution to "; printsq integrand; prin2 "using substitution "; prin2 var; prin2 " -> "; printsq simp list('difference, list('quotient, ss, sqmn), b); >>; realdom := not smember('(sqrt -1),integrand); % print integrand; print realdom; res := integratesq(integrand, newvar, nil, nil); powlis!* := cdr powlis!*; if !*hyperbolic then << ss := list(if do_acosh then 'acosh else 'asinh, list('times,list('plus,var,b), sqmn)); >> else << ss := list('times,list('plus,var,b), sqmn); ss := if do_acosh then subst(ss,'ss, '(log (plus ss (sqrt (difference (times ss ss) 1))))) else subst(ss,'ss,'(log (plus ss (sqrt (plus (times ss ss) 1))))) >>; ss := list(newvar . ss); res := sqrt2top subsq(car res, ss) . sqrt2top subsq(quotsq(cdr res, onemth), ss); if !*trint then << printc "Transforming back..."; printsq car res; prin2 " plus a bad part of "; printsq cdr res >>; if (car res = '(nil . 1)) then return nil; if realdom and smember('(sqrt -1),res) then << if !*trint then print "Wrong sheet"; return nil; % Wrong sheet? >>; return res end; 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>>; return if (v := opmtch u) then simp v else symint u % FJW: symbolic integral end; mkop 'int!*; put('int!*,'simpfn,'simpint!*); symbolic procedure simpint!* u; begin scalar x; return if (x := opmtch('int . u)) then simp x else simpiden('int!* . u) end; symbolic procedure remakesf u; %remakes standard form U, substituting operator INT for INT!*; if domainp u then u else addf(multpf(if eqcar(mvar u,'int!*) then mksp('int . cdr mvar u,ldeg u) else lpow u,remakesf lc u), remakesf red u); symbolic procedure allowedfns u; if null u then t else if atom car u then (car u=intvar) or not depends(car u,intvar) else if (caar u = 'expt and not (cadar u = 'e) and not depends(cadar u, intvar) and depends(caddar u, intvar)) then nil else if flagp(caar u,'transcendental) then allowedfns cdr u else nil; symbolic procedure look_for_power(integrand, var); begin scalar nn, dd, mm; nn := numr integrand; dd := denr integrand; if denr(mm :=quotsq(nn ./ 1, !*kk2q var)) = 1 and even_power(numr mm, var) and even_power(dd, var) then << % substitute x -> sqrt(y) return sqrt_substitute(numr mm, dd, var); >>; if denr(mm :=quotsq(dd ./ 1, !*kk2q var)) = 1 and even_power(nn, var) and even_power(numr mm, var) then << % substitute x -> sqrt(y) return sqrt_substitute(nn, numr mm, var); >>; return nil; end; symbolic procedure even_power(xpr, var); if atom xpr then t else if mvar xpr = var then << if evenp pdeg lpow xpr then even_power(lc xpr, var) and even_power(red xpr, var) else nil >> else if eqcar(mvar xpr, 'expt) and cadr mvar xpr = var and evenp caddr mvar xpr then t else if atom mvar xpr then even_power(lc xpr, var) and even_power(red xpr, var) else if even_power(red xpr, var) and even_power(lc xpr, var) then even_prep(mvar xpr, var); symbolic procedure even_prep(xpr,var); if xpr = var then nil else if atom xpr then t else if eqcar(xpr, 'expt) and cadr xpr = var and evenp caddr xpr then t else if even_prep(car xpr, var) then even_prep(cdr xpr, var); symbolic procedure sqrt_substitute(nn, dd, var); begin scalar newvar, integrand, res, ss, !*keepsqrts; newvar := gensym(); integrand := subst(list('sqrt,newvar), var, list('quotient, prepsq (nn ./ dd), 2)); integrand := prepsq simp integrand; integrand := simp integrand; begin scalar intvar; intvar := newvar; % To circumvent an algint bug/oddity res := integratesq(integrand, newvar, nil, nil); end; ss := list(newvar . list('expt, var, 2)); res := subsq(car res, ss) . multsq((((var .^ 1) .* 2) .+ nil) ./ 1, subsq(cdr res, ss)); if !*trint then << printc "Transforming back..."; printsq car res; prin2 " plus a bad part of "; printsq cdr res >>; return res end; % The following rules probably belong in other places. %----------------------------------------------------------------------- algebraic; operator ci,si; % ei. % FJW: ci,si also defined in specfn(sfint.red), so ... symbolic((algebraic operator ci,si) where !*msg=nil); intrules := {e^(~n*acosh(~x)) => (sqrt(x^2-1)+x)^n when numberp n, e^(~n*asinh(~x)) => (sqrt(x^2+1)+x)^n when numberp n, e^(acosh(~x)) => (sqrt(x^2-1)+x), e^(asinh(~x)) => (sqrt(x^2+1)+x), cosh(log(~x)) => (x^2+1)/(2*x), sinh(log(~x)) => (x^2-1)/(2*x), % These next two are rather uncertain. int(log(~x)/(~b-x),x) => dilog(x/b), int(log(~x)/(~b*x-x^2),x) => dilog(x/b)/b + log(x)^2/(2b), %% FJW: Next 2 rules replaced by ~~ rules below %% int(e^(~x^2),x) => erf(i*x)*sqrt(pi)/(2i), %% int(1/e^(~x^2),x) => erf(x) * sqrt(pi)/2, %% FJW: Missing sqrt(b): %% int(e^(~b*~x^2),x) => erf(i*x)*sqrt(pi)/(2i*sqrt(b)), int(e^(~~b*~x^2),x) => erf(i*sqrt(b)*x)*sqrt(pi)/(2i*sqrt(b)), %% FJW: Rule missing: int(e^(~x^2/~b),x) => erf(i*x/sqrt(b))*sqrt(pi)*sqrt(b)/(2i), %% FJW: Missing sqrt(b): %% int(1/e^(~b*~x^2),x) => erf(x)*sqrt(pi)/(2sqrt(b)), int(1/e^(~~b*~x^2),x) => erf(sqrt(b)*x)*sqrt(pi)/(2sqrt(b)), %% FJW: Rule missing: int(1/e^(~x^2/~b),x) => erf(x/sqrt(b))*sqrt(pi)*sqrt(b)/2, df(ei(~x),x) => exp(x)/x, int(e^(~~b*~x)/x,x) => ei(b*x), % FJW int(e^(~x/~b)/x,x) => ei(x/b), int(1/(exp(~x*~~b)*x),x) => ei(-x*b), % FJW int(1/(exp(~x/~b)*x),x) => ei(-x/b), %% FJW: Next 2 rules replaced by ~~ rules above %% int(e^~x/x,x) => ei(x), %% int(1/(e^~x*x),x) => ei(-x), int(~a^~x/x,x) => ei(x*log(a)), int(1/((~a^~x)*x),x) => ei(-x*log(a)), df(si(~x),x) => sin(x)/x, int(sin(~~b*~x)/x,x) => si(b*x), % FJW int(sin(~x/~b)/x,x) => si(x/b), % FJW %% int(sin(~x)/x,x) => si(x), % FJW int(sin(~x)/x^2,x) => -sin(x)/x +ci(x), int(sin(~x)^2/x,x) =>(log(x)-ci(2x))/2, df(ci(~x),x) => cos(x)/x, int(cos(~~b*~x)/x,x) => ci(b*x), % FJW int(cos(~x/~b)/x,x) => ci(x/b), % FJW %% int(cos(~x)/x,x) => ci(x), % FJW int(cos(~x)/x^2,x) => -cos(x)/x -si(x), int(cos(~x)^2/x,x) =>(log(x)+ci(2x)/2), int(1/log(~~b*~x),x) => ei(log(b*x))/b, % FJW int(1/log(~x/~b),x) => ei(log(x/b))*b, % FJW %% int(1/log(~x),x) => ei(log(x)), % FJW %% int(1/log(~x+~b),x) => ei(log(x+b)), % FJW int(1/log(~~a*~x+~b),x) => ei(log(a*x+b))/b, % FJW int(1/log(~x/~a+~b),x) => ei(log(x/a+b))/b, % FJW int(~x/log(~x),x) => ei(2*log(x)), int(~x^~n/log(x),x) => ei((n+1)*log(x)) when fixp n, int(1/(~x^~n*log(x)),x) => ei((-n+1)*log(x)) when fixp n}; let intrules; endmodule; end;