Artifact cb4f5e8eb41c93bdd2ed9b3a93de7930c57ae799289992f79af4d5ccc96dce1f:
- Executable file
r38/packages/matrix/resultnt.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: 8719) [annotate] [blame] [check-ins using] [more...]
module resultnt; % Author: Eberhard Schruefer. % Modifications by: Anthony C. Hearn, Winfried Neun. %********************************************************************** % * % The resultant function defined here has the following properties: * % * % degr(p1,x)*degr(p2,x) * % resultant(p1,p2,x) = (-1) *resultant(p2,p1,x) * % * % degr(p2,x) * % resultant(p1,p2,x) = p1 if p1 free of x * % * % resultant(p1,p2,x) = 1 if p1 free of x and p2 free of x * % * %********************************************************************** %exports resultant; %imports reorder,setkorder,degr,addf,negf,multf,multpf; load_package polydiv; fluid '(!*bezout !*exp kord!*); switch bezout; put('resultant,'simpfn,'simpresultant); symbolic procedure simpresultant u; if length u neq 3 then rerror(matrix,19, "Resultant called with wrong number of arguments") else resultantsq(simp!* car u,simp!* cadr u,!*a2k caddr u) where !*exp = t; symbolic procedure resultant(u,v,var); % Kept for compatibility with old code. if domainp u and domainp v then 1 else begin scalar x; kord!* := var . kord!*; % updkorder can't be used here. % See sum test. if null domainp u and null(mvar u eq var) then u := reorder u; if null domainp v and null(mvar v eq var) then v := reorder v; x := if !*bezout then bezout_resultant(u,v,var) else polyresultantf(u,v,var); setkorder cdr kord!*; return x end; symbolic procedure resultantsq(u,v,var); % Uses resultant(a*P,b*Q,var) = a^ldeg Q*b^ldeg P*resultant(P,Q,var). if domainp numr u and domainp numr v and denr u = 1 and denr v = 1 then 1 ./ 1 else begin scalar x,y,z; kord!* := var . kord!*; % updkorder can't be used here. % See sum test. if null domainp numr u and null(mvar numr u eq var) then u := reordsq u; if null domainp numr v and null(mvar numr v eq var) then v := reordsq v; if (y := denr u) neq 1 and smember(var,y) then typerr(prepf y,'polynomial) else if (z := denr v) neq 1 and smember(var,z) then typerr(prepf z,'polynomial); u := numr u; v := numr v; if smember(var,coefflist(u,var)) then typerr(prepf u,'polynomial) else if smember(var,coefflist(v,var)) then typerr(prepf v,'polynomial); x := if !*bezout then bezout_resultant(u,v,var) else polyresultantf(u,v,var); if y neq 1 then y := exptf(y,degr(v,var)); if z neq 1 then y := multf(y,exptf(z,degr(u,var))); setkorder cdr kord!*; return x ./ y end; symbolic procedure coefflist(u,var); % Returns list of pairs of degrees and coefficients of var in u. begin scalar z; while not domainp u and mvar u=var do <<z := (ldeg u . lc u) . z; u := red u>>; return if null u then z else (0 . u) . z end; symbolic procedure polyresultantf(u,v,var); % Algorithm is from M. Bronstein, Symbolic Integration I - % Transcendental Functions Algorithms and Computation in Mathematics, % Vol. 1 Springer, Heidelberg, ISBN 3-540-60521-5. % Note var is assumed to be the leading term in kord!* and the main % variable in u and v. begin scalar beta,cd,cn,delta,gam,r,s,temp,x; cd := cn := r := s := 1; gam := -1; if domainp u or domainp v then return 1 else if ldeg u<ldeg v then <<s := (-1)^(ldeg u*ldeg v); temp := u; u := v; v := temp>>; while v do <<delta := ldeg u-ldegr(v,var); beta := negf(multf(r,exptf(gam,delta))); r := lcr(v,var); if delta neq 0 then gam := quotf!*(exptf(negf r,delta),exptf(gam,delta-1)); temp := u; u := v; if not evenp ldeg temp and not evenp ldegr(u,var) then s := -s; v := quotf(pseudo_remf(temp,v,var),beta); if v % We assume that u has var as its main variable. then <<cn := multf(cn,exptf(beta,ldeg u)); cd := multf(cd, exptf(r,(1+delta)*ldeg u-ldeg temp+ldegr(v,var))); if (x := quotf(cd,cn)) then <<cn := 1; cd := x>>>>>>; return if not domainp u and mvar u eq var then nil else if ldeg temp neq 1 then quotf(multf(s,multf(cn,exptf(u,ldeg temp))),cd) else u end; symbolic procedure lcr(u,var); if domainp u or mvar u neq var then u else lc u; symbolic procedure ldegr(u,var); if domainp u or mvar u neq var then 0 else ldeg u; symbolic procedure pseudo_remf(u,v,var); !*q2f simp pseudo!-remainder {mk!*sq(u ./ 1),mk!*sq(v ./ 1),var}; symbolic procedure bezout_resultant(u,v,w); % U and v are standard forms. Result is resultant of u and v % w.r.t. kernel w. Method is Bezout's determinant using exterior % multiplication for its calculation. begin integer n,nm; scalar ap,ep,uh,ut,vh,vt; if domainp u or null(mvar u eq w) then return if not domainp v and mvar v eq w then exptf(u,ldeg v) else 1 else if domainp v or null(mvar v eq w) then return if mvar u eq w then exptf(v,ldeg u) else 1; n := ldeg v - ldeg u; if n < 0 then return multd((-1)**(ldeg u*ldeg v), bezout_resultant(v,u,w)); ep := 1; nm := ldeg v; uh := lc u; vh := lc v; ut := if n neq 0 then multpf(w to n,red u) else red u; vt := red v; ap := addf(multf(uh,vt),negf multf(vh,ut)); ep := b!:extmult(!*sf2exb(ap,w),ep); for j := (nm - 1) step -1 until (n + 1) do <<if degr(ut,w) = j then <<uh := addf(lc ut,multf(!*k2f w,uh)); ut := red ut>> else uh := multf(!*k2f w,uh); if degr(vt,w) = j then <<vh := addf(lc vt,multf(!*k2f w,vh)); vt := red vt>> else vh := multf(!*k2f w,vh); ep := b!:extmult(!*sf2exb(addf(multf(uh,vt), negf multf(vh,ut)),w),ep)>>; if n neq 0 then <<ep := b!:extmult(!*sf2exb(u,w),ep); for j := 1:(n-1) do ep := b!:extmult(!*sf2exb(multpf(w to j,u),w),ep)>>; return if null ep then nil else lc ep end; symbolic procedure !*sf2exb(u,v); %distributes s.f. u with respect to powers in v. if degr(u,v)=0 then if null u then nil else list 0 .* u .+ nil else list ldeg u .* lc u .+ !*sf2exb(red u,v); %**** Support for exterior multiplication **** % Data structure is lpow ::= list of degrees in exterior product % lc ::= standard form symbolic procedure b!:extmult(u,v); %Special exterior multiplication routine. Degree of form v is %arbitrary, u is a one-form. if null u or null v then nil else if v = 1 then u else (if x then cdr x .* (if car x then negf multf(lc u,lc v) else multf(lc u,lc v)) .+ b!:extadd(b!:extmult(!*t2f lt u,red v), b!:extmult(red u,v)) else b!:extadd(b!:extmult(red u,v), b!:extmult(!*t2f lt u,red v))) where x = b!:ordexn(car lpow u,lpow v); symbolic procedure b!:extadd(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 x then y else lpow u .* x .+ y) (addf(lc u,lc v),b!:extadd(red u,red v)) else if b!:ordexp(lpow u,lpow v) then lt u .+ b!:extadd(red u,v) else lt v .+ b!:extadd(u,red v); symbolic procedure b!:ordexp(u,v); if null u then t else if car u > car v then t else if car u = car v then b!:ordexp(cdr u,cdr v) else nil; symbolic procedure b!:ordexn(u,v); %u is a single integer, v a list. Returns nil if u is a member %of v or a dotted pair of a permutation indicator and the ordered %list of u merged into v. begin scalar s,x; a: if null v then return(s . reverse(u . x)) else if u = car v then return nil else if u and u > car v then return(s . append(reverse(u . x),v)) else <<x := car v . x; v := cdr v; s := not s>>; go to a end; endmodule; end;