Artifact 854d5e0afb0e096987da949cd8c2310c997a83985ad066967a4dfc2887d6ee8d:
- Executable file
r37/packages/dipoly/dipoly1.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: 12924) [annotate] [blame] [check-ins using] [more...]
module dipoly1; % Distributive polynomial algorithms. % Authors: R. Gebauer, A. C. Hearn, H. Kredel. % Modification for REDUCE > 3.3: H. Melenk. fluid '(dipvars!* dipzero); symbolic procedure dipconst!? p; not dipzero!? p and dipzero!? dipmred p and evzero!? dipevlmon p; symbolic procedure terprit n; for i:=1:n do terpri(); symbolic procedure dfcprint pl; % h polynomial factor list of distributive polynomials print. for each p in pl do dfcprintin p; symbolic procedure dfcprintin p; % factor with exponent print. ( if cdr p neq 1 then << prin2 " ( "; dipprint1(p1,nil); prin2 " )** "; prin2 cdr p; terprit 2 >> else << prin2 " "; dipprint p1>> ) where p1:= dipmonic a2dip prepf car p; symbolic procedure dfcprin p; % print content, factors and exponents of factorized polynomial p. << terpri(); prin2 " content of factorized polynomials = "; prin2 car p; terprit 2; dfcprint cdr p >>; symbolic procedure diplcm p; % Distributive polynomial least common multiple of denominators. % p is a distributive rational polynomial. diplcm(p) calculates % the least common multiple of the denominators and returns a % base coefficient of the form 1/lcm(denom bc1,.....,denom bci). if dipzero!? p then mkbc(1,1) else bclcmd(diplbc p, diplcm dipmred p); symbolic procedure diprectoint(p,u); % Distributive polynomial conversion rational to integral. % p is a distributive rational polynomial, u is a base coefficient % ( 1/lcm denom p ). diprectoint(p,u) returns a primitive % associate pseudo integral ( denominators are 1 ) distributive % polynomial. if bczero!? u then dipzero else if bcone!? u then p else diprectoint1(p,u); symbolic procedure diprectoint1(p,u); % Distributive polynomial conversion rational to integral internal 1. % diprectoint1 is used in diprectoint. if dipzero!? p then dipzero else dipmoncomp(bclcmdprod(u,diplbc p),dipevlmon p, diprectoint1(dipmred p,u)); symbolic procedure dipresul(p1,p2); % test for fast downwards calculation % p1 and p2 are two bivariate distributive polynomials. % dipresul(p1,p2) returns the resultant of p1 and p2 with respect % respect to the first variable of the variable list (car dipvars!*). begin scalar q1,q2,q,ct; q1:=dip2a diprectoint(p1,diplcm p1); q2:=dip2a diprectoint(p2,diplcm p2); ct := time(); q:= a2dip prepsq simpresultant list(q1,q2,car dipvars!*); ct := time() - ct; prin2 " resultant : "; dipprint dipmonic q; terpri(); prin2 " time resultant : "; prin2 ct; terpri(); end; symbolic procedure dipbcprod (p,a); % /* Distributive polynomial base coefficient product. % p is a distributive polynomial, a is a base coefficient. % dipbcprod(p,a) computes p*a, a distributive polynomial. */ if bczero!? a then dipzero else if bcone!? a then p else dipbcprodin(p,a); symbolic procedure dipbcprodin (p,a); % /* Distributive polynomial base coefficient product internal. % p is a distributive polynomial, a is a base coefficient, % where a is not equal 0 and not equal 1. % dipbcprodin(p,a) computes p*a, a distributive polynomial. */ if dipzero!? p then dipzero else dipmoncomp(bcprod(a, diplbc p), dipevlmon p, dipbcprodin(dipmred p, a)); symbolic procedure dipdif (p1,p2); % /* Distributive polynomial difference. p1 and p2 are distributive % polynomials. dipdif(p1,p2) calculates the difference of the % two distributive polynomials p1 and p2, a distributive polynomial*/ if dipzero!? p1 then dipneg p2 else if dipzero!? p2 then p1 else ( if sl = 1 then dipmoncomp(diplbc p1, ep1, dipdif(dipmred p1, p2) ) else if sl = -1 then dipmoncomp(bcneg diplbc p2, ep2, dipdif(p1,dipmred p2)) else ( if bczero!? al then dipdif(dipmred p1, dipmred p2) else dipmoncomp(al, ep1, dipdif(dipmred p1, dipmred p2) ) ) where al = bcdif(diplbc p1, diplbc p2) ) where sl = evcomp(ep1, ep2) where ep1 = dipevlmon p1, ep2 = dipevlmon p2; symbolic procedure diplength p; % /* Distributive polynomial length. p is a distributive % polynomial. diplength(p) returns the number of terms % of the distributive polynomial p, a digit.*/ if dipzero!? p then 0 else 1 + diplength dipmred p; symbolic procedure diplistsum pl; % /* Distributive polynomial list sum. pl is a list of distributive % polynomials. diplistsum(pl) calculates the sum of all polynomials % and returns a list of one distributive polynomial. */ if null pl or null cdr pl then pl else diplistsum(dipsum(car pl, cadr pl) . diplistsum cddr pl); symbolic procedure diplmerge (pl1,pl2); % /* Distributive polynomial list merge. pl1 and pl2 are lists % of distributive polynomials where pl1 and pl2 are in non % decreasing order. diplmerge(pl1,pl2) returns the merged % distributive polynomial list of pl1 and pl2. */ if null pl1 then pl2 else if null pl2 then pl1 else ( if sl >= 0 then cpl1 . diplmerge(cdr pl1, pl2) else cpl2 . diplmerge(cdr pl2, pl1) ) where sl = evcomp(ep1, ep2) where ep1 = dipevlmon cpl1, ep2 = dipevlmon cpl2 where cpl1 = car pl1, cpl2 = car pl2; symbolic procedure diplsort pl; % /* Distributive polynomial list sort. pl is a list of % distributive polynomials. diplsort(pl) returns the % sorted distributive polynomial list of pl. sort(pl, function dipevlcomp); symbolic procedure dipevlcomp (p1,p2); % /* Distributive polynomial exponent vector leading monomial % compare. p1 and p2 are distributive polynomials. % dipevlcomp(p1,p2) returns a boolean expression true if the % distributive polynomial p1 is smaller or equal the distributive % polynomial p2 else false. */ not evcompless!?(dipevlmon p1, dipevlmon p2); symbolic procedure dipmonic p; % /* Distributive polynomial monic. p is a distributive % polynomial. dipmonic(p) computes p/lbc(p) if p is % not equal dipzero and returns a distributive % polynomial, else dipmonic(p) returns dipzero. */ if dipzero!? p then p else dipbcprod(p, bcinv diplbc p); symbolic procedure dipneg p; % /* Distributive polynomial negative. p is a distributive % polynomial. dipneg(p) returns the negative of the distributive % polynomial p, a distributive polynomial. */ if dipzero!? p then p else dipmoncomp ( bcneg diplbc p, dipevlmon p, dipneg dipmred p ); symbolic procedure dipone!? p; % /* Distributive polynomial one. p is a distributive polynomial. % dipone!?(p) returns a boolean value. If p is the distributive % polynomial one then true else false. */ not dipzero!? p and dipzero!? dipmred p and evzero!? dipevlmon p and bcone!? diplbc p; symbolic procedure dippairsort pl; % /* Distributive polynomial list pair merge sort. pl is a list % of distributive polynomials. dippairsort(pl) returns the % list of merged and in non decreasing order sorted % distributive polynomials. */ if null pl or null cdr pl then pl else diplmerge(diplmerge( car(pl) . nil, cadr(pl) . nil ), dippairsort cddr pl); symbolic procedure dipprod (p1,p2); % /* Distributive polynomial product. p1 and p2 are distributive % polynomials. dipprod(p1,p2) calculates the product of the % two distributive polynomials p1 and p2, a distributive polynomial*/ if diplength p1 <= diplength p2 then dipprodin(p1, p2) else dipprodin(p2, p1); symbolic procedure dipprodin (p1,p2); % /* Distributive polynomial product internal. p1 and p2 are distrib % polynomials. dipprod(p1,p2) calculates the product of the % two distributive polynomials p1 and p2, a distributive polynomial*/ if dipzero!? p1 or dipzero!? p2 then dipzero else ( dipmoncomp(bcprod(bp1, diplbc p2), evsum(ep1, dipevlmon p2), dipsum(dipprodin(dipfmon(bp1, ep1), dipmred p2), dipprodin(dipmred p1, p2) ) ) ) where bp1 = diplbc p1, ep1 = dipevlmon p1; symbolic procedure dipprodls (p1,p2); % /* Distributive polynomial product. p1 and p2 are distributive % polynomials. dipprod(p1,p2) calculates the product of the % two distributive polynomials p1 and p2, a distributive polynomial % using distributive polynomials list sum (diplistsum). */ if dipzero!? p1 or dipzero!? p2 then dipzero else car diplistsum if diplength p1 <= diplength p2 then dipprodlsin(p1, p2) else dipprodlsin(p2, p1); symbolic procedure dipprodlsin (p1,p2); % /* Distributive polynomial product. p1 and p2 are distributive % polynomials. dipprod(p1,p2) calculates the product of the % two distributive polynomials p1 and p2, a distributive polynomial % using distributive polynomials list sum (diplistsum). */ if dipzero!? p1 or dipzero!? p2 then nil else ( dipmoncomp(bcprod(bp1, diplbc p2), evsum(ep1, dipevlmon p2), car dipprodlsin(dipfmon(bp1, ep1), dipmred p2)) . dipprodlsin(dipmred p1, p2) ) where bp1 = diplbc p1, ep1 = dipevlmon p1; %symbolic procedure dipsum (p1,p2); % /* Distributive polynomial sum. p1 and p2 are distributive % polynomials. dipsum(p1,p2) calculates the sum of the % two distributive polynomials p1 and p2, a distributive polynomial*/ % % if dipzero!? p1 then p2 % else if dipzero!? p2 then p1 % else ( if sl = 1 then dipmoncomp(diplbc p1, % ep1, % dipsum(dipmred p1, p2) ) % else if sl = -1 then dipmoncomp(diplbc p2, % ep2, % dipsum(p1,dipmred p2)) % else ( if bczero!? al then dipsum(dipmred p1, % dipmred p2) % else dipmoncomp(al, % ep1, % dipsum(dipmred p1, % dipmred p2) ) % ) where al = bcsum(diplbc p1, diplbc p2) % ) where sl = evcomp(ep1, ep2) % where ep1 = dipevlmon p1, ep2 = dipevlmon p2; symbolic procedure dipsum (p1,p2); % Distributive polynomial sum. p1 and p2 are distributive % polynomials. dipsum(p1,p2) calculates the sum of the % two distributive polynomials p1 and p2. % Iterative version, better suited for very long polynomials. % Warning: this routine uses "dipmred" == "cdr cdr" for a % destructive concatenation. begin scalar w,rw,sl,ep1,ep2,nt,al,done; while not done do <<if dipzero!? p1 then <<nt:=p2; done:=t>> else if dipzero!? p2 then <<nt:=p1; done:=t>> else <<ep1 := dipevlmon p1; ep2 := dipevlmon p2; sl := evcomp(ep1, ep2); % Compute the next term. if sl #= 1 then <<nt := dipmoncomp(diplbc p1, ep1, nil); p1 := dipmred p1>> else if sl #= -1 then <<nt := dipmoncomp(diplbc p2, ep2, nil); p2 := dipmred p2>> else <<al := bcsum(diplbc p1, diplbc p2); nt := if not bczero!? al then dipmoncomp(al,ep1,nil); p1 := dipmred p1; p2 := dipmred p2>>; >>; % Append the term to the sum polynomial. if nt then if null w then w:=rw:=nt else <<cdr cdr rw := nt; rw := nt>>; >>; return w; end; endmodule; end;