Artifact 27821706b0723b9ed3b0675334cc7a80ab6df4034ebaca0799a818cbd8273912:
- Executable file
r37/packages/int/trcase.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: 12295) [annotate] [blame] [check-ins using] [more...]
module trcase; % Driving routine for integration of transcendental fns. % Authors: Mary Ann Moore and Arthur C. Norman. % Modifications by: John P. Fitch. fluid '(!*backtrace !*failhard !*nowarnings !*purerisch !*reverse !*trint badpart ccount cmap cmatrix content cval denbad denominator!* indexlist lhs!* loglist lorder orderofelim power!-list!* rhs!* sillieslist sqfr sqrtflag sqrtlist tanlist svar varlist xlogs zlist); % !*reverse: flag to re-order zlist. % !*nowarnings: flag to lose messages. global '(!*number!* !*ratintspecial !*seplogs !*spsize!* !*statistics gensymcount); switch failhard; exports transcendentalcase; imports backsubst4cs,countz,createcmap,createindices,df2q,dfnumr, difflogs,fsdf,factorlistlist,findsqrts,findtrialdivs,gcdf,mkvect, interr,logstosq,mergin,multbyarbpowers,!*multf, % multsqfree, printdf,printsq,quotf,rationalintegrate,putv, simpint1,solve!-for!-u,sqfree,sqmerge,sqrt2top,substinulist,trialdiv, mergein,negsq,addsq,f2df,mknill,pnth,invsq,multsq,domainp,mk!*sq, mksp,prettyprint; % Note that SEPLOGS keeps logarithmic part of result together as a % kernel form, but this can lead to quite messy results. symbolic procedure transcendentalcase(integrand,svar,xlogs,zlist,varlist); begin scalar divlist,jhd!-content,content,prim,sqfr,dfu,indexlist, % JHD!-CONTENT is local, while CONTENT is free (set in SQFREE). sillieslist,originalorder,wrongway,power!-list!*, sqrtlist,tanlist,loglist,dflogs,eprim,dfun,unintegrand, sqrtflag,badpart,rhs!*,lhs!*,gcdq,cmap,cval,orderofelim,cmatrix; scalar ccount,denominator!*,result,denbad; gensymcount:=0; integrand:=sqrt2top integrand; % Move the sqrts to the numerator. if !*trint then << printc "Extension variables z<i> are"; print zlist>>; % if !*ratintspecial and null cdr zlist then % return rationalintegrate(integrand,svar); % *** now unnormalize integrand, maybe ***. begin scalar w,gg; gg:=1; foreach z in zlist do <<w := subs2 diffsq(simp z,svar); % subs2q? gg := !*multf(gg,quotf(denr w,gcdf(denr w,gg)))>>; gg := quotf(gg,gcdf(gg,denr integrand)); unintegrand := (!*multf(gg,numr integrand) ./ !*multf(gg,denr integrand)); % multf? if !*trint then << printc "After unnormalization the integrand is "; printsq unintegrand >> end; divlist := findtrialdivs zlist; % Also puts some things on loglist sometimes. sqrtlist := findsqrts zlist; divlist := trialdiv(denr unintegrand,divlist); % N.B. the next line also sets 'content' as a free variable. % Since SQFREE may be used later, we copy it into JHD!-CONTENT. prim := sqfree(cdr divlist,zlist); jhd!-content := content; printfactors(prim,nil); eprim := sqmerge(countz car divlist,prim,nil); printfactors(eprim,t); % if !*trint then <<terpri(); % printsf denominator!*; % terpri(); % printc "...content is:"; % printsf jhd!-content>>; sqfr := for each u in eprim collect multup u; % sqfr := multsqfree eprim; % if !*trint then <<printc "...sqfr is:"; % prettyprint sqfr>>; if !*reverse then zlist := reverse zlist; % Alter order function. indexlist := createindices zlist; % if !*trint then << printc "...indices are:"; % prettyprint indexlist>>; dfu:=dfnumr(svar,car divlist); % if !*trint then << terpri(); % printc "************ Derivative of u is:"; % printsq dfu>>; loglist := append(loglist,factorlistlist prim); %%% nconc? loglist := mergein(xlogs,loglist); loglist := mergein(tanlist,loglist); cmap := createcmap(); ccount := length cmap; if !*trint then <<printc "Loglist "; print loglist>>; dflogs := difflogs(loglist,denr unintegrand,svar); if !*trint then <<printc "************ 'Derivative' of logs is:"; printsq dflogs>>; dflogs := addsq((numr unintegrand) ./ 1,negsq dflogs); % Put everything in reduction eqn over common denominator. gcdq := gcdf(denr dflogs,denr dfu); dfun := !*multf(numr dfu,denbad:=quotf(denr dflogs,gcdq)); denbad := !*multf(denr dfu,denbad); denbad := !*multf(denr unintegrand,denbad); dflogs := !*multf(numr dflogs,quotf(denr dfu,gcdq)); dfu := dfun; % Now DFU and DFLOGS are S.F.s. rhs!* := multbyarbpowers f2df dfu; if checkdffail(rhs!*,svar) then <<if !*trint then printsq checkdffail(rhs!*,svar); interr "Simplification fails on above expression">>; if !*trint then << printc "Distributed Form of Numerator is:"; printdf rhs!*>>; lhs!* := f2df dflogs; % if checkdffail(lhs!*,svar) then interr "Simplification failure"; if !*trint then << printc "Distributed Form of integrand is:"; printdf lhs!*; terpri()>>; cval := mkvect(ccount); for i := 0:ccount do putv(cval,i,nil ./ 1); power!-list!* := tansfrom(rhs!*,zlist,indexlist,0); lorder:=maxorder(power!-list!*,zlist,0); originalorder := for each x in lorder collect x; % Must copy as it is overwritten. if !*trint then << printc "Maximum order for variables determined as "; print lorder >>; if !*statistics then << !*number!*:=0; !*spsize!*:=1; foreach xx in lorder do !*spsize!*:=!*spsize!* * (xx+1) >>; % That calculates the largest U that can appear. dfun:=solve!-for!-u(rhs!*,lhs!*,nil); backsubst4cs(nil,orderofelim,cmatrix); % if !*trint then if ccount neq 0 then printvecsq cval; if !*statistics then << prin2 !*number!*; prin2 " used out of "; printc !*spsize!* >>; badpart:=substinulist badpart; %substitute for c<i> still in badpart. dfun:=df2q substinulist dfun; result:= !*multsq(dfun,!*invsq(denominator!* ./ 1)); result:= !*multsq(result,!*invsq(jhd!-content ./ 1)); dflogs:=logstosq(); if not null numr dflogs then << if !*seplogs and (not domainp numr result) then << result:=mk!*sq result; result:=(mksp(result,1) .* 1) .+ nil; result:=result ./ 1 >>; result:=addsq(result,dflogs)>>; if !*trint then << %% prettyprint result; terpri(); printc "*****************************************************"; printc "************ THE INTEGRAL IS : **********************"; printc "*****************************************************"; terpri(); printsq result; terpri()>>; if badpart then begin scalar n,oorder; if !*trint then printc "plus a part which has not been integrated"; lhs!*:=badpart; lorder:=maxorder(power!-list!*,zlist,0); oorder:=originalorder; n:=length lorder; while lorder do << if car lorder > car originalorder then wrongway:=t; if car lorder=car originalorder then n:= n-1; lorder:=cdr lorder; originalorder:=cdr originalorder >>; %% if n=0 then wrongway:=t; % Nothing changed if !*trint and wrongway then printc "Went wrong way"; dfun:=df2q badpart; %% if !*trint %% then <<printsq dfun; printc "Denbad = "; printsf denbad>>; dfun:= !*multsq(dfun,invsq(denbad ./ 1)); badpart := dfun; %%% ***** if wrongway then << if !*trint then printc "Resetting...."; result:=nil ./ 1; dfun := integrand; badpart:=dfun >>; if rootcheckp(unintegrand,svar) then return simpint1(integrand . svar.nil) . (nil ./ 1) else if !*purerisch or allowedfns zlist then << badpart := dfun; dfun := nil ./ 1 >> % JPff else << !*purerisch:=t; if !*trint then <<printc " Applying transformations ..."; printsq dfun>>; denbad:=transform(dfun,svar); if denbad=dfun then << dfun:=nil ./ 1; badpart:=denbad >> else <<denbad:=errorset!*(list('integratesq,mkquote denbad, mkquote svar,mkquote xlogs, nil), !*backtrace); %%% JPF fix for distributive version if not atom denbad then << denbad:=car denbad; % As from an errorset dfun:=untan car denbad; if (dfun neq '(nil . 1)) then badpart:=untan cdr denbad; % There is still a chance that we went the wrong way. % Decide that it is better if there is no bad part if car badpart and not(badpart=denbad) then << wrongway:=nil; lhs!*:=f2df car badpart; lorder:=maxorder(power!-list!*,zlist,0); n:=length lorder; while lorder do << if car lorder > car oorder then wrongway:=t; if car lorder=car oorder then n:= n-1; lorder:=cdr lorder; oorder:=cdr oorder >>; if wrongway or (n=0) then << if !*trint then printc "Still backwards"; dfun := nil ./ 1; badpart := integrand>>>>>> else <<badpart := dfun; dfun := nil ./ 1 >>>>>>; if !*failhard then rerror(int,9,"FAILHARD switch set"); if !*seplogs and not domainp result then << result:=mk!*sq result; if not numberp result then result:=(mksp(result,1) .* 1) .+ nil; result:=result ./ 1>>; result:=addsq(result,dfun) end else badpart:=nil ./ 1; return (sqrt2top result . badpart) % JPff end; symbolic procedure checkdffail(u,v); % Sometimes simplification fails and this gives the integrator the % idea that something is a constant when it is not. Check for this. if null u then nil else if depends(lc u,v) then lc u else checkdffail(red u,v); symbolic procedure printfactors(w,prdenom); % W is a list of factors to each power. If PRDENOM is true % this prints denominator of answer, else prints square-free % decomposition. begin scalar i,wx; i:=1; if prdenom then << denominator!* := 1; if !*trint then printc "Denominator of 1st part of answer is:"; if not null w then w:=cdr w >>; loopx: if w=nil then return; if !*trint then <<prin2 "Factors of multiplicity "; prin2 i; prin2t ":"; terpri()>>; wx:=car w; while not null wx do << if !*trint then printsf car wx; for j:=1 : i do denominator!*:= !*multf(car wx,denominator!*); wx:=cdr wx >>; i:=i+1; w:=cdr w; go to loopx end; % unfluid '(dfun svar xlogs); endmodule; end;