Artifact e7526b01039c09a70270c8883c7f85b432b673451f6409523edf6537a59222e9:
- Executable file
r38/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: 11720) [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. % Modification of the function 'dipprodin' by Arthur Norman (august 2002, % REDUCE 3.7). 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 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); % The following function was observed recursing very deeply indeed when % certain examples were attempted. Automatic recursion to iteration % conversion in the compiler was not applicable in this case, so a hand % adjustment follows. % symbolic procedure dipprodin(p1,p2); % Distributive polynomial product internal. p1 and p2 are distrib % polynomials. dipprodin(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; % This next definition is one that recursion elimination can handle. % As compared to the original code it introduces a slight time % inefficiency. The original version exploited the fact that the leading % monomial in the result was the product of the two input leading % monomials. In this version dipsum will have to do an exponent % comparison to re-discover this. But the assymptotic overhead grows % linearly while the overall cost here grows quadratically (or worse) if % the two input polys are around the same length, so the cost is ok. symbolic procedure dipprodin(p1, p2); % Distributive polynomial product internal. p1 and p2 are distrib % polynomials. dipprodin(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 dipsum(dipprodin1(diplbc p1,dipevlmon p1,p2), dipprodin(dipmred p1,p2)); symbolic procedure dipprodin1(p1lbc,p1lmon,p2); if dipzero!? p2 then dipzero else dipmoncomp(bcprod(p1lbc,diplbc p2), evsum(p1lmon,dipevlmon p2), dipprodin1(p1lbc,p1lmon,dipmred p2)); 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. % Iterative version,better suited for very long polynomials. if null p1 then p2 else if null p2 then p1 else begin scalar al,done,ep1,ep2,nt,rw,sl,w; 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)else 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;