123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111 |
- module intfac; % Interface between integrator and factorizer.
- % Author: Anthony C. Hearn.
- % Based on earlier versions by James Davenport, Mary Ann Moore and
- % Arthur Norman.
- fluid '(!*intfac !*surds kord!* zlist); % clogflag
- exports int!-fac;
- symbolic procedure int!-fac x;
- % X is a primitive, square-free polynomial, except for monomial
- % factors. Result is a list of 'factors' wrt zlist.
- % Throughout most of the integrator we want to add new surds, so we
- % turn surds on. However, we use *intfac to inhibit use of quadratic
- % factorizer in the poly/primfac module, since things don't work
- % properly if this is used.
- begin scalar !*intfac,!*surds;
- !*intfac := !*surds := t;
- return int!-fac!-inner x;
- end;
- symbolic procedure int!-fac!-inner x;
- % X is a primitive, square-free polynomial, except for monomial
- % factors. Result is a list of 'factors' wrt zlist.
- begin scalar factors;
- factors := fctrf x;
- factors := cdr factors; % Ignore monomial factor.
- % Make sure x really square-free.
- factors := for each u in factors
- collect if cdr u=1 then car u
- else interr list(x,"not square free");
- % It seems we need the logs ordered ahead of atans, hence reverse.
- return reversip for each u in factors join fac2int u
- end;
- symbolic procedure fac2int u;
- % Returns a list of all the arctangents and logarithms arising from
- % an attempt to take the one irreducible (but not necessarily the
- % absolutely irreducible) factor u.
- begin scalar degrees,x;
- degrees := for each w in zlist collect (degreef(u,w) . w);
- if assoc(1,degrees) then return list ('log . (u ./ 1))
- % An irreducible polynomial of degree 1 is absolutely irreducible.
- else if x := assoc(2,degrees) then return int!-quadterm(u,cdr x)
- else if assoc(0,degrees) then return list('log . (u ./ 1));
- % This suggests a surd occurs. Should that be an error?
- if !*trint
- then <<printc "*** Polynomial";
- printsf u;
- printc "has not been completely factored">>;
- return list ('log . (u ./ 1))
- end;
- symbolic procedure int!-quadterm(pol,var);
- % Add in logs and atans corresponding to splitting the polynomial pol
- % given it is quadratic wrt var. Does not assume pol is univariate.
- % We need to rootxf!* so that
- % int(1/(x**2*y0+x**2+x*y0**2+3*x*y0+x+y0**2+y0),x) comes out in
- % terms of real functions.
- begin scalar a,b,c,discrim,kord,res,w;
- kord := setkorder(var . kord!*); % It shouldn't matter if
- % var occurs twice.
- c := reorder pol;
- if ldeg c neq 2
- then <<setkorder kord;
- rerror(int,5,"Invalid polynomial in int-quadterm")>>;
- a := lc c;
- c := red c;
- if not domainp c and mvar c = var and ldeg c = 1
- then <<b := lc c; c := red c>>;
- setkorder kord;
- discrim := powsubsf addf(multf(b,b),multd(-4,multf(a,c)));
- if null discrim then interr "discrim is zero in quadterm";
- % A quadratic usually implies an atan term.
- % if not clogflag
- % then <<w := rootxf(negf discrim,2);
- % if not(w eq 'failed) then go to atancase>>;
- w := rootxf!*(negf discrim,2);
- if not(w eq 'failed) then go to atancase;
- w := rootxf!*(discrim,2); % Maybe only rootxf is needed here.
- if w eq 'failed then return list ('log . !*f2q pol);
- % if w eq 'failed then rederr "Integration failure in int-quadterm";
- discrim := w;
- w := multpf(mksp(var,1),a);
- w := addf(multd(2,w),b); % 2*a*x + b.
- a := addf(w,discrim);
- b := addf(w,negf discrim);
- % Remove monomial multipliers.
- a := quotf(a,cdr comfac a);
- b := quotf(b,cdr comfac b);
- return ('log . !*f2q a) . ('log . !*f2q b) . res;
- atancase:
- res := ('log . !*f2q pol) . res; % One part of answer.
- a := multpf(mksp(var,1),a);
- a := addf(b,multd(2,a));
- a := fquotf(a,w);
- return ('atan . a) . res
- end;
- symbolic procedure rootxf!*(u,n);
- (if x eq 'failed or smemq('i,x) and not smemq('i,u)
- then (rootxf(u,n) where !*surds=nil)
- else x)
- where x=rootxf(u,n);
- endmodule;
- end;
|