Artifact f57c141c52933b918c1a3090606a19ab6657e5cf151ba605df00613ee4724da8:
- File
r33/int.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: 154866) [annotate] [blame] [check-ins using] [more...]
module int!-intro; % General support for REDUCE integrator. % Authors: A. C. Norman and P. M. A. Moore. % Modified by: J. Davenport, J. P. Fitch, A. C. Hearn. % Note that at one point, INT had been flagged SIMP0FN. However, that % lead to problems when the arguments of INT contained pattern % variables. fluid '(!*conscount !*noextend !*pvar gaussiani); global '(btrlevel frlis!* gensymcount initl!*); !*conscount:=10000; % default maximum number of conses in certain % operations. !*pvar:='!_a; btrlevel := 5; %default to a reasonably full backtrace. % The next smacro is needed at this point to define gaussiani. symbolic smacro procedure !*kk2f u; !*p2f mksp(u,1); gaussiani := !*kk2f '(sqrt -1); gensymcount := 0; initl!* := append('(!*noextend), initl!*); flag('(interr),'transfer); %For the compiler; flag ('(atan dilog erf expint expt log tan),'transcendental); comment Kludge to define derivative of an integral and integral of a derivative; frlis!* := union('(!=x !=y),frlis!*); put('df,'opmtch,'(((int !=y !=x) !=x) (nil . t) (evl!* !=y) nil) . get('df,'opmtch)); put('int,'opmtch,'(((df !=y !=x) !=x) (nil . t) (evl!* !=y) nil) . get('int,'opmtch)); put('evl!*,'opmtch,'(((!=x) (nil . t) !=x nil))); put('evl!*,'simpfn,'simpiden); %Various functions used throughout the integrator. smacro procedure !*kk2q a; ((mksp(a,1) .* 1) .+ nil) ./ 1; symbolic smacro procedure divsf(u,v); sqrt2top(u ./ v); symbolic procedure flatten u; if null u then nil else if atom u then list u else if atom car u then car u . flatten cdr u else nconc(flatten car u,flatten cdr u); symbolic procedure int!-gensym1 u; << gensymcount:=gensymcount+1; compress append(explode u,explode gensymcount) >>; symbolic smacro procedure maninp(u,v,w); interr "MANINP called -- not implemented"; symbolic procedure mknill n; if n=0 then nil else nil . mknill(n-1); % Various selectors written as macros. smacro procedure argof u; % Argument of a unary function. cadr u; smacro procedure firstsubs u; % The first substitution in a substitution list. car u; smacro procedure lsubs u; car u; smacro procedure rsubs u; cdr u; smacro procedure lfirstsubs u; caar u; smacro procedure rfirstsubs u; cdar u; put('nthroot,'simpfn,'simpiden); % The binary n-th root operator nthroot(x,2)=sqrt(x) % no simplification is used here. % Hope is that pbuild introduces it, and simplog removes it. % Selectors for the taylor series structure. % Format is: %function.((first.last computed so far) . assoc list of computed terms). % ***store-hack-1***: % remove this macro if more store is available. smacro procedure tayshorten u;nil; smacro procedure taylordefn u; car u; symbolic smacro procedure taylorfunction u; caar u; smacro procedure taylornumbers u; cadr u; smacro procedure taylorfirst u; caadr u; smacro procedure taylorlast u; cdadr u; smacro procedure taylorlist u; cddr u; smacro procedure taylormake(fn,nums,alist); fn.(nums.alist); endmodule; module contents; % Authors: Mary Ann Moore and Arthur C. Norman fluid '(clogflag content indexlist sqfr varlist zlist); exports contents,contentsmv,dfnumr,difflogs,factorlistlist,multsqfree, multup,sqfree,sqmerge; imports int!-fac,fquotf,gcdf,interr,!*multf,partialdiff,quotf,ordop, addf,negf,domainp,difff,mksp,negsq,invsq,addsq,!*multsq,diffsq; comment we assume no power substitution is necessary in this module; symbolic procedure contents(p,v); % Find the contents of the polynomial p wrt variable v; % Note that v may not be the main variable of p; if domainp(p) then p else if v=mvar p then contentsmv(p,v,nil) else if ordop(v,mvar p) then p else contentsmv(makemainvar(p,v),v,nil); symbolic procedure contentsmv(p,v,sofar); % Find contents of polynomial P; % V is main variable of P; % SOFAR is partial result; if sofar=1 then 1 else if domainp p then gcdf(p,sofar) else if not v=mvar p then gcdf(p,sofar) else contentsmv(red p,v,gcdf(lc p,sofar)); symbolic procedure makemainvar(p,v); % Bring v up to be the main variable in polynomial p; % Note that the reconstructed p must be used with care since; % it does not conform to the normal reduce ordering rules; if domainp p then p else if v=mvar p then p else mergeadd(mulcoeffsby(makemainvar(lc p,v),lpow p,v), makemainvar(red p,v),v); symbolic procedure mulcoeffsby(p,pow,v); % Multiply each coefficient in p by the standard power pow; if null p then nil else if domainp p or not v=mvar p then ((pow .* p) .+ nil) else (lpow p .* ((pow .* lc p) .+ nil)) .+ mulcoeffsby(red p,pow,v); symbolic procedure mergeadd(a,b,v); % Add polynomials a and b given that they have same main variable v; if domainp a or not v=mvar a then if domainp b or not v=mvar b then addf(a,b) else lt b .+ mergeadd(a,red b,v) else if domainp b or not v=mvar b then lt a .+ mergeadd(red a,b,v) else (lambda xc; if xc=0 then (lpow a .* addf(lc a,lc b)) .+ mergeadd(red a,red b,v) else if xc>0 then lt a .+ mergeadd(red a,b,v) else lt b .+ mergeadd(a,red b,v)) (tdeg lt a-tdeg lt b); symbolic procedure sqfree(p,vl); if (null vl) or (domainp p) then <<content:=p; nil>> else begin scalar w,v,dp,gg,pg,dpg,p1,w1; w:=contents(p,car vl); % content of p ; p:=quotf(p,w); % make p primitive; w:=sqfree(w,cdr vl); % process content by recursion; if p=1 then return w; v:=car vl; % pick out variable from list; while not (p=1) do << dp:=partialdiff(p,v); gg:=gcdf(p,dp); pg:=quotf(p,gg); dpg:=negf partialdiff(pg,v); p1:=gcdf(pg,addf(quotf(dp,gg),dpg)); w1:=p1.w1; p:=gg>>; return sqmerge(reverse w1,w,t) end; symbolic procedure sqmerge(w1,w,simplew1); % w and w1 are lists of factors of each power. if simplew1 is true % then w1 contains only single factors for each power. ; if null w1 then w else if null w then if car w1=1 then nil.sqmerge(cdr w1,w,simplew1) else (if simplew1 then list car w1 else car w1). sqmerge(cdr w1,w,simplew1) else if car w1=1 then (car w).sqmerge(cdr w1,cdr w,simplew1) else append(if simplew1 then list car w1 else car w1,car w). sqmerge(cdr w1,cdr w,simplew1); symbolic procedure multup l; % l is a list of s.f.'s. result is s.f. for product of elements of l; begin scalar res; res:=1; while not null l do << res:=multf(res,car l); l:=cdr l >>; return res end; symbolic procedure diflist(l,cl,x,rl); % Differentiates l (list of s.f.'s) wrt x to produce the sum of % terms for the derivative of numr of 1st part of answer. cl is % coefficient list (s.f.'s) & rl is list of derivatives we have % dealt with so far. Result is s.q.; if null l then nil ./ 1 else begin scalar temp; temp:=!*multf(multup rl,multup cdr l); temp:=!*multsq(difff(car l,x),!*f2q temp); temp:=!*multsq(temp,(car cl) ./ 1); return addsq(temp,diflist(cdr l,cdr cl,x,(car l).rl)) end; symbolic procedure multsqfree w; % W is list of sqfree factors. result is product of each LIST IN W % to give one polynomial for each sqfree power; if null w then nil else (multup car w).multsqfree cdr w; symbolic procedure l2lsf l; % L is a list of kernels. result is a list of same members as s.f.'s; if null l then nil else ((mksp(car l,1) .* 1) .+ nil).l2lsf cdr l; symbolic procedure dfnumr(x,dl); % Gives the derivative of the numr of the 1st part of answer. % dl is list of any exponential or 1+tan**2 that occur in integrand % denr. these are divided out from result before handing it back. % result is s.q., ready for printing. begin scalar temp1,temp2,coeflist,qlist,count; if not null sqfr then << count:=0; qlist:=cdr sqfr; coeflist:=nil; while not null qlist do << count:=count+1; coeflist:=count.coeflist; qlist:=cdr qlist >>; coeflist:=reverse coeflist >>; temp1:=!*multsq(diflist(l2lsf zlist,l2lsf indexlist,x,nil), !*f2q multup sqfr); if not null sqfr and not null cdr sqfr then << temp2:=!*multsq(diflist(cdr sqfr,coeflist,x,nil), !*f2q multup l2lsf zlist); temp2:=!*multsq(temp2,(car sqfr) ./ 1) >> else temp2:=nil ./ 1; temp1:=addsq(temp1,negsq temp2); temp2:=cdr temp1; temp1:=car temp1; qlist:=nil; while not null dl do << if not car dl member qlist then qlist:=(car dl).qlist; dl:=cdr dl >>; while not null qlist do << temp1:=quotf(temp1,car qlist); qlist:=cdr qlist >>; return temp1 ./ temp2 end; symbolic procedure difflogs(ll,denm1,x); % LL is list of log terms (with coeffts), den is common denominator % over which they are to be put. Result is s.q. for derivative of all % these wrt x. if null ll then nil ./ 1 else begin scalar temp,qu,cvar,logoratan,arg; logoratan:=caar ll; cvar:=cadar ll; arg:=cddar ll; temp:=!*multsq(cvar ./ 1,diffsq(arg,x)); if logoratan='iden then qu:=1 ./ 1 else if logoratan='log then qu:=arg else if logoratan='atan then qu:=addsq(1 ./ 1,!*multsq(arg,arg)) else interr "Logoratan=? in difflogs"; %Note call to special division routine; qu:=fquotf(!*multf(!*multf(denm1,numr temp), denr qu),numr qu); %*MUST* GO EXACTLY; temp:=!*multsq(!*invsq (denr temp ./ 1),qu); %result of fquotf is a s.q; return !*addsq(temp,difflogs(cdr ll,denm1,x)) end; symbolic procedure factorlistlist (w,clogflag); % W is list of lists of sqfree factors in s.f. result is list of log % terms required for integral answer. the arguments for each log fn % are in s.q.; begin scalar res,x,y; while not null w do << x:=car w; while not null x do << y:=facbypp(car x,varlist); while not null y do << res:=append(int!-fac car y,res); y:=cdr y >>; x:=cdr x >>; w:=cdr w >>; return res end; symbolic procedure facbypp(p,vl); % Use contents/primitive parts to try to factor p. if null vl then list p else begin scalar princilap!-part,co; co:=contents(p,car vl); vl:=cdr vl; if co=1 then return facbypp(p,vl); %this var no help. princilap!-part:=quotf(p,co); %primitive part. if princilap!-part=1 then return facbypp(p,vl); %again no help; return nconc(facbypp(princilap!-part,vl),facbypp(co,vl)) end; endmodule; module csolve; % routines to do with the C constants. % Author: John P. Fitch. fluid '(ccount cmap cmatrix cval loglist neweqn); global '(!*trint); exports backsubst4cs,createcmap,findpivot,printspreadc,printvecsq, spreadc,subst4eliminateds; imports nth,interr,!*multf,printsf,printsq,quotf,putv,negf,invsq, negsq,addsq,multsq,mksp,addf,domainp,pnth; symbolic procedure findpivot cvec; % Finds first non-zero element in CVEC and returns its cell number. % If no such element exists, result is nil. begin scalar i,x; i:=1; x:=getv(cvec,i); while i<ccount and null x do << i:=i+1; x:=getv(cvec,i) >>; if null x then return nil; return i end; symbolic procedure subst4eliminatedcs(neweqn,substorder,ceqns); % Substitutes into NEWEQN for all the C's that have been eliminated so % far. These are given by CEQNS. SUBSTORDER gives the order of % substitution as well as the constant multipliers. Result is the % transformed NEWEQN. if null substorder then neweqn else begin scalar nxt,row,cvar,temp; row:=car ceqns; nxt:=car substorder; if null (cvar:=getv(neweqn,nxt)) then return subst4eliminatedcs(neweqn,cdr substorder,cdr ceqns); nxt:=getv(row,nxt); for i:=0 : ccount do << temp:=!*multf(nxt,getv(neweqn,i)); temp:=addf(temp,negf !*multf(cvar,getv(row,i))); putv(neweqn,i,temp) >>; return subst4eliminatedcs(neweqn,cdr substorder,cdr ceqns) end; symbolic procedure backsubst4cs(cs2subst,cs2solve,cmatrix); % Solves the C-eqns and sets vector CVAL to the C-constant values % CMATRIX is a list of matrix rows for C-eqns after Gaussian % elimination has been performed. CS2SOLVE is a list of the remaining % C's to evaluate and CS2SUBST are the C's we have evaluated already. if null cmatrix then nil else begin scalar eqnn,cvar,already,substlist,temp,temp2; eqnn:=car cmatrix; cvar:=car cs2solve; already:=nil ./ 1; % The S.Q. nil. substlist:=cs2subst; % Now substitute for previously evaluated c's: while not null substlist do << temp:=car substlist; if not null getv(eqnn,temp) then already:=addsq(already,multsq(getv(eqnn,temp) ./ 1, getv(cval,temp))); substlist:=cdr substlist >>; % Now solve for the c given by cvar (any remaining c's assumed zero). temp:=negsq addsq(getv(eqnn,0) ./ 1,already); if not null (temp2:=quotf(numr temp,getv(eqnn,cvar))) then temp:=temp2 ./ denr temp else temp:=multsq(temp,invsq(getv(eqnn,cvar) ./ 1)); if not null numr temp then putv(cval,cvar, resimp rootextractsq subs2q temp); backsubst4cs(reversewoc(cvar . reversewoc cs2subst), cdr cs2solve,cdr cmatrix) end; %********************************************************************** % Routines to deal with linear equations for the constants C. %********************************************************************** symbolic procedure createcmap; %Sets LOGLIST to list of things of form (LOG C-constant f), where f is % function linear in one of the z-variables and C-constant is in S.F. % When creating these C-constant names, the CMAP is also set up and % returned as the result. begin scalar i,l,c; l:=loglist; i:=1; while not null l do << c:=(int!-gensym1('c) . i) . c; i:=i+1; rplacd(car l,((mksp(caar c,1) .* 1) .+ nil) . cdar l); l:=cdr l >>; if !*trint then printc ("Constants Map" . c); return c end; symbolic procedure spreadc(eqnn,cvec1,w); % Sets a vector 'cvec1' to coefficients of c<i> in eqnn. if domainp eqnn then putv(cvec1,0,addf(getv(cvec1,0), !*multf(eqnn,w))) else begin scalar mv,t1,t2; spreadc(red eqnn,cvec1,w); mv:=mvar eqnn; t1:=assoc(mv,cmap); %tests if it is a c var. if not null t1 then return << t1:=cdr t1; %loc in vector for this c. if not (tdeg lt eqnn=1) then interr "Not linear in c eqn"; t2:=addf(getv(cvec1,t1),!*multf(w,lc eqnn)); putv(cvec1,t1,t2) >>; t1:=((lpow eqnn) .* 1) .+ nil; %this main var as sf. spreadc(lc eqnn,cvec1,!*multf(w,t1)) end; symbolic procedure printspreadc cvec1; begin for i:=0 : ccount do << prin2 i; printc ":"; printsf(getv(cvec1,i)) >>; printc "End of printspreadc output" end; %symbolic procedure printvecsq cvec; %% Print contents of cvec which contains s.q.'s (not s.f.'s). %% Starts from cell 1 not 0 as above routine (printspreadc). % begin % for i:=1 : ccount do << % prin2 i; % printc ":"; % if null getv(cvec,i) then printc "0" % else printsq(getv(cvec,i)) >>; % printc "End of printvecsq output" % end; endmodule; module cuberoot; % Cube roots of standard forms. % Authors: Mary Ann Moore and Arthur C. Norman. fluid '(cuberootflag); exports cuberootdf; imports contentsmv,gcdf,!*multf,nrootn,partialdiff,printdf,quotf,vp2, mksp,mk!*sq,domainp; symbolic procedure cuberootsq a; cuberootf numr a ./ cuberootf denr a; symbolic procedure cuberootf p; begin scalar ip,qp; if null p then return nil; ip:=cuberootf1 p; qp:=cdr ip; ip:=car ip; %respectable and nasty parts of the cuberoot. if numberp qp and onep qp then return ip; %exact root found. qp:=list('expt,prepf qp,'(quotient 1 3)); cuberootflag:=t; %symbolic cube-root introduced. qp:=(mksp(qp,1).* 1) .+ nil; return !*multf(ip,qp) end; symbolic procedure cuberootf1 p; % Returns a . b with p=a**2*b. % Does this need power reduction? if domainp p then nrootn(p,3) else begin scalar co,ppp,g,pg; co:=contentsmv(p,mvar p,nil); %contents of p. ppp:=quotf(p,co); %primitive part. % now consider ppp=p1*p2**2*p3**3*p4**4*... co:=cuberootf1(co); %process contents via recursion. g:=gcdf(ppp,partialdiff(ppp,mvar ppp)); % g=p2*p3**2*p4**3*... if not domainp g then << pg:=quotf(ppp,g); %pg=p1*p2*p3*p4*... g:=gcdf(g,partialdiff(g,mvar g)); % g=g3*g4**2*g5**3*... g:=gcdf(g,pg)>>; %a triple factor of ppp. if domainp g then pg:=1 . ppp else << pg:=quotf(ppp,!*multf(g,!*multf(g,g))); %what's left. pg:=cuberootf1 pg; %split that up. rplaca(pg,!*multf(car pg,g))>>; %put in the thing found here. rplaca(pg,!*multf(car pg,car co)); rplacd(pg,!*multf(cdr pg,cdr co)); return pg end; endmodule; module idepend; % Routines for considering dependency among variables. % Authors: Mary Ann Moore and Arthur C. Norman. fluid '(taylorvariable); exports dependspl,dependsp,involvesq,involvsf; imports taylorp,domainp; symbolic procedure dependsp(x,v); if null v then t else if depends(x,v) then x else if atom x then if x eq v then x else nil else if car x = '!*sq then involvesq(cadr x,v) else if taylorp x then if v eq taylorvariable then taylorvariable else nil else begin scalar w; if x=v then return v; % Check if a prefix form expression depends on the variable v. % Note this assumes the form x is in normal prefix notation; w := x; % preserve the dependency; x := cdr x; % ready to recursively check arguments; scan: if null x then return nil; % no dependency found; if dependsp(car x,v) then return w; x:=cdr x; go to scan end; symbolic procedure involvesq(sq,term); involvesf(numr sq,term) or involvesf(denr sq,term); symbolic procedure involvesf(sf,term); if domainp sf or null sf then nil else dependsp(mvar sf,term) or involvesf(lc sf,term) or involvesf(red sf,term); symbolic procedure dependspl(dep!-list,var); % True if any member of deplist (a list of prefix forms) depends on % var. dep!-list and (dependsp(car dep!-list,var) or dependspl(cdr dep!-list,var)); symbolic procedure taylorp exxpr; % Sees if a random entity is a taylor expression. not atom exxpr and not atom car exxpr and flagp(taylorfunction exxpr,'taylor); endmodule; module df2q; % Conversion from distributive to standard forms. % Authors: Mary Ann Moore and Arthur C. Norman. fluid '(indexlist zlist); exports df2q; imports addf,gcdf,mksp,!*multf,quotf; comment We assume that results already have reduced powers, so that no power substitution is necessary; symbolic procedure df2q p; % Converts distributed form P to standard quotient; begin scalar n,d,gg,w; if null p then return nil ./ 1; d:=denr lc p; w:=red p; while not null w do << gg:=gcdf(d,denr lc w); %get denominator of answer... d:=!*multf(d,quotf(denr lc w,gg)); %..as lcm of denoms in input w:=red w >>; n:=nil; %place to build numerator of answer while not null p do << n:=addf(n,!*multf(xl2f(lpow p,zlist,indexlist), !*multf(numr lc p,quotf(d,denr lc p)))); p:=red p >>; return n ./ d end; symbolic procedure xl2f(l,z,il); % L is an exponent list from a D.F., Z is the Z-list, % IL is the list of indices. % Value is L converted to standard form. ; if null z then 1 else if car l=0 then xl2f(cdr l,cdr z,cdr il) else if not atom car l then begin scalar temp; if caar l=0 then temp:= car il else temp:=list('plus,car il,caar l); temp:=mksp(list('expt,car z,temp),1); return !*multf(((temp .* 1) .+ nil), xl2f(cdr l,cdr z,cdr il)) end % else if minusp car l then ; % multsq(invsq (((mksp(car z,-car l) .* 1) .+ nil)), ; % xl2f(cdr l,cdr z,cdr il)) ; else !*multf((mksp(car z,car l) .* 1) .+ nil, xl2f(cdr l,cdr z,cdr il)); endmodule; module distrib; % Routines for manipulating distributed forms. % Authors: Mary Ann Moore and Arthur C. Norman. fluid '(indexlist sqrtlist zlist); exports dfprintform,multbyarbpowers,negdf,quotdfconst,sub1ind,vp1, vp2,plusdf,multdf,multdfconst,orddf; imports interr,addsq,negsq,exptsq,simp,domainp,mk!*sq,addf, multsq,invsq,minusp,mksp,sub1; %*********************************************************************** % NOTE: The expressions lt,red,lc,lpow have been used on distributed % forms as the latter's structure is sufficiently similar to % s.f.'s. However lc df is a s.q. not a s.f. and lpow df is a % list of the exponents of the variables. This also makes % lt df different. Red df is d.f. as expected. %**********************************************************************; symbolic procedure plusdf(u,v); % U and V are D.F.'s. Value is D.F. for U+V; if null u then v else if null v then u else if lpow u=lpow v then (lambda(x,y); if null numr x then y else (lpow u .* x) .+ y) (!*addsq(lc u,lc v),plusdf(red u,red v)) else if orddf(lpow u,lpow v) then lt u .+ plusdf(red u,v) else (lt v) .+ plusdf(u,red v); symbolic procedure orddf(u,v); % U and V are the LPOW of a D.F. - i.e. the list of exponents ; % Value is true if LPOW U '>' LPOW V and false otherwise ; if null u then if null v then interr "Orddf = case" else interr "Orddf v longer than u" else if null v then interr "Orddf u longer than v" else if exptcompare(car u,car v) then t else if exptcompare(car v,car u) then nil else orddf(cdr u,cdr v); symbolic procedure exptcompare(x,y); if atom x then if atom y then x>y else nil else if atom y then t else car x > car y; symbolic procedure negdf u; if null u then nil else (lpow u .* negsq lc u) .+ negdf red u; symbolic procedure multdf(u,v); % U and V are D.F.'s. Value is D.F. for U*V; % reduces squares of square-roots as it goes; if null u or null v then nil else begin scalar y; %use (a+b)*(c+d) = (a*c) + a*(c+d) + b*(c+d); y:=multerm(lt u,lt v); %leading terms; y:=plusdf(y,multdf(red u,v)); y:=plusdf(y,multdf((lt u) .+ nil,red v)); return y end; symbolic procedure multerm(u,v); %multiply two terms to give a D.F.; begin scalar coef; coef:=!*multsq(cdr u,cdr v); %coefficient part; return multdfconst(coef,mulpower(car u,car v)) end; symbolic procedure mulpower(u,v); % u and v are exponent lists. multiply corresponding forms; begin scalar r,s; r:=addexptsdf(u,v); if not null sqrtlist then s:=reduceroots(r,zlist); r:=(r .* (1 ./ 1)) .+ nil; if not (s=nil) then r:=multdf(r,s); return r end; symbolic procedure reduceroots(r,zl); begin scalar s; while not null r do << if eqcar(car zl,'sqrt) then s:=tryreduction(r,car zl,s); r:=cdr r; zl:=cdr zl >>; return s end; symbolic procedure tryreduction(r,var,s); begin scalar x; x:=car r; %current exponent if not atom x then << r:=x; x:=car r >>; %numeric part if (x=0) or (x=1) then return s; %no reduction possible x:=divide(x,2); rplaca(r,cdr x); %reduce exponent as redorded x:=car x; var:=simp cadr var; %sqrt arg as a s q var:=!*exptsq(var,x); x:=multdfconst(1 ./ denr var,f2df numr var); %distribute if s=nil then s:=x else s:=multdf(s,x); return s end; symbolic procedure addexptsdf(x,y); % X and Y are LPOW's of D.F. Value is list of sum of exponents; if null x then if null y then nil else interr "X too long" else if null y then interr "Y too long" else exptplus(car x,car y).addexptsdf(cdr x,cdr y); symbolic procedure exptplus(x,y); if atom x then if atom y then x+y else list (x+car y) else if atom y then list (car x +y) else interr "Bad exponent sum"; symbolic procedure multdfconst(x,u); % X is S.Q. not involving Z variables of D.F. U. Value is D.F.; % for X*U; if (null u) or (null numr x) then nil else lpow u .* !*multsq(x,lc u) .+ multdfconst(x,red u); %symbolic procedure quotdfconst(x,u); % multdfconst(!*invsq x,u); symbolic procedure f2df p; % P is standard form. Value is P in D.F.; if domainp p then dfconst(p ./ 1) else if mvar p member zlist then plusdf(multdf(vp2df(mvar p,tdeg lt p,zlist),f2df lc p), f2df red p) else plusdf(multdfconst(((lpow p .* 1) .+ nil) ./ 1,f2df lc p), f2df red p); symbolic procedure vp1(var,degg,z); % Takes VAR and finds it in Z (=list), raises it to power DEGG and puts % the result in exponent list form for use in a distributed form. if null z then interr "Var not in z-list after all" else if var=car z then degg.vp2 cdr z else 0 . vp1(var,degg,cdr z); symbolic procedure vp2 z; % Makes exponent list of zeroes. if null z then nil else 0 . vp2 cdr z; symbolic procedure vp2df(var,exprn,z); % Makes VAR**EXPRN into exponent list and then converts the resulting % power into a distributed form. % Special care with square-roots. if eqcar(var,'sqrt) and (exprn>1) then mulpower(vp1(var,exprn,z),vp2 z) else (vp1(var,exprn,z) .* (1 ./ 1)) .+ nil; symbolic procedure dfconst q; % Makes a distributed form from standard quotient constant Q; if numr q=nil then nil else ((vp2 zlist) .* q) .+ nil; %df2q moved to a section of its own. symbolic procedure df2printform p; %Convert to a standard form good enough for printing. if null p then nil else begin scalar mv,co; mv:=xl2q(lpow p,zlist,indexlist); if mv=(1 ./ 1) then << co:=lc p; if denr co=1 then return addf(numr co, df2printform red p); co:=mksp(mk!*sq co,1); return (co .* 1) .+ df2printform red p >>; co:=lc p; if not (denr co=1) then mv:=!*multsq(mv,1 ./ denr co); mv:=mksp(mk!*sq mv,1) .* numr co; return mv .+ df2printform red p end; symbolic procedure xl2q(l,z,il); % L is an exponent list from a D.F.,Z is the Z-list, IL is the list of % indices. Value is L converted to standard quotient. ; if null z then 1 ./ 1 else if car l=0 then xl2q(cdr l,cdr z,cdr il) else if not atom car l then begin scalar temp; if caar l=0 then temp:= car il else temp:=list('plus,car il,caar l); temp:=mksp(list('expt,car z,temp),1); return !*multsq(((temp .* 1) .+ nil) ./ 1, xl2q(cdr l,cdr z,cdr il)) end else if minusp car l then !*multsq(!*invsq(((mksp(car z,-car l) .* 1) .+ nil) ./ 1), xl2q(cdr l,cdr z,cdr il)) else !*multsq(((mksp(car z,car l) .* 1) .+ nil) ./ 1, xl2q(cdr l,cdr z,cdr il)); %symbolic procedure sub1ind power; % if atom power then power-1 % else list sub1 car power; symbolic procedure multbyarbpowers u; % Multiplies the ordinary D.F., U, by arbitrary powers % of the z-variables; % i-1 j-1 k-1 % i.e. x z z ... so result is D.F. with the exponent list % 1 2 %appropriately altered to contain list elements instead of numeric ones. if null u then nil else ((addarbexptsdf lpow u) .* lc u) .+ multbyarbpowers red u; symbolic procedure addarbexptsdf x; % Adds the arbitrary powers to powers in exponent list, X, to produce % new exponent list. e.g. 3 -> (2) to represent x**3 now becoming : % 3 i-1 i+2 ; % x * x = x . ; if null x then nil else list exptplus(car x,-1) . addarbexptsdf cdr x; endmodule; module divide; % Exact division of standard forms to give a S Q. % Authors: Mary Ann Moore and Arthur C. Norman. fluid '(residue sqrtlist zlist); global '(!*trdiv !*trint); exports fquotf,testdivdf,dfquotdf; imports df2q,f2df,gcdf,interr,multdf,negdf,plusdf,printdf,printsf, quotf,multsq,invsq,negsq; % Intended for dividing out known factors as produced by the % integration program. horrible and slow, i expect!! symbolic procedure dfquotdf(a,b); begin scalar residue; if (!*trint or !*trdiv) then << printc "Dfquotdf called on "; printdf a; printdf b>>; a:=dfquotdf1(a,b); if (!*trint or !*trdiv) then << printc "Quotient given as "; printdf a >>; if not null residue then begin scalar gres,w; if !*trint or !*trdiv then << printc "Residue in dfquotdf ="; printdf residue; printc "Which should be zero"; w:=residue; gres:=numr lc w; w:=red w; while not null w do << gres:=gcdf(gres,numr lc w); w:=red w >>; printc "I.e. the following vanishes"; printsf gres>>; interr "Non-exact division due to a log term" end; return a end; symbolic procedure fquotf(a,b); % Input: a and b standard quotients with (a/b) an exact % division with respect to the variables in zlist, % but not necessarily obviously so. the 'non-obvious' problems % will be because of (e.g.) square-root symbols in b % output: standard quotient for (a/b) % (prints message if remainder is not 'clearly' zero. % A must not be zero. begin scalar t1; if null a then interr "A=0 in fquotf"; t1:=quotf(a,b); %try it the easy way if not null t1 then return t1 ./ 1; %ok return df2q dfquotdf(f2df a,f2df b) end; symbolic procedure dfquotdf1(a,b); begin scalar q; if null b then interr "Attempt to divide by zero"; q:=sqrtlist; %remove sqrts from denominator, maybe. while not null q do begin scalar conj; conj:=conjsqrt(b,car q); %conjugate wrt given sqrt if not (b=conj) then << a:=multdf(a,conj); b:=multdf(b,conj) >>; q:=cdr q end; q:=dfquotdf2(a,b); residue:=reversewoc residue; return q end; symbolic procedure dfquotdf2(a,b); % As above but a and b are distributed forms, as is the result. if null a then nil else begin scalar xd,lcd; xd:=xpdiff(lpow a,lpow b); if xd='failed then << xd:=lt a; a:=red a; residue:=xd .+ residue; return dfquotdf2(a,b) >>; lcd:= !*multsq(lc a,!*invsq lc b); if null numr lcd then return dfquotdf2(red a,b); % Should not be necessary; lcd := xd .* lcd; xd:=plusdf(a,multdf(negdf (lcd .+ nil),b)); if xd and (lpow xd = lpow a % Again, should not be necessary; or xpdiff(lpow xd,lpow b) = 'failed) then <<if !*trint or !*trdiv then <<printc "Dfquotdf trouble:"; printdf xd>>; xd := rootextractdf xd; if !*trint or !*trdiv then printdf xd>>; return lcd .+ dfquotdf2(xd,b) end; symbolic procedure rootextractdf u; if null u then nil else begin scalar v; v := resimp rootextractsq lc u; return if null numr v then rootextractdf red u else (lpow u .* v) .+ rootextractdf red u end; symbolic procedure rootextractsq u; if null numr u then u else rootextractf numr u ./ rootextractf denr u; symbolic procedure rootextractf v; if domainp v then v else begin scalar u,r,c,x,p; u := mvar v; p := ldeg v; r := rootextractf red v; c := rootextractf lc v; if null c then return r else if atom u then return (lpow v .* c) .+ r else if car u eq 'sqrt or car u eq 'expt and eqcar(caddr u,'quotient) and car cdaddr u = 1 and numberp cadr cdaddr u then <<p := divide(p,if car u eq 'sqrt then 2 else cadr cdaddr u); if car p = 0 then return if null c then r else (lpow v .* c) .+ r else if numberp cadr u then <<c := multd(cadr u ** car p,c); p := cdr p>> else <<x := simpexpt list(cadr u,car p); if denr x = 1 then <<c := multf(numr x,c); p := cdr p>>>>>>; return if p=0 then addf(c,r) else if null c then r else ((u to p) .* c) .+ r end; % The following hack makes sure that the results of differentiation % gets passed through ROOTEXTRACT % a) This should not be done this way, since the effect is global % b) Should this be done via TIDYSQRT? put('df,'simpfn,'simpdf!*); symbolic procedure simpdf!* u; begin scalar v,v1; v:=simpdf u; v1:=rootextractsq v; if not(v1=v) then return resimp v1 else return v end; symbolic procedure xpdiff(a,b); %Result is list a-b, or 'failed' if a member of this would be negative. if null a then if null b then nil else interr "B too long in xpdiff" else if null b then interr "A too long in xpdiff" else if car b>car a then 'failed else (lambda r; if r='failed then 'failed else (car a-car b) . r) (xpdiff(cdr a,cdr b)); symbolic procedure conjsqrt(b,var); % Subst(var=-var,b). if null b then nil else conjterm(lpow b,lc b,var) .+ conjsqrt(red b,var); symbolic procedure conjterm(xl,coef,var); % Ditto but working on a term. if involvesp(xl,var,zlist) then xl .* negsq coef else xl .* coef; symbolic procedure involvesp(xl,var,zl); % Check if exponent list has non-zero power for variable. if null xl then interr "Var not found in involvesp" else if car zl=var then (not zerop car xl) else involvesp(cdr xl,var,cdr zl); endmodule; module driver; % Driving routines for integration program. % Author: Mary Ann Moore and Arthur C. Norman. % Modifications by: John P. Fitch. fluid '(!*backtrace !*exp !*gcd !*keepsqrts !*mcd !*nolnr !*purerisch !*rationalize !*sqrt !*structure !*uncached basic!-listofnewsqrts basic!-listofallsqrts expression gaussiani intvar listofnewsqrts listofallsqrts loglist sqrt!-intvar sqrt!-places!-alist variable varlist xlogs zlist); global '(!*algint !*failhard !*trint); exports integratesq,simpint,purge,simpint1; imports algebraiccase,algfnpl,findzvars,getvariables,interr,printsq, transcendentalcase,varsinlist,kernp,simpcar,prepsq,mksq,simp, opmtch,formlnr; switch algint,nolnr,trint; % 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; begin scalar ans,expression,variable,loglist,w, !*purerisch,intvar,listofnewsqrts,listofallsqrts, sqrtfn,sqrt!-intvar,sqrt!-places!-alist, basic!-listofallsqrts,basic!-listofnewsqrts; if atom u or null cdr u then rederr "Not enough arguments for INT"; variable := !*a2k cadr u; w := cddr u; if w then rederr "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. begin scalar dmode!*,!*exp,!*gcd,!*keepsqrts,!*mcd,!*sqrt, !*rationalize,!*structure,!*uncached; !*keepsqrts := !*sqrt := t; !*exp := !*gcd := !*mcd := !*structure := !*uncached := t; dmode!* := nil; if !*algint then <<intvar:=variable; % until fix JHD code % 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))>>; expression := int!-simp car u; % loglist := for each x in w collect int!-simp x; ans := errorset('(integratesq expression variable loglist), !*backtrace,!*backtrace); end; 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; % We now need to check that all simplifications have been done % but we have to make sure INT is not resimplified. put('int,'simpfn,'simpiden); ans := errorset('(resimp expression),t,!*backtrace); put('int,'simpfn,'simpint); put('sqrt,'simpfn,sqrtfn); return if errorp ans then error1() else car ans end; 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 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); begin scalar varlist,zlist; if !*trint then << printc "Integrand is..."; printsq integrand >>; varlist:=getvariables integrand; varlist:=varsinlist(xlogs,varlist); %in case more exist in xlogs 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 "with 'new' functions :"; print zlist >>; if !*purerisch and not allowedfns zlist then return simpint1 (integrand . var.nil); % If it is not suitable for Risch; 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 nil else if car a memq '(atan equal log plus quotient sqrt times) then begin scalar aa,bb; foreach zz in cdr a do << bb:=pseudodiff(zz,var); if aa then aa:=bb . aa else bb >>; return aa end else if car a eq 'expt then if depends(cadr a,var) then prepsq simp list('log,cadr a) . cadr a . caddr a . append(pseudodiff(cadr a,var),pseudodiff(caddr a,var)) else caddr a . pseudodiff(caddr a,var) else list prepsq simpdf(list(a,var)); symbolic procedure simpint1 u; begin scalar v,!*sqrt; 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; u:=errorset(list('simp,mkquote v), !*backtrace,!*backtrace); if pairp u then v:=car u; !*nolnr:= cdr !*nolnr; return v>>; return if (v := opmtch u) then simp v else if !*failhard then rederr "FAILHARD switch set" else mksq(u,1) 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 or flagp(caar u,'transcendental) then allowedfns cdr u else nil; symbolic procedure purge(a,b); if null a then b else if null b then nil else purge(cdr a,delete(car a,b)); endmodule; module d3d4; % Splitting of cubics and quartics. % Authors: Mary Ann Moore and Arthur C. Norman. fluid '(knowndiscrimsign zlist); global '(!*trint); exports cubic,quartic; imports covecdf,cuberootf,nth,forceazero,makepolydf,multdf,multdfconst, !*multf,negdf,plusdf,printdf,printsf,quadratic,sqrtf,vp1,vp2,addf, negf; symbolic procedure cubic(pol,var,res); %Split the univariate (wrt z-vars) cubic pol, at least if a %change of origin puts it in the form (x-a)**3-b=0; begin scalar a,b,c,d,v,shift,p; v:=covecdf(pol,var,3); shift:=forceazero(v,3); %make coeff x**2 vanish. %also checks univariate. % if shift='failed then go to prime; a:=getv(v,3); b:=getv(v,2); %=0, I hope!; c:=getv(v,1); d:=getv(v,0); if !*trint then << printc "Cubic has coefficients"; printsf a; printsf b; printsf c; printsf d >>; if not null c then << if !*trint then printc "Cubic too hard to split"; go to exit >>; a:=cuberootf(a); %can't ever fail; d:=cuberootf(d); if !*trint then << printc "Cube roots of a and d are"; printsf a; printsf d>>; %now a*(x+shift)+d is a factor of pol; %create x+shift in p; p:=(vp2 zlist .* shift) .+ nil; p:=(vp1(var,1,zlist) .* (1 ./ 1)) .+ p; %(x+shift); b:=nil; b:=(vp2 zlist .* (d ./ 1)) .+ b; b:=plusdf(b,multdfconst(a ./ 1,p)); b:=makepolydf b; %get rid of denominator. if !*trint then << printc "One factor of the cubic is"; printdf b >>; res:=('log . b) . res; %now form the (quadratic) cofactor; b:=(vp2 zlist .* (!*multf(d,d) ./ 1)) .+ nil; b:=plusdf(b,multdfconst(negf !*multf(a,d) ./ 1,p)); b:=plusdf(b,multdfconst(!*multf(a,a) ./ 1, multdf(p,p))); return quadratic(makepolydf b,var,res); %deal with what is left; prime: if !*trint then printc "The following cubic does not split"; exit: if !*trint then printdf pol; return ('log . pol) . res end; symbolic procedure quartic(pol,var,res); %Splits univariate (wrt z-vars) quartics that can be written %in the form (x-a)**4+b*(x-a)**2+c; begin scalar a,b,c,d,ee,v,shift,p,q,p1,p2,dsc; v:=covecdf(pol,var,4); shift:=forceazero(v,4); %make coeff x**3 vanish; % if shift='failed then go to prime; a:=getv(v,4); b:=getv(v,3); %=0, I hope. c:=getv(v,2); d:=getv(v,1); ee:=getv(v,0); if !*trint then << printc "Quartic has coefficients"; printsf a; printsf b; printsf c; printsf d; printsf ee >>; if d then <<if !*trint then printc "Quartic too hard to split"; go to exit >>; b:=c; c:=ee; %squash up the notation; if knowndiscrimsign eq 'negative then go to complex; dsc := addf(!*multf(b,b),multf(-4,!*multf(a,c))); p2 := minusf c; if not p2 and minusf dsc then go to complex; p1 := null b or minusf b; if not p1 then if p2 then p1 := t else p2 := t; p1 := if p1 then 'positive else 'negative; p2 := if p2 then 'negative else 'positive; a := sqrtf a; dsc := sqrtf dsc; if a eq 'failed or dsc eq 'failed then go to prime; ee := invsq(addf(a,a) ./ 1); d := multsq(addf(b,negf dsc) ./ 1,ee); ee := multsq(addf(b,dsc) ./ 1,ee); if !*trint then <<printc "Quadratic factors will have coefficients"; printsf a; print 0; printsq d; printc "or"; printsq ee>>; p := (vp2 zlist .* shift) .+ nil; p := (vp1(var,1,zlist) .* (1 ./ 1)) .+ p; %(x+shift); q := multdf(p,p); %square of same; q := multdfconst(a ./ 1,q); p := plusdf(q,(vp2 zlist .* d) .+ nil); q := plusdf(q,(vp2 zlist .* ee) .+ nil); if !*trint then <<printc "Allowing for change of origin:"; printdf p; printdf q>>; knowndiscrimsign := p1; res := quadratic(p,var,res); knowndiscrimsign := p2; res := quadratic(q,var,res); go to quarticdone; complex: a:=sqrtf(a); c:=sqrtf(c); if a eq 'failed or c eq 'failed then go to prime; b:=addf(!*multf(2,!*multf(a,c)),negf b); b:=sqrtf b; if b eq 'failed then go to prime; %now a*(x+shift)**2 (+/-) b*(x+shift) + c is a factor. if !*trint then << printc "Quadratic factors will have coefficients"; printsf a; printsf b; printsf c>>; p:=(vp2 zlist .* shift) .+ nil; p:=(vp1(var,1,zlist) .* (1 ./ 1)) .+ p; %(x+shift); q:=multdf(p,p); %square of same; p:=multdfconst(b ./ 1,p); q:=multdfconst(a ./ 1,q); q:=plusdf(q,(vp2 zlist .* (c ./ 1)) .+ nil); if !*trint then << printc "Allowing for change of origin, p (+/-) q with p,q="; printdf p; printdf q>>; %now p+q and p-q are the factors of the quartic; knowndiscrimsign := 'negative; res:=quadratic(plusdf(q,p),var,res); res:=quadratic(plusdf(q,negdf p),var,res); quarticdone: knowndiscrimsign := nil; if !*trint then printc "Quartic done"; return res; prime: if !*trint then printc "The following quartic does not split"; exit: if !*trint then printdf pol; return ('log . pol) . res end; endmodule; module factr; % Crude factorization routine for integrator. % Authors: Mary Ann Moore and Arthur C. Norman. fluid '(zlist); global '(!*trint); exports int!-fac,var2df; imports cubic,df2q,f2df,interr,multdf,printdf,quadratic,quartic,unifac, uniform,vp1,vp2,sub1; symbolic procedure int!-fac x; % Input: primitive, square-free polynomial (s.form). %output: % list of 'factors' wrt zlist % each item in this list is either % log . sq % or atan . sq % and these logs and arctans are all that is needed in the % integration of 1/(argument). begin scalar res,pol,dset,var,degree,vars; pol:=f2df x; %convert to distributed form dset:=degreeset(pol); %now extract factors of the form 'x' or 'log(x)' etc; %these correspond to items in dset with a non-zero cdr. begin scalar zl,ds; zl:=zlist; ds:=dset; while not null ds do << if onep cdar ds then << res:=('log . var2df(car zl,1,zlist)) . res; %record in answer. pol:=multdf(var2df(car zl,-1,zlist),pol); %divide out. if !*trint then << printc "Trivial factor found"; printdf cdar res>>; rplaca(ds,sub1 caar ds . cdar ds) >> else if null zerop cdar ds then interr "Repeated trivial factor in arg to factor"; zl:=cdr zl; ds:=cdr ds >>; end; %single term factors all removed now. dset:=mapcar(dset,function car); %get lower bounds. if !*trint then printc ("Upper bounds of remaining factors are now: " . dset); if dset=vp2 zlist then go to finished; %thing left is constant. begin scalar ds,zl; var:=car zlist; degree:=car dset; if not zerop degree then vars:=var . vars; ds:=cdr dset; zl:=cdr zlist; while not null ds do << if not zerop car ds then << vars:=car zl . vars; if zerop degree or degree>car ds then << var:=car zl; degree:=car ds >> >>; zl:=cdr zl; ds:=cdr ds >> end; % Now var is variable that this poly involves to lowest degree. % Degree is the degree of the poly in same variable. if !*trint then printc ("Best var is " . var . "with exponent " . degree); if onep degree then << res:=('log . pol) . res; %certainly irreducible. if !*trint then << printc "The following is certainly irreducible"; printdf pol>>; go to finished >>; if degree=2 then << if !*trint then << printc "Quadratic"; printdf pol>>; res:=quadratic(pol,var,res); go to finished >>; dset:=uniform(pol,var); if not (dset='failed) then << if !*trint then << printc "Univariate polynomial"; printdf pol >>; res:=unifac(dset,var,degree,res); go to finished >>; if not null cdr vars then go to nasty; %only try univariate now. if degree=3 then << if !*trint then << printc "Cubic"; printdf pol>>; res:=cubic(pol,var,res); % if !*overlaymode then excise 'd3d4; go to finished >>; if degree=4 then << if !*trint then << printc "Quartic"; printdf pol>>; res:=quartic(pol,var,res); % if !*overlaymode then excise 'd3d4; go to finished>>; %else abandon hope and hand back some rubbish. nasty: res:=('log . pol) . res; printc "The following polynomial has not been properly factored"; printdf pol; go to finished; finished: %res is a list of d.f. s as required pol:=nil; %convert back to standard forms while not null res do begin scalar type,arg; type:=caar res; arg:=cdar res; arg:=df2q arg; if type='log then rplacd(arg,1); pol:=(type . arg) . pol; res:=cdr res end; return pol end; symbolic procedure var2df(var,n,zlist); ((vp1(var,n,zlist) .* (1 ./ 1)) .+ nil); symbolic procedure degreeset pol; % Finds degree bounds for all vars in distributed form poly. degreesub(dbl lpow pol,red pol); symbolic procedure dbl x; % Converts list of x into list of (x . x). if null x then nil else (car x . car x) . dbl cdr x; symbolic procedure degreesub(cur,pol); % Update degree bounds 'cur' to include info about pol. << while not null pol do << cur:=degreesub1(cur,lpow pol); pol:=red pol >>; cur >>; symbolic procedure degreesub1(cur,nxt); %Merge information from exponent set next into cur. if null cur then nil else degreesub2(car cur,car nxt) . degreesub1(cdr cur,cdr nxt); symbolic procedure degreesub2(two,one); max(car two,one) . min(cdr two,one); endmodule; module ibasics; % Some basic support routines for integrator. % Authors: Mary Ann Moore and Arthur C. Norman. fluid '(!*backtrace !*gcd !*sqfree indexlist sqrtflag sqrtlist varlist zlist); global '(!*trint); exports partialdiff,printdf,rationalintegrate,interr; imports df2printform,printsf,varsinsf,addsq,multsq,multd,mksp; symbolic procedure printdf u; % Print distributed form via cheap conversion to reduce structure. begin scalar !*gcd; printsf df2printform u; end; symbolic procedure !*n2sq(u1); if u1=0 then nil . 1 else u1 . 1; symbolic procedure indx(n); if n<2 then (list 1) else(n . indx(isub1 n)); symbolic procedure interr mess; <<(!*trint or !*backtrace) and <<prin2 "***** INTEGRATION PACKAGE ERROR: "; printc mess>>; error1()>>; symbolic procedure rationalintegrate(x,var); begin scalar n,d; n:=numr x; d:=denr x; if not var member varsinsf(d,nil) then return !*multsq(polynomialintegrate(n,var),1 ./ d); rederr "Rational integration not coded yet" end; symbolic procedure polynomialintegrate(x,v); % Integrate standard form. result is standard quotient. if null x then nil ./ 1 else if atom x then ((mksp(v,1) .* 1) .+ nil) ./ 1 else begin scalar r; r:=polynomialintegrate(red x,v); % deal with reductum if v=mvar x then begin scalar degree,newlt; degree:=1+tdeg lt x; newlt:=((mksp(v,degree) .* lc x) .+ nil) ./ 1; % up exponent r:=addsq(!*multsq(newlt,1 ./ degree),r) end else begin scalar newterm; newterm:=(((lpow x) .* 1) .+ nil) ./ 1; newterm:=!*multsq(newterm,polynomialintegrate(lc x,v)); r:=addsq(r,newterm) end; return r end; symbolic procedure partialdiff(p,v); % Partial differentiation of p wrt v - p is s.f. as is result. if domainp p then nil else if v=mvar p then (lambda x; if x=1 then lc p else ((mksp(v,x-1) .* multd(x,lc p)) .+ partialdiff(red p,v))) (tdeg lt p) else (lambda x; if null x then partialdiff(red p,v) else ((lpow p .* x) .+ partialdiff(red p,v))) (partialdiff(lc p,v)); put('pdiff,'simpfn,'simppdiff); symbolic procedure ncdr(l,n); % we can use small integer arithmetic here. if n=0 then l else ncdr(cdr l,isub1 n); symbolic procedure mkilist(old,term); if null old then nil else term.mkilist(cdr old,term); %symbolic procedure addin(lista,first,listb); %if null lista % then nil % else ((first.car listb).car lista).addin(cdr lista,first,cdr listb); symbolic procedure removeduplicates(u); % Purges duplicates from the list passed to it. if null u then nil else if (atom u) then u.nil else if member(car u,cdr u) then removeduplicates cdr u else (car u).removeduplicates cdr u; symbolic procedure jsqfree(sf,var); begin varlist:=getvariables(sf ./ 1); zlist:=findzvars(varlist,list var,var,nil); sqrtlist:=findsqrts varlist; % before the purge sqrtflag:=not null sqrtlist; varlist:=purge(zlist,varlist); if sf eq !*sqfree then return list list sf else return sqfree(sf,zlist) end; symbolic procedure jfactor(sf,var); begin scalar varlist,zlist,indexlist,sqrtlist,sqrtflag; scalar prim,answer,u,v; % scalar var2 prim:=jsqfree(sf,var); indexlist:=createindices zlist; prim:=factorlistlist (prim,t); while prim do << if caar prim eq 'log then << u:=cdar prim; u:=!*multsq(numr u ./ 1,1 ./ cdr stt(numr u,var)); v:=denr u; while not atom v do v:=lc v; if (numberp v) and (0> v) then u:=(negf numr u) ./ (negf denr u); answer:=u.answer >> else if caar prim eq 'atan then << % We can ignore this, since we also get LOG (U**2+1) as another term. >> else interr "Unexpected term in jfactor"; prim:=cdr prim >>; return answer end; symbolic procedure stt(u,x); if domainp u then if u eq nil then ((-1) . nil) else (0 . u) else if mvar u eq x then ldeg u . lc u else if ordop(x,mvar u) then (0 . u) else begin scalar ltlc,ltrest; ltlc:=stt(lc u,x); ltrest:= stt(red u,x); if car ltlc = car ltrest then go to merge; if car ltlc > car ltrest then return car ltlc . !*multf(cdr ltlc,(lpow u .* 1) .+ nil) else return ltrest; merge: return car ltlc.addf(cdr ltrest, !*multf(cdr ltlc,(lpow u .* 1) .+ nil)) end; symbolic procedure gcdinonevar(u,v,x); % Gcd of u and v, regarded as polynnmials in x. if null u then v else if null v then u else begin scalar u1,v1,z,w; u1:=stt(u,x); v1:=stt(v,x); loop: if (car u1) > (car v1) then go to ok; z:=u1;u1:=v1;v1:=z; z:=u;u:=v;v:=z; ok: if car v1 iequal 0 then interr "Coprimeness in gcd"; z:=gcdf(cdr u1,cdr v1); w:=quotf(cdr u1,z); if (car u1) neq (car v1) then w:=!*multf(w,!*p2f mksp(x,(car u1)-(car v1))); u:=addf(!*multf(v,w), !*multf(u,negf quotf(cdr v1,z))); if null u then return v; u1:=stt(u,x); go to loop end; symbolic procedure mapply(funct,l); if null l then rederr "Empty list to mapply" else if null cdr l then car l else apply(funct,list(car l,mapply(funct,cdr l))); symbolic procedure !*lcm!*(u,v); !*multf(u,quotf(v,gcdf(u,v))); symbolic procedure multsql(u,l); % Multiplies (!*multsq) each element of l by u. if null l then nil else !*multsq(u,car l).multsql(u,cdr l); symbolic procedure intersect(x,y); if null x then nil else if member(car x,y) then car(x) . intersect(cdr x,y) else intersect(cdr x,y); symbolic procedure mapvec(v,f); begin scalar n; n:=upbv v; for i:=0:n do apply(f,list getv(v,i)) end; endmodule; module jpatches; % Routines for manipulating sf's with power folding. % Author: James H. Davenport. fluid '(sqrtflag); exports !*addsq,!*multsq,!*invsq,!*multf,squashsqrtsq,!*exptsq,!*exptf; % !*MULTF(A,B) multiplies the polynomials (standard forms) U and V % in much the same way as MULTF(U,V) would, EXCEPT... % (1) !*MULTF inhibits the action of OFF EXP and of non-commutative % multiplications % (2) Within !*MULTF powers of square roots, and powers of % exponential kernels are reduced as if substitution rules % such as FOR ALL X LET SQRT(X)**2=X were being applied; % Note that !*MULTF comes between MULTF and !*Q2F SUBS2F MULTF in its % behaviour, and that it is the responsibility of the user to call it % in sensible places where its services are needed; % similarly for the other functions defined here; %symbolic procedure !*addsq(u,v); %U and V are standard quotients. % %Value is canonical sum of U and V; % if null numr u then v % else if null numr v then u % else if denr u=1 and denr v=1 then addf(numr u,numr v) ./ 1 % else begin scalar nu,du,nv,dv,x; % x := gcdf(du:=denr u,dv:=denr v); % du:=quotf(du,x); dv:=quotf(dv,x); % nu:=numr u; nv:=numr v; % u:=addf(!*multf(nu,dv),!*multf(nv,du)); % if u=nil then return nil ./ 1; % v:=!*multf(du,denr v); % return !*ff2sq(u,v) % end; %symbolic procedure !*multsq(a,b); % begin % scalar n,d; % n:=!*multf(numr a,numr b); % d:=!*multf(denr a,denr b); % return !*ff2sq(n,d) % end; %symbolic procedure !*ff2sq(a,b); % begin % scalar gg; % if null a then return nil ./ 1; % gg:=gcdf(a,b); % if not (gg=1) then << % a:=quotf(a,gg); % b:=quotf(b,gg) >>; % if minusf b then << % a:=negf a; % b:=negf b >>; % return a ./ b % end; symbolic procedure !*addsq(u,v); %U and V are standard quotients. %Value is canonical sum of U and V; if null numr u then v else if null numr v then u else if denr u=1 and denr v=1 then addf(numr u,numr v) ./ 1 else begin scalar du,dv,x,y,z; x := gcdf(du:=denr u,dv:=denr v); du:=quotf(du,x); dv:=quotf(dv,x); y:=addf(!*multf(dv,numr u),!*multf(du,numr v)); if null y then return nil ./ 1; z:=!*multf(denr u,dv); if minusf z then <<y := negf y; z := negf z>>; if x=1 then return y ./ z; x := gcdf(y,x); return if x=1 then y ./ z else quotf(y,x) ./ quotf(z,x) end; symbolic procedure !*multsq(u,v); %U and V are standard quotients. Result is the canonical product of %U and V with surd powers suitably reduced. if null numr u or null numr v then nil ./ 1 else if denr u=1 and denr v=1 then !*multf(numr u,numr v) ./ 1 else begin scalar w,x,y; x := gcdf(numr u,denr v); y := gcdf(numr v,denr u); w := !*multf(quotf(numr u,x),quotf(numr v,y)); x := !*multf(quotf(denr u,y),quotf(denr v,x)); if minusf x then <<w := negf w; x := negf x>>; y := gcdf(w,x); % another factor may have been generated. return if y=1 then w ./ x else quotf(w,y) ./ quotf(x,y) end; symbolic procedure !*invsq a; % Note that several examples (e.g., int(1/(x**8+1),x)) give a more % compact result when SQRTFLAG is true if SQRT2TOP is not called. if sqrtflag then sqrt2top invsq a else invsq a; symbolic procedure !*multf(u,v); % U and V are standard forms % Value is SF for U*V; begin scalar x,y; if null u or null v then return nil else if u = 1 then return squashsqrt v else if v = 1 then return squashsqrt u else if domainp u then return multd(u,squashsqrt v) else if domainp v then return multd(v,squashsqrt u); x:=mvar u; y:=mvar v; if x eq y then go to c else if ordop(x,y) then go to b; x:=!*multf(u,lc v); y:=!*multf(u,red v); return if null x then y else if not domainp lc v and mvar u eq mvar lc v and not atom mvar u and car mvar u memq '(expt sqrt) then addf(!*multf(x,!*p2f lpow v),y) % what about noncom? else makeupsf(lpow v,x,y); b: x:=!*multf(lc u,v); y:=!*multf(red u,v); return if null x then y else if not domainp lc u and mvar lc u eq mvar v and not atom mvar v and car mvar v memq '(expt sqrt) then addf(!*multf(!*p2f lpow u,x),y) else makeupsf(lpow u,x,y); c: y:=addf(!*multf(list lt u,red v),!*multf(red u,v)); if eqcar(x,'sqrt) then return addf(squashsqrt y,!*multfsqrt(x, !*multf(lc u,lc v),ldeg u + ldeg v)) else if eqcar(x,'expt) and prefix!-rational!-numberp caddr x then return addf(squashsqrt y,!*multfexpt(x, !*multf(lc u,lc v),ldeg u + ldeg v)); x:=mkspm(x,ldeg u + ldeg v); return if null x or null (u:=!*multf(lc u,lc v)) then y else x .* u .+ y end; symbolic procedure makeupsf(u,x,y); % Makes u .* x .+ y except when u is not a valid lpow (because of % sqrts). if atom car u or cdr u = 1 then u .* x .+ y else if caar u eq 'sqrt then addf(!*multfsqrt(car u,x,cdr u),y) else if <<begin scalar v; v:=car u; if car v neq 'expt then return nil; v:=caddr v; if atom v then return nil; return (car v eq 'quotient and numberp cadr v and numberp caddr v) end >> then addf(!*multfexpt(car u,x,cdr u),y) else u .* x .+ y; symbolic procedure !*multfsqrt(x,u,w); % This code (Due to Norman a& Davenport) squashes SQRT(...)**2. begin scalar v; w:=divide(w,2); v:=!*q2f simp cadr x; u:=!*multf(u,exptf(v,car w)); if not zerop cdr w then u:=!*multf(u,!*p2f mksp(x,1)); return u end; symbolic procedure !*multfexpt(x,u,w); begin scalar expon,v; expon:=caddr x; x:=cadr x; w:=w * cadr expon; expon:=caddr expon; v:=gcdn(w,expon); w:=w/v; v:=expon/v; if not (w > 0) then rederr "Invalid exponent" else if v = 1 then return !*multf(u,exptf(if numberp x then x else if atom x then !*k2f x else !*q2f if car x eq '!*sq then argof x else simp x, w)); expon:=0; while not (w < v) do <<expon:=expon + 1; w:=w-v>>; if expon>0 then u:=!*multf(u,exptf(!*q2f simp x,expon)); if w = 0 then return u; x:=list('expt,x,list('quotient,1,v)); return multf(squashsqrt u,!*p2f mksp(x,w)) end; symbolic procedure prefix!-rational!-numberp u; % Tests for m/n in prefix representation. eqcar(u,'quotient) and numberp cadr u and numberp caddr u; symbolic procedure squashsqrtsq sq; !*multsq(squashsqrt numr sq ./ 1, 1 ./ squashsqrt denr sq); symbolic procedure squashsqrt sf; if (not sqrtflag) or atom sf or atom mvar sf then sf else if car mvar sf eq 'sqrt and ldeg sf > 1 then addf(squashsqrt red sf,!*multfsqrt(mvar sf,lc sf,ldeg sf)) else if car mvar sf eq 'expt and prefix!-rational!-numberp caddr mvar sf and ldeg sf >= caddr caddr mvar sf then addf(squashsqrt red sf,!*multfexpt(mvar sf,lc sf,ldeg sf)) else (lpow sf .* squashsqrt lc sf) .+ squashsqrt red sf; %remd 'simpx1; %symbolic procedure simpx1(u,m,n); % %u,m and n are prefix expressions; % %value is the standard quotient expression for u**(m/n); % begin scalar flg,z; % if null frlis!* or null xn(frlis!*,flatten (m . n)) % then go to a; % exptp!* := t; % return !*k2q list('expt,u,if n=1 then m % else list('quotient,m,n)); % a: if numberp m and fixp m then go to e % else if atom m then go to b % else if car m eq 'minus then go to mns % else if car m eq 'plus then go to pls % else if car m eq 'times and numberp cadr m and fixp cadr m % and numberp n % then go to tms; % b: z := 1; % c: if atom u and not numberp u then flag(list u,'used!*); % u := list('expt,u,if n=1 then m else list('quotient,m,n)); % if not u member exptl!* then exptl!* := u . exptl!*; % d: return mksq(u,if flg then -z else z); %u is already in lowest %% %terms; % e: if numberp n and fixp n then go to int; % z := m; % m := 1; % go to c; % mns: m := cadr m; % if !*mcd then return invsq simpx1(u,m,n); % flg := not flg; % go to a; % pls: z := 1 ./ 1; % pl1: m := cdr m; % if null m then return z; % z := multsq(simpexpt list(u, % list('quotient,if flg then list('minus,car m) % else car m,n)), % z); % go to pl1; % tms: z := gcdn(n,cadr m); % n := n/z; % z := cadr m/z; % m := retimes cddr m; % go to c; % int:z := divide(m,n); % if cdr z<0 then z:= (car z - 1) . (cdr z+n); % if 0 = cdr z % then return simpexpt list(u,car z); % if n = 2 % then return multsq(simpexpt list(u,car z), % simpsqrti u); % return multsq(simpexpt list(u,car z), % mksq(list('expt,u,list('quotient,1,n)),cdr z)) % end; symbolic procedure !*exptsq(a,n); % raise A to the power N using !*MULTSQ; if n=0 then 1 ./ 1 else if n=1 then a else if n<0 then !*exptsq(invsq a,-n) else begin scalar q,r; q:=divide(n,2); r:=cdr q; q:=car q; q:=!*exptsq(!*multsq(a,a),q); if r=0 then return q else return !*multsq(a,q) end; symbolic procedure !*exptf(a,n); % raise A to the power N using !*MULTF; if n=0 then 1 else if n=1 then a else begin scalar q,r; q:=divide(n,2); r:=cdr q; q:=car q; q:=!*exptf(!*multf(a,a),q); if r=0 then return q else return !*multf(a,q) end; endmodule; module hacksqrt; % Routines for manipulation of sqrt expressions. % Author: James H. Davenport. fluid '(nestedsqrts thisplace); exports sqrtsintree,sqrtsinsq,sqrtsinsql,sqrtsinsf,sqrtsign; exports degreenest,sortsqrts; imports mkvect,interr,getv,dependsp,union; symbolic procedure sqrtsintree(u,var,slist); % Adds to slist all the sqrts in the prefix-type tree u. if atom u then slist else if car u eq '!*sq then union(slist,sqrtsinsq(cadr u,var)) else if car u eq 'sqrt then if dependsp(argof u,var) then << slist:=sqrtsintree(argof u,var,slist); % nested square roots if member(u,slist) then slist else u.slist >> else slist else sqrtsintree(car u,var,sqrtsintree(cdr u,var,slist)); symbolic procedure sqrtsinsq(u,var); % Returns list of all sqrts in sq. sqrtsinsf(denr u,sqrtsinsf(numr u,nil,var),var); symbolic procedure sqrtsinsql(u,var); % Returns list of all sqrts in sq list. if null u then nil else sqrtsinsf(denr car u, sqrtsinsf(numr car u,sqrtsinsql(cdr u,var),var),var); symbolic procedure sqrtsinsf(u,slist,var); % Adds to slist all the sqrts in sf. if domainp u or null u then slist else << if eqcar(mvar u,'sqrt) and dependsp(argof mvar u,var) and not member(mvar u,slist) then begin scalar slist2; slist2:=sqrtsintree(argof mvar u,var,nil); if slist2 then << nestedsqrts:=t; slist:=union(slist2,slist) >>; slist:=(mvar u).slist end; sqrtsinsf(lc u,sqrtsinsf(red u,slist,var),var) >>; symbolic procedure easysqrtsign(slist,things); % This procedure builds a list of all substitutions for all possible % combinations of square roots in list. if null slist then things else easysqrtsign(cdr slist, nconc(mapcons(things,(car slist).(car slist)), mapcons(things, list(car slist,'minus,car slist)))); symbolic procedure hardsqrtsign(slist,things); % This procedure fulfils the same role for nested sqrts % ***assumption: the simpler sqrts come further up the list. if null slist then things else begin scalar thisplace,answers,pos,neg; thisplace:=car slist; answers:=mapcar(things,function (lambda u;sublis(u,thisplace).u)); pos:=mapcar(answers,function (lambda u;(thisplace.car u).(cdr u))); % pos is sqrt(f) -> sqrt(innersubst f) neg:=mapcar(answers, function (lambda u;list(thisplace,'minus,car u).(cdr u))); % neg is sqrt(f) -> -sqrt(innersubst f) return hardsqrtsign(cdr slist,nconc(pos,neg)) end; symbolic procedure degreenest(pf,var); % Returns the maximum degree of nesting of var % inside sqrts in the prefix form pf. if atom pf then 0 else if car pf eq 'sqrt then if dependsp(cadr pf,var) then iadd1 degreenest(cadr pf,var) else 0 else if car pf eq 'expt then if dependsp(cadr pf,var) then if eqcar(caddr pf,'quotient) then iadd1 degreenest(cadr pf,var) else degreenest(cadr pf,var) else 0 else degreenestl(cdr pf,var); symbolic procedure degreenestl(u,var); %Returns max degreenest from list of pfs u. if null u then 0 else max(degreenest(car u,var), degreenestl(cdr u,var)); symbolic procedure sortsqrts(u,var); % Sorts list of sqrts into order required by hardsqrtsign % (and many other parts of the package). begin scalar i,v; v:=mkvect(10); %should be good enough! while u do << i:=degreenest(car u,var); if i iequal 0 then interr "Non-dependent sqrt found"; if i > 10 then interr "Degree of nesting exceeds 10 (recompile with 10 increased)"; putv(v,i,(car u).getv(v,i)); u:=cdr u >>; u:=getv(v,10); for i :=9 step -1 until 1 do u:=nconc(getv(v,i),u); return u end; symbolic procedure sqrtsign(sqrts,x); if nestedsqrts then hardsqrtsign(sortsqrts(sqrts,x),list nil) else easysqrtsign(sqrts,list nil); endmodule; module kron; % Kronecker factorization of univ. polys over integers. % Authors: Mary Ann Moore and Arthur C. Norman. global '(!*trint); exports linfac,quadfac; imports zfactor; % Only linear and quadratic factors are found. symbolic procedure linfac(w); trykr(w,'(0 1)); symbolic procedure quadfac(w); trykr(w,'(-1 0 1)); symbolic procedure trykr(w,points); %Look for factor of w by evaluation at (points) and use of % interpolate. Return (fac . cofac) with fac=nil if none % found and cofac=nil if nothing worthwhile is left. begin scalar values,attempt; if null w then return nil . nil; if (length points > car w) then return w . nil; %that says if w is already tiny, it is already factored. values:=mapcar(points,function (lambda x; evalat(w,x))); if !*trint then << printc ("At x= " . points); printc ("p(x)= " . values)>>; if 0 member values then go to lucky; %(x-1) is a factor! values:=mapcar(values,function zfactors); rplacd(values,mapcar(cdr values,function (lambda y; append(y,mapcar(y,function minus))))); if !*trint then <<printc "Possible factors go through some of"; print values>>; attempt:=search4fac(w,values,nil); if null attempt then attempt:=nil . w; return attempt; lucky: %here (x-1) is a factor because p(0) or p(1) or p(-1) %vanished and cases p(0), p(-1) will have been removed %elsewhere. attempt:='(1 1 -1); %the factor return attempt . testdiv(w,attempt) end; symbolic procedure zfactors n; % Produces a list of all (positive) integer factors of the integer n. if n=0 then list 0 else if (n:=abs n)=1 then list 1 else combinationtimes zfactor n; symbolic procedure search4fac(w,values,cv); % Combinatorial search. cv gets current selected value-set. % Returns nil if fails, else factor . cofactor. if null values then tryfactor(w,cv) else begin scalar ff,q; ff:=car values; %try all values here loop: if null ff then return nil; %no factor found q:=search4fac(w,cdr values,(car ff) . cv); if null q then << ff:=cdr ff; go to loop>>; return q end; symbolic procedure tryfactor(w,cv); % Tests if cv represents a factor of w. begin scalar ff,q; if null cddr cv then ff:=linethrough(cadr cv,car cv) else ff:=quadthrough(caddr cv,cadr cv,car cv); if ff='failed then return nil; %it does not interpolate q:=testdiv(w,ff); if q='failed then return nil; %not a factor return ff . q end; symbolic procedure evalat(p,n); % Evaluate polynomial at integer point n. begin scalar r; r:=0; p:=cdr p; while not null p do << r:=n*r+car p; p:=cdr p >>; return r end; symbolic procedure testdiv(a,b); % Quotient a/b or failed. begin scalar q; q:=testdiv1(cdr a,car a,cdr b,car b); if q='failed then return q; return (car a-car b) . q end; symbolic procedure testdiv1(a,da,b,db); if da<db then begin check0: if null a then return nil else if not zerop car a then return 'failed; a:=cdr a; go to check0 end else begin scalar q; q:=divide(car a,car b); if zerop cdr q then q:=car q else return 'failed; a:=testdiv1(ambq(cdr a,cdr b,q),da-1,b,db); if a='failed then return a; return q . a end; symbolic procedure ambq(a,b,q); % A-B*Q with Q an integer. if null b then a else ((car a)-(car b)*q) . ambq(cdr a,cdr b,q); symbolic procedure linethrough(y0,y1); begin scalar a; a:=y1-y0; if zerop a then return 'failed; if a<0 then <<a:=-a; y0:=-y0 >>; if onep gcdn(a,y0) then return list(1,a,y0); return 'failed end; symbolic procedure quadthrough(ym1,y0,y1); begin scalar a,b,c; a:=divide(ym1+y1,2); if zerop cdr a then a:=(car a)-y0 else return 'failed; if zerop a then return 'failed; %linear things already done. c:=y0; b:=divide(y1-ym1,2); if zerop cdr b then b:=car b else return 'failed; if not onep gcdn(a,gcd(b,c)) then return 'failed; if a<0 then <<a:=-a; b:=-b; c:=-c>>; return list(2,a,b,c) end; endmodule; module lowdeg; % Splitting of low degree polynomials. % Author: To be determined. fluid '(clogflag knowndiscrimsign zlist); global '(!*trint); exports forceazero,makepolydf,quadratic,covecdf,exponentdf; imports dfquotdf,gcdf,interr,minusdfp,multdf,multdfconst,!*multf, negsq,minusp,printsq,multsq,invsq,pnth,nth,mknill, negdf,plusdf,printdf,printsq,quotf,sqrtdf,var2df,vp2,addsq,sub1; symbolic procedure covecdf(pol,var,degree); % Extract coefficients of polynomial wrt var, given a degree-bound % degree. Result is a lisp vector. begin scalar v,x,w; w:=pol; v:=mkvect(degree); while not null w do << x:=exponentof(var,lpow w,zlist); if (x<0) or (x>degree) then interr "Bad degree in covecdf"; putv(v,x,lt w . getv(v,x)); w:=red w >>; for i:=0:degree do putv(v,i,multdf(reversewoc getv(v,i), var2df(var,-i,zlist))); return v end; symbolic procedure quadratic(pol,var,res); % Add in to res logs or arctans corresponding to splitting the % polynomial. Pol given that it is quadratic wrt var. % Does not assume pol is univariate. begin scalar a,b,c,w,discrim; w:=covecdf(pol,var,2); a:=getv(w,2); b:=getv(w,1); c:=getv(w,0); % that split the quadratic up to find the coefficients a,b,c. if !*trint then << printc "a="; printdf a; printc "b="; printdf b; printc "c="; printdf c>>; discrim:=plusdf(multdf(b,b), multdfconst((-4) . 1,multdf(a,c))); if !*trint then << printc "Discriminant is"; printdf discrim>>; if null discrim then interr "Discrim=0 in quadratic"; if knowndiscrimsign then <<if knowndiscrimsign eq 'negative then go to atancase>> else if (not clogflag) and (minusdfp discrim) then go to atancase; discrim:=sqrtdf(discrim); if discrim='failed then go to nofactors; if !*trint then << printc "Square root is"; printdf discrim>>; w:=var2df(var,1,zlist); w:=multdf(w,a); b:=multdfconst(1 ./ 2,b); discrim:=multdfconst(1 ./ 2,discrim); w:=plusdf(w,b); %a*x+b/2. a:=plusdf(w,discrim); b:=plusdf(w,negdf(discrim)); if !*trint then << printc "Factors are"; printdf a; printdf b>>; return ('log . a) . ('log . b) . res; atancase: discrim:=sqrtdf negdf discrim; %sqrt(4*a*c-b**2) this time! if discrim='failed then go to nofactors; %sqrt did not exist? res := ('log . pol) . res; %one part of the answer a:=multdf(a,var2df(var,1,zlist)); a:=plusdf(b,multdfconst(2 ./ 1,a)); a:=dfquotdf(a,discrim); %assumes division is exact return ('atan . a) . res; nofactors: if !*trint then <<printc "The following quadratic does not seem to factor"; printdf pol>>; return ('log . pol) . res end; symbolic procedure exponentof(var,l,zl); if null zl then interr "Var not found in exponentof" else if var=car zl then car l else exponentof(var,cdr l,cdr zl); symbolic procedure df2sf a; if null a then nil else if ((null red a) and (denr lc a = 1) and (lpow a=vp2 zlist)) then numr lc a else interr "Nasty cubic or quartic"; symbolic procedure makepolydf p; % Multiply df by lcm of denominators of all coefficient denominators. begin scalar h,w; if null(w:=p) then return nil; %poly is zero already. h:=denr lc w; %a good start. w:=red w; while not null w do << h:=quotf(!*multf(h,denr lc w),gcdf(h,denr lc w)); w:=red w >>; %h is now lcm of denominators. return multdfconst(h ./ 1,p) end; symbolic procedure forceazero(p,n); %Shift polynomial p so that coeff of x**(n-1) vanishes. %Return the amount of the shift, update (vector) p. begin scalar r,i,w; for i:=0:n do putv(p,i,df2sf getv(p,i)); %convert to polys. r:=getv(p,n-1); if null r then return nil ./ 1; %already zero. r:= !*multsq(r ./ 1,invsq(!*multf(n,getv(p,n)) ./ 1)); % Used to be subs2q multsq %the shift amount. %now I have to set p:=subst(x-r,x,p) and then reduce to sf again. if !*trint then << printc "Shift is by "; printsq r>>; w:=mkvect(n); %workspace vector. for i:=0:n do putv(w,i,nil ./ 1); %zero it. i:=n; while not minusp i do << mulvecbyxr(w,negsq r,n); %W:=(X-R)*W; putv(w,0,addsq(getv(w,0),getv(p,i) ./ 1)); i:=i-1 >>; if !*trint then << printc "SQ shifted poly is"; print w>>; for i:=0:n do putv(p,i,getv(w,i)); w:=denr getv(p,0); for i:=1:n do w:=quotf(!*multf(w,denr getv(p,i)), gcdf(w,denr getv(p,i))); for i:=0:n do putv(p,i,numr !*multsq(getv(p,i),w ./ 1)); % Used to be subs2q multsq w:=getv(p,0); for i:=1:n do w:=gcdf(w,getv(p,i)); if not (w=1) then for i:=0:n do putv(p,i,quotf(getv(p,i),w)); if !*trint then << printc "Final shifted poly is "; print p>>; return r end; symbolic procedure mulvecbyxr(w,r,n); % W is a vector representing a poly of degree n. % Multiply it by (x+r). begin scalar i,im1; i:=n; im1:=sub1 i; while not minusp im1 do << putv(w,i,!*addsq(getv(w,im1),!*multsq(r,getv(w,i)))); % Used to be subs2q addsq/multsq i:=im1; im1:=sub1 i >>; putv(w,0,!*multsq(getv(w,0),r)); % Used to be subs2q multsq return w end; endmodule; module reform; % Reformulate expressions using C-constant substitution. % Authors: Mary Ann Moore and Arthur C. Norman. fluid '(cmap cval loglist ulist); global '(!*trint); exports logstosq,substinulist; imports prepsq,mksp,nth,multsq,addsq,domainp,invsq,plusdf; symbolic procedure substinulist ulist; % Substitutes for the C-constants in the values of the U's given in % ULIST. Result is a D.F. if null ulist then nil else begin scalar temp,lcu; lcu:=lc ulist; temp:=evaluateuconst numr lcu; if null numr temp then temp:=nil else temp:=((lpow ulist) .* !*multsq(temp,!*invsq(denr lcu ./ 1))) .+ nil; return plusdf(temp,substinulist red ulist) end; symbolic procedure evaluateuconst coefft; % Substitutes for the C-constants into COEFFT (=S.F.). Result is S.Q.; if null coefft or domainp coefft then coefft ./ 1 else begin scalar temp; if null(temp:=assoc(mvar coefft,cmap)) then temp:=(!*p2f lpow coefft) ./ 1 else temp:=getv(cval,cdr temp); temp:=!*multsq(temp,evaluateuconst(lc coefft)); % Next line had addsq previously return !*addsq(temp,evaluateuconst(red coefft)) end; symbolic procedure logstosq; % Converts LOGLIST to sum of the log terms as a S.Q.; begin scalar lglst,logsq,i,temp; i:=1; lglst:=loglist; logsq:=nil ./ 1; loop: if null lglst then return logsq; temp:=cddr car lglst; if !*trint then <<printc "SF arg for log etc ="; printc temp>>; if not (caar lglst='iden) then << temp:=prepsq temp; %convert to prefix form. temp:=list(caar lglst,temp); %function name. temp:=((mksp(temp,1) .* 1) .+ nil) ./ 1 >>; temp:=!*multsq(temp,getv(cval,i)); % Next line had addsq previously logsq:=!*addsq(temp,logsq); lglst:=cdr lglst; i:=i+1; go to loop end; endmodule; module simplog; % Simplify logarithms. % Authors: Mary Ann Moore and Arthur C. Norman. exports simplog,simplogsq; imports quotf,prepf,mksp,simp!*,!*multsq,simptimes,addsq,minusf,negf, addf,comfac,negsq,mk!*sq,carx; symbolic procedure simplog(exxpr); simplogi carx(exxpr,'simplog); symbolic procedure simplogi(sq); begin if atom sq then go to simplify; if car sq eq 'times then return simpplus(for each u in cdr sq collect mk!*sq simplogi u); if car sq eq 'quotient then return addsq(simplogi cadr sq, negsq simplogi caddr sq); if car sq eq 'expt then return simptimes list(caddr sq, mk!*sq simplogi cadr sq); if car sq eq 'nthroot then return !*multsq(1 ./ caddr sq,simplogi cadr sq); % we had (nthroot of n). if car sq eq 'sqrt then return !*multsq(1 ./ 2,simplogi argof sq); if car sq = '!*sq then return simplogsq cadr sq; simplify: sq:=simp!* sq; return simplogsq sq end; symbolic procedure simplogsq sq; addsq((simplog2 numr sq),negsq(simplog2 denr sq)); symbolic procedure simplog2(sf); if atom sf then if null sf then rederr "Log 0 formed" else if numberp sf then if sf iequal 1 then nil ./ 1 else if sf iequal 0 then rederr "Log 0 formed" else((mksp(list('log,sf),1) .* 1) .+ nil) ./ 1 else formlog(sf) else begin scalar form; form:=comfac sf; if not null car form then return addsq(formlog(form .+ nil), simplog2 quotf(sf,form .+ nil)); % we have killed common powers. form:=cdr form; if form neq 1 then return addsq(simplog2 form, simplog2 quotf(sf,form)); % remove a common factor from the sf. return (formlog sf) end; symbolic procedure formlogterm(sf); begin scalar u; u:=mvar sf; if not atom u and (car u member '(times sqrt expt nthroot)) then u:=addsq(simplog2 lc sf, !*multsq(simplogi u,simp!* ldeg sf)) else if (lc sf iequal 1) and (ldeg sf iequal 1) then u:=((mksp(list('log,u),1) .* 1) .+ nil) ./ 1 else u:=addsq(simptimes list(list('log,u),ldeg sf), simplog2 lc sf); return u end; symbolic procedure formlog sf; if null red sf then formlogterm sf else if minusf sf then addf((mksp(list('log,-1),1) .* 1) .+ nil, formlog2 negf sf) ./ 1 else (formlog2 sf) ./ 1; symbolic procedure formlog2 sf; ((mksp(list('log,prepf sf),1) .* 1) .+ nil); endmodule; module simpsqrt; % Simplify square roots. % Authors: Mary Ann Moore and Arthur C. Norman. % Heavily modified J.H. Davenport for algebraic functions. fluid '(!*backtrace !*conscount !*galois !*pvar basic!-listofallsqrts gaussiani basic!-listofnewsqrts kord!* knowntobeindep listofallsqrts listofnewsqrts sqrtflag sqrtlist sqrt!-places!-alist variable varlist zlist); global '(!*tra !*trint); % This module should be rewritten in terms of the REDUCE function % SIMPSQRT; % remd 'simpsqrt; exports proper!-simpsqrt,simpsqrti,simpsqrtsq,simpsqrt2,sqrtsave, newplace,actualsimpsqrt,formsqrt; symbolic procedure proper!-simpsqrt(exprn); simpsqrti carx(exprn,'proper!-simpsqrt); symbolic procedure simpsqrti sq; begin scalar u; if atom sq then if numberp sq then return (simpsqrt2 sq) ./ 1 else if (u:=get(sq,'avalue)) then return simpsqrti cadr u % BEWARE!!! This is VERY system dependent. else return simpsqrt2((mksp(sq,1) .* 1) .+ nil) ./ 1; % If it doesnt have an AVALUE then it is itself; if car sq eq 'times then return mapply(function multsq, mapcar(cdr sq,function simpsqrti)); if car sq eq 'quotient then return multsq(simpsqrti cadr sq, invsq simpsqrti caddr sq); if car sq eq 'expt and numberp caddr sq then if evenp caddr sq then return simpexpt list(cadr sq,caddr sq / 2) else return simpexpt list(mk!*sq simpsqrti cadr sq,caddr sq); if car sq = '!*sq then return simpsqrtsq cadr sq; return simpsqrtsq tidysqrt simp!* sq end; symbolic procedure simpsqrtsq sq; (simpsqrt2 numr sq) ./ (simpsqrt2 denr sq); symbolic procedure simpsqrt2 sf; if minusf sf then if sf iequal -1 then gaussiani else begin scalar u; u:=negf sf; if numberp u then return multf(gaussiani,simpsqrt3 u); % we cannot negate general expressions for the following reason: % (%%%thesis remark%%%) % sqrt(x*x-1) under x->1/x gives sqrt(1-x*x)/x=i*sqrt(x*x-1)/x % under x->1/x gives x*i*sqrt(-1+1/x*x)=i**2*sqrt(x*x-1) % hence an abysmal catastrophe; return simpsqrt3 sf end else simpsqrt3 sf; symbolic procedure simpsqrt3 sf; begin scalar u; u:=assoc(sf,listofallsqrts); if u then return cdr u; % now see if 'knowntobeindep'can help. u:=atsoc(listofnewsqrts,knowntobeindep); if null u then go to no; u:=assoc(sf,cdr u); if u then << listofallsqrts:=u.listofallsqrts; return cdr u >>; no: u:=actualsimpsqrt sf; listofallsqrts:=(sf.u).listofallsqrts; return u end; symbolic procedure sqrtsave(u,v,place); begin scalar a; %u is new value of listofallsqrts, v of new. a:=assoc(place,sqrt!-places!-alist); if null a then sqrt!-places!-alist:=(place.(listofnewsqrts.listofallsqrts)) .sqrt!-places!-alist else rplacd(a,listofnewsqrts.listofallsqrts); listofnewsqrts:=v; % throw away things we are not going to need in future. if not !*galois then listofallsqrts:=u; % we cannot guarantee the validity of our calculations. if listofallsqrts eq u then return nil; v:=listofallsqrts; while not (cdr v eq u) do v:=cdr v; rplacd(v,nil); % listofallsqrts is all those added since routine was entered. v:=atsoc(listofnewsqrts,knowntobeindep); if null v then knowntobeindep:=(listofnewsqrts.listofallsqrts) . knowntobeindep else rplacd(v,union(cdr v,listofallsqrts)); listofallsqrts:=u; return nil end; symbolic procedure newplace(u); % Says to restart algebraic bases at a new place u. begin scalar v; v:=assoc(u,sqrt!-places!-alist); if null v then << listofallsqrts:=basic!-listofallsqrts; listofnewsqrts:=basic!-listofnewsqrts >> else << v:=cdr v; listofnewsqrts:=car v; listofallsqrts:=cdr v >>; return if v then v else listofnewsqrts.listofallsqrts end; symbolic procedure mknewsqrt u; % U is prefix form. begin scalar v,w; if not !*galois then go to new; % no checking required. v:=addf(!*p2f mksp(!*pvar,2),negf !*q2f tidysqrt simp u); % count !*conscount; w:=errorset(list('afactor,mkquote v,mkquote !*pvar),t,!*backtrace); % if !*tra then << % princ "*** That took "; % princ (!*conscount - count nil); % printc " conses" >>; % count 0; if atom w then go to new else w:=car w; %the actual result of afactor. if cdr w then go to notnew; new: w:=sqrtify u; listofnewsqrts:=w . listofnewsqrts; return !*kk2f w; notnew: w:=car w; v:=stt(w,!*pvar); if car v neq 1 then rederr "HELP ..."; w:=addf(w,multf(cdr v,(mksp(!*pvar,car v) .* -1) .+nil)); w:=sqrt2top(w ./ cdr v); w:=quotf(numr w,denr w); if null w then rederr "Division failure"; return w end; symbolic procedure actualsimpsqrt(sf); if sf iequal -1 then gaussiani else actualsqrtinner(sf,listofnewsqrts); symbolic procedure actualsqrtinner(sf,l); if null l then actualsimpsqrt2 sf else begin scalar z; % z:=simp argof car l; % if denr z neq 1 or (z := numr z) iequal -1 z:=!*q2f simp argof car l; if z iequal -1 then return actualsqrtinner(sf,cdr l); z:=quotf(sf,z); if null z then return actualsqrtinner(sf,cdr l); return !*multf(!*kk2f car l,actualsimpsqrt z) end; symbolic procedure actualsimpsqrt2(sf); if atom sf then if null sf then nil else if numberp sf then if sf < 0 then multf(gaussiani,actualsimpsqrt2(- sf)) %Above 2 lines inserted JHD 4 Sept 80; % test case: SQRT(B*X**2-C)/SQRT(X); else begin scalar n; n:=int!-sqrt sf; % Changed for conformity with DEC20 LISP JHD July 1982; if not fixp n then return mknewsqrt sf else return n end else mknewsqrt(sf) else begin scalar form; form:=comfac sf; if car form then return multf((if null cdr sf and (car sf = form) then formsqrt(form .+ nil) else simpsqrt2(form .+ nil)), %The above 2 lines changed by JHD; %(following suggestions of Morrison); %to conform to Standard LISP 4 Sept 80; simpsqrt2 quotf(sf,form .+ nil)); % we have killed common powers. form:=cdr form; if form neq 1 then return multf(simpsqrt2 form, simpsqrt2 quotf(sf,form)); % remove a common factor from the sf. return formsqrt sf end; symbolic procedure int!-sqrt n; % Return sqrt of n if same is exact, or something non-numeric % otherwise. if not numberp n then 'nonnumeric else if n<0 then 'negative else if floatp n then sqrt!-float n else if n<2 then n else int!-nr(n,(n+1)/2); symbolic procedure int!-nr(n,root); % root is an overestimate here. nr moves downwards to root; begin scalar w; w:=root*root; if n=w then return root; w:=(root+n/root)/2; if w>=root then return !*q2f simpsqrt list n; return int!-nr(n,w) end; symbolic procedure formsqrt(sf); if (null red sf) then if (lc sf iequal 1) and (ldeg sf iequal 1) then mknewsqrt mvar sf else multf(if evenp ldeg sf then !*p2f mksp(mvar sf,ldeg sf / 2) else exptf(mknewsqrt mvar sf,ldeg sf),simpsqrt2 lc sf) else begin scalar varlist,zlist,sqrtlist,sqrtflag; scalar v,l,n,w; % This returns a list, the i-th member of which is % a list of the factors of multiplicity i (as s.f's); v:=jsqfree(sf,if variable and involvesf(sf,variable) then variable else findatom mvar sf); % VARIABLE is the best thing to do square-free % decompositions with respect to, but anything % else will do if VARIABLE is not set; if null cdr v and null cdar v then return mknewsqrt prepf sf; % The JSQFREE did nothing; l:=nil; n:=0; while v do << n:=n+1; w:=car v; while w do << l:=list('expt,mk!*sq !*f2q car w,n) . l; w:=cdr w >>; v:=cdr v >>; if null cdr l then l:=car l else l:='times.l; % makes L into a valid prefix form; return !*q2f simpsqrti l end; symbolic procedure findatom pf; if atom pf then pf else findatom argof pf; symbolic procedure sqrtify u; % Actually creates the sqrt. begin scalar v,n,prinlist; n:=degreenest(u,nil); % v:=list('sqrt,u); v := mksqrt u; % To ensure sqrt**2 rule set. prinlist:=member(v,kord!*); if prinlist then return car prinlist; % This might improve sharing. % anyway, we mustn't re-add this object to KORD!*; while kord!* and eqcar(car kord!*,'sqrt) and (degreenest(argof car kord!*,nil) > n) do << prinlist:=(car kord!*) . prinlist; kord!*:=cdr kord!* >>; kord!*:=v . kord!*; while prinlist do << kord!*:=(car prinlist) . kord!*; prinlist:=cdr prinlist >>; return v end; % deflist ('((sqrt (((x) quotient (sqrt x) (times 2 x))))),'dfn); % put('sqrt,'simpfn,'proper!-simpsqrt); % Now in driver. endmodule; module isolve; % Routines for solving the final reduction equation. % Author: Mary Ann Moore and Arthur C. Norman. % Modifications by: John P. Fitch. fluid '(badpart ccount cmap cmatrix cval indexlist lhs!* lorder nnn orderofelim pt rhs!* sillieslist tanlist ulist zlist); global '(!*number!* !*statistics !*trint); exports solve!-for!-u; imports nth,findpivot,gcdf,int!-gensym1,mkvect,interr,multdfconst, !*multf!*,negdf,orddf,plusdf,printdf,printsf,printspreadc,printsq, quotf,putv,spreadc,subst4eliminatedcs,mknill,pnth,domainp,addf, invsq,multsq; symbolic procedure uterm(powu,rhs!*); % Finds the contribution from RHS!* of reduction equation, of the; % U-coefficient given by POWU. Result is in D.F.; if null rhs!* then nil else begin scalar coef,power; power:=addinds(powu,lpow rhs!*); coef:=evaluatecoeffts(numr lc rhs!*,powu); if null coef then return uterm(powu,red rhs!*); coef:=coef ./ denr lc rhs!*; return plusdf((power .* coef) .+ nil,uterm(powu,red rhs!*)) end; symbolic procedure solve!-for!-u(rhs!*,lhs!*,ulist); % Solves the reduction eqn LHS!*=RHS!*. Returns list of U-coefficients % and their values (ULIST are those we have so far), and a list of % C-equations to be solved (CLIST are the eqns we have so far). if null lhs!* then ulist else begin scalar u,lpowlhs; lpowlhs:=lpow lhs!*; begin scalar ll,mm,chge; ll:=maxorder(rhs!*,zlist,0); mm:=lorder; while mm do << if car ll < car mm then << chge:=t; rplaca(mm,car ll) >>; ll:=cdr ll; mm:=cdr mm >>; if !*trint and chge then << print ("Maxorder now ".lorder) >> end; u:=pickupu(rhs!*,lpow lhs!*,t); if null u then << if !*trint then << printc "***** C-equation to solve:"; printsf numr lc lhs!*; printc " = 0"; printc " ">>; % Remove a zero constant from the lhs, rather than use % Gauss Elim; if gausselimn(numr lc lhs!*,lt lhs!*) then lhs!*:=squashconstants(red lhs!*) else lhs!*:=red lhs!* >> else << ulist:=(car u . !*multsq(coefdf(lhs!*,lpowlhs),invsq cdr u)).ulist; % used to be subs2q multsq if !*statistics then !*number!*:=!*number!*+1; if !*trint then <<prin2 "***** U"; prin2 car u; prin2t " ="; printsq multsq(coefdf(lhs!*,lpowlhs),invsq cdr u); printc " ">>; lhs!*:=plusdf(lhs!*, negdf multdfconst(cdar ulist,uterm(car u,rhs!*))) >>; if !*trint then << printc ".... LHS is now:"; printdf lhs!*; printc " ">>; return solve!-for!-u(rhs!*,lhs!*,ulist) end; symbolic procedure squashconstants(express); begin scalar constlst,ii,xp,cl,subby,cmt,xx; constlst:=reverse cmap; cmt:=cmatrix; xxx: xx:=car cmt; % Look at next row of Cmatrix; cl:=constlst; % and list of the names; ii:=1; % will become index of removed constant; while not getv(xx,ii) do << ii:=ii+1; cl:=cdr cl >>; subby:=caar cl; %II is now index, and SUBBY the name; if member(subby,sillieslist) then <<cmt:=cdr cmt; go to xxx>>; %This loop must terminate; % This is because at least one constant remains; xp:=prepsq !*f2q getv(xx,0); % start to build up the answer; cl:=cdr cl; if not (ccount=ii) then for jj:=ii+1:ccount do << if getv(xx,jj) then xp:=list('plus,xp, list('times,caar cl, prepsq !*f2q getv(xx,jj))); cl:=cdr cl >>; xp:=list('quotient,list('minus,xp), prepsq !*f2q getv(xx,ii)); if !*trint then << prin2 "Replace "; prin2 subby; prin2 " by "; printsq simp xp >>; sillieslist:=subby . sillieslist; return subdf(express,xp,subby) end; symbolic procedure checku(ulist,u); % Checks that U is not already in ULIST - ie. that this u-coefficient; % has not already been given a value; if null ulist then nil else if (car u) = caar ulist then t else checku(cdr ulist,u); symbolic procedure checku1(powu,rhs!*); %Checks that use of a particular U-term will not cause trouble; %by introducing negative exponents into lhs when it is used; begin top: if null rhs!* then return nil; if negind(powu,lpow rhs!*) then if not null evaluatecoeffts(numr lc rhs!*,powu) then return t; rhs!*:=red rhs!*; go to top end; symbolic procedure negind(pu,pr); %check if substituting index values in power gives rise to -ve % exponents; if null pu then nil else if (car pu+caar pr)<0 then t else negind(cdr pu,cdr pr); symbolic procedure evaluatecoeffts(coefft,indlist); % Substitutes the values of the i,j,k,...'s that appear in the S.F. ; % COEFFT (=coefficient of r.h.s. of reduction equation). Result is S.F.; if null coefft or domainp coefft then if coefft=0 then nil else coefft else begin scalar temp; if mvar coefft member indexlist then temp:=valuecoefft(mvar coefft,indlist,indexlist) else temp:=!*p2f lpow coefft; temp:=!*multf(temp,evaluatecoeffts(lc coefft,indlist)); return addf(temp,evaluatecoeffts(red coefft,indlist)) end; symbolic procedure valuecoefft(var,indvalues,indlist); % Finds the value of VAR, which should be in INDLIST, given INDVALUES; % - the corresponding values of INDLIST variables; if null indlist then interr "Valuecoefft - no value" else if var eq car indlist then if car indvalues=0 then nil else car indvalues else valuecoefft(var,cdr indvalues,cdr indlist); symbolic procedure addinds(powu,powrhs); % Adds indices in POWU to those in POWRHS. Result is LPOW of D.F.; if null powu then if null powrhs then nil else interr "Powrhs too long" else if null powrhs then interr "Powu too long" else (car powu + caar powrhs).addinds(cdr powu,cdr powrhs); symbolic procedure pickupu(rhs!*,powlhs,flg); % Picks up the 'lowest' U coefficient from RHS!* if it exists and % returns it in the form of LT of D.F.. % Returns NIL if no legal term in RHS!* can be found. % POWLHS is the power we want to match (LPOW of D.F). % and COEFFU is the list of previous coefficients that must be zero; begin scalar coeffu,u; pt:=rhs!*; top: if null pt then return nil; %no term found - failed; u:=nextu(lt pt,powlhs); %check this term...; if null u then go to notthisone; if not testord(car u,lorder) then go to neverthisone; if not checkcoeffts(coeffu,car u) then go to notthisone; %that inhibited clobbering things already passed over; if checku(ulist,u) then go to notthisone; %that avoided redefining a u value; if checku1(car u,rhs!*) then go to neverthisone; %avoid introduction of negative exponents; if flg then u:=patchuptan(list u,powlhs,red pt,rhs!*); return u; neverthisone: coeffu:=(lc pt) . coeffu; notthisone: pt:=red pt; go to top end; symbolic procedure patchuptan(u,powlhs,rpt,rhs!*); begin scalar uu,cc,dd,tanlist,redu,redu1; pt:=rpt; while pt do << if (uu:=pickupu(pt,powlhs,nil)) and testord(car uu,lorder) then << % Nasty found, patch it up; cc:=(int!-gensym1('!C).caar u).cc; % CC is an alist of constants; if !*trint then <<prin2 "***** U"; prin2 caar u; prin2t " ="; print caar cc >>; redu:=plusdf(redu, multdfconst(!*k2q caar cc,uterm(caar u,rhs!*))); u:=uu.u >>; if pt then pt:=red pt >>; redu1:=redu; while redu1 do begin scalar xx; xx:=car redu1; if !*trint then << prin2 "Introduced residue "; print xx >>; if (not testord(car xx,lorder)) then << if !*trint then << printsq cdr xx; printc " = 0" >>; if dd:=killsingles(cadr xx,cc) then << redu:=subdf(redu,0,car dd); redu1:=subdf(redu1,0,car dd); ulist:=((cdr dd).(nil ./ 1)).ulist; u:=rmve(u,cdr dd); cc:=purgeconst(cc,dd) >> else redu1:=cdr redu1 >> else redu1:=cdr redu1 end; for each xx in redu do << if (not testord(car xx,lorder)) then << while cc do << addctomap(caar cc); ulist:=((cdar cc).(!*k2q caar cc)) . ulist; if !*statistics then !*number!*:=!*number!*+1; cc:=cdr cc >>; gausselimn(numr lc redu,lt redu)>> >>; if redu then << while cc do << addctomap(caar cc); ulist:=((cdar cc).(!*k2q caar cc)).ulist; if !*statistics then !*number!*:=!*number!*+1; cc:=cdr cc >>; lhs!*:=plusdf(lhs!*,negdf redu) >>; return car u end; symbolic procedure killsingles(xx,cc); if atom xx then nil else if not (cdr xx eq nil) then nil else begin scalar dd; dd:=assoc(caaar xx,cc); if dd then return dd; return killsingles(cdar xx,cc) end; symbolic procedure rmve(l,x); if caar l=x then cdr l else cons(car l,rmve(cdr l,x)); symbolic procedure subdf(a,b,c); % SUBSTITUTE B FOR C INTO THE DF A; % Used to get rid of silly constants introduced; if a=nil then nil else begin scalar x; x:=subs2q subf(numr lc a,list (c . b)) ; if x=(nil . 1) then return subdf(red a,b,c) else return plusdf( list ((lpow a).((car x).!*multf(cdr x,denr lc a))), subdf(red a,b,c)) end; symbolic procedure testord(a,b); % Test order of two DF's in recursive fashion; if null a then t else if car a leq car b then testord(cdr a,cdr b) else nil; symbolic procedure tanfrom(rhs!*,z,nnn); % We notice that in all bad cases we have (j-num)tan**j...; % Extract the num; begin scalar n,zz,r,rr; r:=rhs!*; n:=0; zz:=zlist; while car zz neq z do << n:=n+1; zz:=cdr zz >>; a: if null r then go to b; rr:=caar r; % The list of powers; for i:=1:n do rr:=cdr rr; if fixp caar rr then if caar rr>0 then << rr:=numr cdar r; if null red rr then rr:=nil ./ 1 else if fixp (rr:=quotf(red rr,lc rr)) then rr:=-rr else rr:=0>>; if atom rr then go to b; r:=cdr r; go to a; b: if null r then return maxfrom(lhs!*,nnn)+1; return max(rr,maxfrom(lhs!*,nnn)+1) end; symbolic procedure coefdf(y,u); if y=nil then nil else if lpow y=u then lc y else coefdf(red y,u); symbolic procedure purgeconst(a,b); % Remove a const from and expression. May be the same as DELETE?; if null a then nil else if car a=b then purgeconst(cdr a,b) else cons(car a,purgeconst(cdr a,b)); symbolic procedure maxorder(rhs!*,z,n); % Find a limit on the order of terms, theis is ad hoc; if null z then nil else if eqcar(car z,'sqrt) then cons(1,maxorder(rhs!*,cdr z,n+1)) else if (atom car z) or (caar z neq 'tan) then cons(maxfrom(lhs!*,n)+1,maxorder(rhs!*,cdr z,n+1)) else cons(tanfrom(rhs!*,car z,n),maxorder(rhs!*,cdr z,n+1)); symbolic procedure maxfrom(l,n); % Largest order in the nth varable; if null l then 0 else max(nth(caar l,n+1),maxfrom(cdr l,n)); symbolic procedure copy u; if atom u then u else cons(copy car u,copy cdr u); symbolic procedure addctomap cc; begin scalar ncval; ccount:=ccount+1; ncval:=mkvect(ccount); for i:=0:(ccount-1) do putv(ncval,i,getv(cval,i)); putv(ncval,ccount,nil ./ 1); cval:=ncval; cmap:=(cc . ccount).cmap; if !*trint then << prin2 "Constant map changed to "; print cmap >>; cmatrix:=mapcar(cmatrix,function addtovector); end; symbolic procedure addtovector v; begin scalar vv; vv:=mkvect(ccount); for i:=0:(ccount-1) do putv(vv,i,getv(v,i)); putv(vv,ccount,nil); return vv end; symbolic procedure checkcoeffts(cl,indv); % checks to see that the coefficients in CL (coefficient list - S.Q.s); % are zero when the i,j,k,... are given values in INDV (LPOW of; % D.F.). if so the result is true else NIL=false; if null cl then t else begin scalar res; res:=evaluatecoeffts(numr car cl,indv); if not(null res or res=0) then return nil else return checkcoeffts(cdr cl,indv) end; symbolic procedure nextu(ltrhs,powlhs); % picks out the appropriate U coefficients for term: LTRHS to match the % powers of the z-variables given in POWLHS (= exponent list of D.F.). % return this coefficient in form LT of D.F. If U coefficient does % not exist then result is NIL. If it is multiplied by a zero then % result is NIL; if null ltrhs then nil else begin scalar indlist,ucoefft; indlist:=subtractinds(powlhs,car ltrhs,nil); if null indlist then return nil; ucoefft:=evaluatecoeffts(numr cdr ltrhs,indlist); if null ucoefft or ucoefft=0 then return nil; return indlist .* (ucoefft ./ denr cdr ltrhs) end; symbolic procedure subtractinds(powlhs,l,sofar); % subtract the indices in list L from those in POWLHS to find; % appropriate values for i,j,k,... when equating coefficients of terms; % on lhs of reduction eqn. SOFAR is the resulting value list we; % have constructed so far. if any i,j,k,... value is -ve then result; % is NIL; if null l then reversewoc sofar else if ((car powlhs)-(caar l))<0 then nil else subtractinds(cdr powlhs,cdr l, ((car powlhs)-(caar l)) . sofar); symbolic procedure gausselimn(equation,tokill); % Performs Gaussian elimination on the matrix for the c-equations; % as each c-equation is found. EQUATION is the next one to deal with; begin scalar newrow,pivot; if zerop ccount then go to noway; %failure newrow:=mkvect(ccount); spreadc(equation,newrow,1); subst4eliminatedcs(newrow,reverse orderofelim,reverse cmatrix); pivot:=findpivot newrow; if null pivot then go to nopivotfound; orderofelim:=pivot . orderofelim; newrow:=makeprim newrow; %remove hcf from new equation cmatrix:=newrow . cmatrix; % if !*trint then printspreadc newrow; return t; nopivotfound: if null getv(newrow,0) then << if !*trint then printc "Already included"; return nil>>; %equation was 0=0 noway: badpart:=tokill . badpart; %non-integrable term. if !*trint then printc "Inconsistent"; return nil end; symbolic procedure makeprim row; begin scalar g; g:=getv(row,0); for i:=1:ccount do g:=gcdf(g,getv(row,i)); if g neq 1 then for i:=0:ccount do putv(row,i,quotf(getv(row,i),g)); for i := 0:ccount do <<g := getv(row,i); if g and not domainp g then putv(row,i,numr resimp((rootextractf g) ./ 1))>>; return row end; endmodule; module sqrtf; % Square root of standard forms. % Authors: Mary Ann Moore and Arthur C. Norman. fluid '(!*noextend zlist); exports minusdfp,sqrtdf,nrootn,domainp,minusf; imports contentsmv,gcdf,interr,!*multf,partialdiff,printdf,quotf, simpsqrt2,vp2; symbolic procedure minusdfp a; % Test sign of leading coedd of d.f. if null a then interr "Minusdfp 0 illegal" else minusf numr lc a; symbolic procedure sqrtdf l; % Takes square root of distributive form. "Failed" usually means % that the square root is not among already existing objects. if null l then nil else begin scalar c; if lpow l=vp2 zlist then go to ok; c:=sqrtsq df2q l; if numr c eq 'failed then return 'failed; if denr c eq 'failed then return 'failed; return for each u in f2df numr c collect (car u).!*multsq(cdr u,1 ./ denr c); ok: c:=sqrtsq lc l; if numr c eq 'failed or denr c eq 'failed then return 'failed else return (lpow l .* c) .+nil end; symbolic procedure sqrtsq a; sqrtf numr a ./ sqrtf denr a; symbolic procedure sqrtf p; begin scalar ip,qp; if null p then return nil; ip:=sqrtf1 p; qp:=cdr ip; ip:=car ip; %respectable and nasty parts of the sqrt. if qp=1 then return ip; %exact root found. if !*noextend then return 'failed; % We cannot add new square roots in this case, since it is % then impossible to determine if one square root depends % on another if new ones are being added all the time. if zlistp qp then return 'failed; % Liouville's theorem tells you that you never need to add % new algebraics depending on the variable of integration. qp:=simpsqrt2 qp; return !*multf(ip,qp) end; symbolic procedure zlistp qp; if atom qp then member(qp,zlist) else or(member(mvar qp,zlist),zlistp lc qp,zlistp red qp); symbolic procedure sqrtf1 p; % Returns a . b with p=a**2*b. if domainp p then if fixp p then nrootn(p,2) else !*q2f simpsqrt list prepf p . 1 else begin scalar co,pp,g,pg; co:=contentsmv(p,mvar p,nil); %contents of p. pp:=quotf(p,co); %primitive part. co:=sqrtf1(co); %process contents via recursion. g:=gcdf(pp,partialdiff(pp,mvar pp)); pg:=quotf(pp,g); g:=gcdf(g,pg); %a repeated factor of pp. if g=1 then pg:=1 . pp else << pg:= quotf(pp,!*multf(g,g)); %what is still left. pg:=sqrtf1(pg); %split that up. rplaca(pg,!*multf(car pg,g))>>; %put in the thing found here. rplaca(pg, !*multf(car pg,car co)); rplacd(pg, !*multf(cdr pg,cdr co)); return pg end; % NROOTN removed as in REDUCE base. endmodule; module tdiff; % Differentiation routines for integrator. % Authors: Mary Ann Moore and Arthur C. Norman. exports !-!-simpdf; imports simpcar,kernp,diffsq,prepsq,msgpri; flag('(!-!-simpdf),'lose); 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,tt; tt := time(); %start the clock; v := cdr u; u := simpcar u; a: if null v or null numr u then go to exit; x := if null y or y=0 then simpcar v else y; if null kernp x then go to e; x := caaaar x; v := cdr v; if null v then go to c; y := simpcar v; if null numr y then go to d else if not denr y=1 or not numberp numr y then go to c; y := car y; v := cdr v; b: if y=0 then go to a; u := diffsq(u,x); y := y-1; go to b; c: u := diffsq(u,x); go to a; d: y := nil; v := cdr v; go to a; exit: print list('time,time()-tt); return u; e: msgpri("Differentiation wrt",prepsq x,"not allowed",nil,t) end; put('tdf,'simpfn,'!-!-simpdf); endmodule; module tidysqrt; % General tidying up of square roots. % Authors: Mary Ann Moore and Arthur C. Norman. % Modifications by J.H. Davenport. exports sqrt2top,tidysqrt; %symbolic procedure tidysqrtdf a; % if null a then nil % else begin scalar tt,r; % tt:=tidysqrt lc a; % r:=tidysqrtdf red a; % if null numr tt then return r; % return ((lpow a) .* tt) .+ r % end; % symbolic procedure tidysqrt q; begin scalar nn,dd; nn:=tidysqrtf numr q; if null nn then return nil ./ 1; %answer is zero. dd:=tidysqrtf denr q; return multsq(nn,invsq dd) end; symbolic procedure tidysqrtf p; %Input - standard form. %Output - standard quotient. %Simplifies sqrt(a)**n with n>1. if domainp p then p ./ 1 else begin scalar v,w; v:=lpow p; if car v='i then v:=mksp('(sqrt -1),cdr v); %I->sqrt(-1); if eqcar(car v,'sqrt) and not onep cdr v then begin scalar x; %here we have a reduction to apply. x:=divide(cdr v,2); %halve exponent. w:=exptsq(simp cadar v,car x); %rational part of answer. if not zerop cdr x then w:=multsq(w, ((mksp(car v,1) .* 1) .+ nil) ./ 1); %the next line allows for the horrors of nested sqrts. w:=tidysqrt w end else w:=((v .* 1) .+ nil) ./ 1; v:=multsq(w,tidysqrtf lc p); return addsq(v,tidysqrtf red p) end; symbolic procedure multoutdenr q; % Move sqrts in a sq to the numerator. begin scalar n,d,root,conj; n:=numr q; d:=denr q; while (root:=findsquareroot d) do << conj:=conjugatewrt(d,root); n:=!*multf(n,conj); d:=!*multf(d,conj) >>; while (root:=findnthroot d) do << conj:=conjugateexpt(d,root,kord!*); n:=!*multf(n,conj); d:=!*multf(d,conj) >>; return (n . d); end; symbolic procedure conjugateexpt(d,root,kord!*); begin scalar ord,ans,repl,xi; ord:=caddr caddr root; % the denominator of the exponent; ans:=1; kord!*:= (xi:=gensym()) . kord!*; % XI is an ORD'th root of unity; for i:=1:ord-1 do << ans:=!*multf(ans,numr subf(d, list(root . list('times,root,list('explt,xi,i))))); while (mvar ans eq xi) and ldeg ans > ord do ans:=addf(red ans,(xi) to (ldeg ans - ord) .* lc ans .+ nil); if (mvar ans eq xi) and ldeg ans = ord then ans:=addf(red ans,lc ans) >>; if (mvar ans eq xi) and ldeg ans = ord-1 then << repl:=-1; for i:=1:ord-2 do repl:=(xi) to i .* -1 .+ repl; ans:=addf(red ans,!*multf(lc ans,repl)) >>; if not domainp ans and mvar ans eq xi then interr "Conjugation failure"; return ans; end; symbolic procedure sqrt2top q; begin scalar n,d; n:=multoutdenr q; d:=denr n; n:=numr n; if d eq denr q then return q;%no change. if d iequal 1 then return (n ./ 1); q:=gcdcoeffsofsqrts n; if q iequal 1 then if minusf d then return (negf n ./ negf d) else return (n ./ d); q:=gcdf(q,d); n:=quotf(n,q); d:=quotf(d,q); if minusf d then return (negf n ./ negf d) else return (n ./ d) end; %symbolic procedure denrsqrt2top q; %begin % scalar n,d; % n:=multoutdenr q; % d:=denr n; % n:=numr n; % if d eq denr q % then return d; % no changes; % if d iequal 1 % then return 1; % q:=gcdcoeffsofsqrts n; % if q iequal 1 % then return d; % q:=gcdf(q,d); % if q iequal 1 % then return d % else return quotf(d,q) % end; symbolic procedure findsquareroot p; % Locate a sqrt symbol in poly p. if domainp p then nil else begin scalar w; w:=mvar p; %check main var first. if atom w then return nil; %we have passed all sqrts. if eqcar(w,'sqrt) then return w; w:=findsquareroot lc p; if null w then w:=findsquareroot red p; return w end; symbolic procedure findnthroot p; nil; % Until corrected. symbolic procedure x!-findnthroot p; % Locate an n-th root symbol in poly p. if domainp p then nil else begin scalar w; w:=mvar p; %check main var first. if atom w then return nil; %we have passed all sqrts. if eqcar(w,'expt) and eqcar(caddr w,'quotient) then return w; w:=findnthroot lc p; if null w then w:=findnthroot red p; return w end; symbolic procedure conjugatewrt(p,var); % Var -> -var in form p. if domainp p then p else if mvar p=var then begin scalar x,c,r; x:=tdeg lt p; %degree c:=lc p; %coefficient r:=red p; %reductum x:=remainder(x,2); %now just 0 or 1. if x=1 then c:=negf c; %-coefficient. return (lpow p .* c) .+ conjugatewrt(r,var) end else if ordop(var,mvar p) then p else (lpow p .* conjugatewrt(lc p,var)) .+ conjugatewrt(red p,var); symbolic procedure gcdcoeffsofsqrts u; if atom u then if numberp u and minusp u then -u else u else if eqcar(mvar u,'sqrt) then begin scalar v; v:=gcdcoeffsofsqrts lc u; if v iequal 1 then return v else return gcdf(v,gcdcoeffsofsqrts red u) end else begin scalar root; root:=findsquareroot u; if null root then return u; u:=makemainvar(u,root); root:=gcdcoeffsofsqrts lc u; if root iequal 1 then return 1 else return gcdf(root,gcdcoeffsofsqrts red u) end; endmodule; module trcase; % Driving routine for integration of transcendental fns. % Authors: Mary Ann Moore and Arthur C. Norman. % Modifications by: John P. Fitch. fluid '(!*backtrace !*nowarnings !*purerisch !*reverse badpart ccount cmap cmatrix content cuberootflag cval denbad denominator indexlist lhs!* loglist lorder orderofelim rhs!* sillieslist sqfr sqrtflag sqrtlist tanlist svar varlist xlogs zlist); % !*reverse: flag to re-order zlist. % !*nowarnings: flag to lose messages. global '(!*failhard !*number!* !*ratintspecial !*seplogs !*spsize!* !*statistics !*trint 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,prepsq; % 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, sqrtlist,tanlist,loglist,dflogs,eprim,dfun,unintegrand, sqrtflag,badpart,rhs!*,lhs!*,gcdq,cmap,cval,orderofelim,cmatrix; scalar cuberootflag,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); 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)); if !*trint then << printc "Unnormalized integrand ="; printsq unintegrand >> end; divlist:=findtrialdivs zlist; %also puts some things on loglist sometimes. % if !*trint % then << printc "Exponentials and tans to try dividing:"; % print divlist>>; sqrtlist:=findsqrts zlist; % if !*trint then << printc "Square-root z-variables"; % print sqrtlist >>; divlist:=trialdiv(denr unintegrand,divlist); % if !*trint then << printc "Divisors:"; % print car divlist; % print cdr 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:"; % superprint sqfr>>; if !*reverse then zlist:=reverse zlist; %ALTER ORDER FUNCTION; indexlist:=createindices zlist; % if !*trint then << printc "...indices are:"; % superprint indexlist>>; dfu:=dfnumr(svar,car divlist); % if !*trint then << terpri(); % printc "************ Derivative of u is:"; % printsq dfu>>; loglist:=append(loglist,factorlistlist (prim,nil)); 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 interr "Simplification failure"; if !*trint then << printc "Distributed Form of U is:"; printdf rhs!*>>; lhs!*:=f2df dflogs; if checkdffail(lhs!*,svar) then interr "Simplification failure"; if !*trint then << printc "Distributed Form of l.h.s. is:"; printdf lhs!*; terpri()>>; cval:=mkvect(ccount); for i:=0 : ccount do putv(cval,i,nil ./ 1); lorder:=maxorder(rhs!*,zlist,0); originalorder:=lorder; if !*trint then << printc "Maximum order 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 not (ccount=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; % if !*trint then superprint dfun; result:= !*multsq(dfun,!*invsq(denominator ./ 1)); result:= !*multsq(result,!*invsq(jhd!-content ./ 1)); % if !*trint then superprint result; 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 << superprint result; terpri(); printc "*****************************************************"; printc "************ THE INTEGRAL IS : **********************"; printc "*****************************************************"; terpri(); printsq result; terpri()>>; if not null badpart then << if !*trint then printc "plus a bad part"; lhs!*:=badpart; lorder:=maxorder(rhs!*,zlist,0); while lorder do << if car lorder > car originalorder then wrongway:=t; lorder:=cdr lorder; originalorder:=cdr originalorder >>; dfun:=df2q badpart; if !*trint then <<printsq dfun; printc "Denbad = "; printsf denbad>>; dfun:= !*multsq(dfun,invsq(denbad ./ 1)); if wrongway then << result:= nil ./ 1; dfun:=integrand >>; if rootcheckp(unintegrand,svar) then return simpint1(integrand . svar.nil) else if !*purerisch or allowedfns zlist then dfun:=simpint1 (dfun . svar.nil) else << !*purerisch:=t; if !*trint then <<printc " [Transforming ..."; printsq dfun>>; denbad:=transform(dfun,svar); if denbad=dfun then dfun:=simpint1(dfun . svar.nil) else <<denbad:=errorset('(integratesq denbad svar xlogs), !*backtrace,!*backtrace); if not atom denbad then dfun:=untan car denbad else dfun:=simpint1(dfun . svar.nil) >> >>; if !*trint then printsq dfun; if !*failhard then rederr "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) >>; % if !*overlaymode then excise transcode; return sqrt2top result end; symbolic procedure checkdffail(u,v); u and (depends(lc u,v) or 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 "; print i>>; 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; module halfangle; % Routines for conversion to half angle tangents. % Author: Steve Harrington. % Modifications by: John P. Fitch. fluid '(!*gcd); exports halfangle,untan; symbolic procedure transform(u,x); % Transform the SQ U to remove the 'bad' functions sin, cos, cot etc % in favor of half angles; halfangle(u,x); symbolic procedure quotqq(u1,v1); multsq(u1, invsq(v1)); symbolic procedure !*subtrq(u1,v1); addsq(u1, negsq(v1)); symbolic procedure !*int2qm(u1); if u1=0 then nil . 1 else u1 . 1; symbolic procedure halfangle(r,x); % Top level procedure for converting; % R is a rational expression to be converted, % X the integration variable. % A rational expression is returned. quotqq(hfaglf(numr(r),x), hfaglf(denr(r),x)); symbolic procedure hfaglf(p,x); % Converting polynomials, a rational expression is returned. if domainp(p) then !*f2q(p) else subs2q addsq(multsq(exptsq(hfaglk(mvar(p),x), ldeg(p)), hfaglf(lc(p),x)), hfaglf(red(p),x)); symbolic procedure hfaglk(k,x); % Converting kernels, a rational expression is returned. begin scalar kt; if atom k or not member(x,flatten(cdr(k))) then return !*k2q k; k := car(k) . hfaglargs(cdr(k), x); kt := simp list('tan, list('quotient, cadr(k), 2)); return if car(k) = 'sin then quotqq(multsq(!*int2qm(2),kt), addsq(!*int2qm(1), exptsq(kt,2))) else if car(k) = 'cos then quotqq(!*subtrq(!*int2qm(1),exptsq(kt,2)),addsq(!*int2qm(1), exptsq(kt,2))) else if car(k) = 'tan then quotqq(multsq(!*int2qm(2),kt), !*subtrq(!*int2qm(1), exptsq(kt,2))) else if car(k) = 'sinh then quotqq(!*subtrq(exptsq(!*k2q('expt.('e. cdr k)),2), !*int2qm(1)), multsq(!*int2qm(2), !*k2q('expt . ('e . cdr(k))))) else if car(k) = 'cosh then quotqq(addsq(exptsq(!*k2q('expt.('e. cdr k)),2), !*int2qm(1)), multsq(!*int2qm(2), !*k2q('expt . ('e . cdr(k))))) else if car(k) = 'tanh then quotqq(!*subtrq(exptsq(!*k2q('expt.('e. cdr k)),2), !*int2qm(1)), addsq(exptsq(!*k2q ('expt.('e.cdr(k))),2), !*int2qm(1))) else !*k2q(k); % additional transformation might be added here. end; symbolic procedure hfaglargs(l,x); % Conversion of argument list. if null l then nil else prepsq(hfaglk(car(l),x)) . hfaglargs(cdr(l), x); symbolic procedure untanf x; % This should be done by a table. begin scalar y,z,w; if domainp x then return x . 1; y := mvar x; if eqcar(y,'int) then error1(); % assume all is hopeless. z := ldeg x; w := 1 . 1; y := if atom y then !*k2q y else if car y eq 'tan then if evenp z then <<z := z/2; simp list('quotient, list('plus, list('minus, list('cos, 'times . (2 . cdr y))), 1),list('plus, list('cos, 'times . (2 . cdr y)), 1))>> else if z=1 then simp list('quotient, list('plus, list('minus, list('cos, 'times . (2 . cdr y))), 1),list('sin, 'times . (2 . cdr y))) else <<z := (z - 1)/2; w := simp list('quotient, list('plus, list('minus, list('cos, 'times . (2 . cdr y))), 1),list('sin, 'times . (2 . cdr y))); simp list('quotient, list('plus, list('minus, list('cos, 'times . (2 . cdr y))), 1),list('plus, list('cos, 'times . (2 . cdr y)), 1))>> else simp y; return addsq(multsq(multsq(exptsq(y,z),untanf lc x),w), untanf red x) end; symbolic procedure untanlist(y); if null y then nil else (prepsq (untan(simp car y)) . untanlist(cdr y)); symbolic procedure untan(x); % Expects x to be canonical quotient. begin scalar y; y:=cossqchk sinsqrdchk multsq(untanf(numr x), invsq untanf(denr x)); return if length flatten y>length flatten x then x else y end; symbolic procedure sinsqrdchk(x); multsq(sinsqchkf(numr x), invsq sinsqchkf(denr x)); symbolic procedure sinsqchkf(x); begin scalar y,z,w; if domainp x then return x . 1; y := mvar x; z := ldeg x; w := 1 . 1; y := if eqcar(y,'sin) then if evenp z then <<z := quotient(z,2); simp list('plus,1,list('minus, list('expt,('cos . cdr(y)),2)))>> else if z = 1 then !*k2q y else << z := quotient(difference(z,1),2); w := !*k2q y; simp list('plus,1,list('minus, list('expt,('cos . cdr(y)),2)))>> else !*k2q y; return addsq(multsq(multsq(exptsq(y,z),sinsqchkf(lc x)),w), sinsqchkf(red x)); end; symbolic procedure cossqchkf(x); begin scalar y,z,w,x1,x2; if domainp x then return x . 1; y := mvar x; z := ldeg x; w := 1 . 1; x1 := cossqchkf(lc x); x2 := cossqchkf(red x); x := addsq(multsq(!*p2q lpow x,x1),x2); y := if eqcar(y,'cos) then if evenp z then <<z := quotient(z,2); simp list('plus,1,list('minus, list('expt,('sin . cdr(y)),2)))>> else if z = 1 then !*k2q y else << z := quotient(difference(z,1),2); w := !*k2q y; simp list('plus,1,list('minus, list('expt,('sin . cdr(y)),2)))>> else !*k2q y; y := addsq(multsq(multsq(exptsq(y,z),w),x1),x2); return if length(y) > length(x) then x else y; end; symbolic procedure cossqchk(x); begin scalar !*gcd; !*gcd := t; return multsq(cossqchkf(numr x), invsq cossqchkf(denr x)) end; symbolic procedure lrootchk(l,x); % Checks each member of list l for a root. if null l then nil else krootchk(car l, x) or lrootchk(cdr l, x); symbolic procedure krootchk(f,x); % Checks a kernel to see if it is a root. if atom f then nil else if car(f) = 'sqrt and member(x, flatten cdr f) then t else if car(f) = 'expt and not atom caddr(f) and caaddr(f) = 'quotient and member(x, flatten cadr f) then t else lrootchk(cdr f, x); symbolic procedure rootchk1p(f,x); % Checks polynomial for a root. if domainp f then nil else krootchk(mvar f,x) or rootchk1p(lc f,x) or rootchk1p(red f,x); symbolic procedure rootcheckp(f,x); % Checks rational (standard quotient) for a root. rootchk1p(numr f,x) or rootchk1p(denr f,x); endmodule; module trialdiv; % Trial division routines. % Authors: Mary Ann Moore and Arthur C. Norman. fluid '(denominator loglist tanlist); global '(!*trint); exports countz,findsqrts,findtrialdivs,trialdiv,simp,mksp; imports !*multf,printsf,quotf; symbolic procedure countz dl; % DL is a list of S.F.s; begin scalar s,n,rl; loop2: if null dl then return arrangelistz rl; n:=1; loop1: n:=n+1; s:=car dl; dl:=cdr dl; if not null dl and (s eq car dl) then go to loop1 else rl:=(s.n).rl; go to loop2 end; symbolic procedure arrangelistz d; begin scalar n,s,rl,r; n:=1; if null d then return rl; loopd: if (cdar d)=n then s:=(caar d).s else r:=(car d).r; d:=cdr d; if not null d then go to loopd; d:=r; rl:=s.rl; s:=nil; r:=nil; n:=n+1; if not null d then go to loopd; return reversewoc rl end; symbolic procedure findtrialdivs zl; %zl is list of kernels found in integrand. result is a list %giving things to be treated specially in the integration %viz: exps and tans. %Result is list of form ((a . b) ...) % with a a kernel and car a=expt or tan % and b a standard form for either expt or (1+tan**2). begin scalar dlists1,args1; while not null zl do << if exportan car zl then << if caar zl='tan then << args1:=(mksp(car zl,2) .* 1) .+ 1; tanlist:=(args1 ./ 1) . tanlist>> else args1:=!*k2f car zl; dlists1:=(car zl . args1) . dlists1>>; zl:=cdr zl >>; return dlists1 end; symbolic procedure exportan dl; if atom dl then nil else begin % extract exp or tan fns from the z-list. if eq(car dl,'tan) then return t; nxt: if not eq(car dl,'expt) then return nil; dl:=cadr dl; if atom dl then return t; % Make sure we find nested exponentials? go to nxt end; symbolic procedure findsqrts z; begin scalar r; while not null z do << if eqcar(car z,'sqrt) then r:=(car z) . r; z:=cdr z >>; return r end; symbolic procedure trialdiv(x,dl); begin scalar qlist,q; while not null dl do if not null(q:=quotf(x,cdar dl)) then << if (caaar dl='tan) and not eqcar(qlist,cdar dl) then loglist:=('iden . simp cadr caar dl) . loglist; %tan fiddle! qlist:=(cdar dl).qlist; x:=q >> else dl:=cdr dl; return qlist.x end; endmodule; module unifac; % Univariate factorization for integration. % Authors: Mary Ann Moore and Arthur C. Norman. fluid '(knowndiscrimsign zlist); global '(!*trint); exports unifac; imports cubic,linfac,printdf,quadfac,quadratic,quartic,vp1, gcd,minusp,prettyprint; symbolic procedure unifac(pol,var,degree,res); begin scalar w,q,c; w:=pol; if !*trint then superprint w; %now try looking for linear factors. trylin: q:=linfac(w); if null car q then go to nomorelin; res := ('log . back2df(car q,var)) . res; w:=cdr q; go to trylin; nomorelin: q:=quadfac(w); if null car q then go to nomorequad; res := quadratic(back2df(car q,var),var,res); w:=cdr q; go to nomorelin; nomorequad: if null w then return res; %all done. degree:=car w; %degree of what is left. c:=back2df(w,var); if degree=3 then res:=cubic(c,var,res) else if degree=4 then res:=quartic(c,var,res) else if evenp degree and pairp (q := halfpower cddr w) then <<w := (degree/2) . (cadr w . q); w := unifac(w,var,car w,nil); res := pluckfactors(w,var,res)>> else << printc "The following has not been split"; printdf c; res:=('log . c) . res>>; return res end; symbolic procedure halfpower w; if null w then nil else if car w=0 then (lambda r; if r eq 'failed then r else cadr w . r) halfpower cddr w else 'failed; symbolic procedure pluckfactors(w,var,res); begin scalar p,q,knowndiscrimsign; while w do <<p := car w; if car p eq 'atan then nil else if car p eq 'log then <<q := doublepower cdr p . q; %prin2 "q="; %printdf car q; >> else interr "Bad form"; w := cdr w>>; while q do <<p := car q; if caaar p=4 then <<knowndiscrimsign := 'negative; res := quartic(p,var,res); knowndiscrimsign := nil>> else if caaar p=2 then res := quadratic(p,var,res) else res := ('log . p) . res; q := cdr q>>; return res end; symbolic procedure doublepower r; if null r then nil else ((for each j in caar r collect 2*j) . cdar r) . doublepower cdr r; symbolic procedure back2df(p,v); % Undo the effect of uniform. begin scalar r,n; n:=car p; p:=cdr p; while not minusp n do << if not zerop car p then r:= (vp1(v,n,zlist) .* (car p ./ 1)) .+ r; p:=cdr p; n:=n-1 >>; return reversewoc r end; endmodule; module uniform; % Convert from D.F. to list of coefficients. % Authors: Mary Ann Moore and Arthur C. Norman. fluid '(zlist); exports uniform; imports exponentof; symbolic procedure uniform(p,v); %Convert from d.f. in one variable (v) to a simple list of %coeffs (with degree consed onto front). %Fails if coefficients are not all simple integers. if null p then 0 . (0 . nil) else begin scalar a,b,c,d; a:=exponentof(v,lpow p,zlist); b:=lc p; if not(denr b=1) then return 'failed; b:=numr b; if null b then b:=0 else if not numberp b then return 'failed; if a=0 then return a . (b . nil); %constant term. c:=uniform(red p,v); if c='failed then return 'failed; d:=car c; c:=cdr c; d:=d+1; while not (a=d) do << c:=0 . c; d:=d+1>>; return a . (b . c) end; endmodule; module makevars; % Make dummy variables for integration process. % Authors: Mary Ann Moore and Arthur C. Norman. fluid '(!*gensymlist!* !*purerisch); exports getvariables,varsinlist,varsinsq,varsinsf,findzvars, createindices,mergein; imports dependsp,union; % Note that 'i' is already maybe committed for sqrt(-1), % also 'l' and 'o' are not used as they print badly on certain % terminals etc and may lead to confusion. !*gensymlist!* := '(! j ! k ! m ! n ! p ! q ! r ! s ! t ! u ! v ! w ! x ! y ! z); %mapc(!*gensymlist!*,function remob); %REMOB protection; symbolic procedure varsinlist(l,vl); % L is a list of s.q. - find all variables mentioned, % given thal vl is a list already known about. begin while not null l do << vl:=varsinsf(numr car l,varsinsf(denr car l,vl)); l:=cdr l >>; return vl end; symbolic procedure getvariables sq; varsinsf(numr sq,varsinsf(denr sq,nil)); symbolic procedure varsinsq(sq,vl); varsinsf(numr sq,varsinsf(denr sq,vl)); symbolic procedure varsinsf(form,l); if domainp form then l else begin while not domainp form do << l:=varsinsf(lc form,union(l,list mvar form)); form:=red form >>; return l end; symbolic procedure findzvars(vl,zl,var,flg); begin scalar v; % VL is the crude list of variables found in the original integrand; % ZL must have merged into it all EXP, LOG etc terms from this. % If FLG is true then ignore DF as a function. scan: if null vl then return zl; v:=car vl; % next variable. vl:=cdr vl; % at present items get put onto ZL if they are non-atomic % and they depend on the main variable. The arguments of % functions are decomposed by recursive calls to findzvar. %give up if V has been declared dependent on other things; if atom v and v neq var and depends(v,var) then rederr "Can't integrate in the presence of side-relations" else if not atom v and (not v member zl) and dependsp(v,var) then if car v='!*sq then zl:=findzvarssq(cadr v,zl,var) else if car v memq '(times quotient plus minus difference int) or (((car v) eq 'expt) and fixp caddr v) then zl:=findzvars(cdr v,zl,var,flg) else if flg and car v eq 'df then <<!*purerisch := t; return zl>> % try and stop it else zl:=v . findzvars(cdr v,zl,var,flg); % scan arguments of fn. %ACH: old code used to look only at CADR if a DF involved. go to scan end; symbolic procedure findzvarssq(sq,zl,var); findzvarsf(numr sq,findzvarsf(denr sq,zl,var),var); symbolic procedure findzvarsf(sf,zl,var); if domainp sf then zl else findzvarsf(lc sf, findzvarsf(red sf, findzvars(list mvar sf,zl,var,nil), var), var); symbolic procedure createindices zl; % Produces a list of unique indices, each associated with a ; % different Z-variable; reversewoc crindex1(zl,!*gensymlist!*); symbolic procedure crindex1(zl,gl); begin if null zl then return nil; if null gl then << gl:=list int!-gensym1 'i; %new symbol needed; nconc(!*gensymlist!*,gl) >>; return (car gl) . crindex1(cdr zl,cdr gl) end; symbolic procedure rmember(a,b); if null b then nil else if a=cdar b then car b else rmember(a,cdr b); symbolic procedure mergein(dl,ll); % Adjoin logs of things in dl to existing list ll. if null dl then ll else if rmember(car dl,ll) then mergein(cdr dl,ll) else mergein(cdr dl,('log . car dl) . ll); endmodule; module vect; % Vector support routines. % Authors: Mary Ann Moore and Arthur C. Norman. % Modified by: James H. Davenport. exports mkuniquevect,mkvec,mkvecf2q,mkidenm,copyvec,vecsort,swap, non!-null!-vec,mkvect2; symbolic procedure mkuniquevect v; begin scalar u,n; n:=upbv v; for i:=0:n do begin scalar uu; uu:=getv(v,i); if not (uu member u) then u:=uu.u end; return mkvec u end; symbolic procedure mkvec(l); begin scalar v,i; v:=mkvect(isub1 length l); i:=0; while l do <<putv(v,i,car l); i:=iadd1 i; l:=cdr l>>; return v end; symbolic procedure mkvecf2q(l); begin scalar v,i,ll; v:=mkvect(isub1 length l); i:=0; while l do << ll:=car l; if ll = 0 then ll:=nil; putv(v,i,!*f2q ll); i:=iadd1 i; l:=cdr l >>; return v end; symbolic procedure mkidenm n; begin scalar ans,u; scalar c0,c1; c0:=nil ./ 1; c1:= 1 ./ 1; % constants. ans:=mkvect(n); for i:=0 step 1 until n do << u:=mkvect n; for j:=0 step 1 until n do if i iequal j then putv(u,j,c1) else putv(u,j,c0); putv(ans,i,u) >>; return ans end; symbolic procedure copyvec(v,n); begin scalar new; new:=mkvect(n); for i:=0:n do putv(new,i,getv(v,i)); return new end; symbolic procedure vecsort(u,l); % Sorts vector v of numbers into decreasing order. % Performs same interchanges of all vectors in the list l. begin scalar j,k,n,v,w; n:=upbv u;% elements 0...n exist. % algorithm used is a bubble sort. for i:=1:n do begin v:=getv(u,i); k:=i; loop: j:=k; k:=isub1 k; w:=getv(u,k); if v<=w then goto ordered; putv(u,k,v); putv(u,j,w); mapc(l,function (lambda u;swap(u,j,k))); if k>0 then goto loop; ordered: end; return nil end; symbolic procedure swap(u,j,k); if null u then nil else begin scalar v; %swaps elements i,j of vector u. v:=getv(u,j); putv(u,j,getv(u,k)); putv(u,k,v) end; symbolic procedure non!-null!-vec v; begin scalar cnt; cnt := 0; for i:=0:upbv v do if getv(v,i) then cnt:=iadd1 cnt; return cnt end; symbolic procedure mkvect2(n,initial); begin scalar u; u:=mkvect n; for i:=0:n do putv(u,i,initial); return u end; endmodule; end;