Artifact 6be40125bf650255e362718a7f4da94bd393b410e4ea50afb326c1591636e393:
- Executable file
r37/packages/sum/prod.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: 9391) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/sum/prod.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: 9391) [annotate] [blame] [check-ins using]
module prod; % Module for production of finite terms. % Author: F.Kako Hiroshima Univ. % Mon Sep. 7, 1987. % % usage: % prod(expression,variable[,lower[,upper]]); % lower and upper are optional. fluid '(!*trsum prod_last_attempt_rules!*); symbolic procedure simp!-prod u; %ARGUMENT CAR U: expression of prefix form; % CADR U: kernel; % CADDR U: lower bound; % CADDDR U: upper bound; %value : expression of sq form; begin scalar v,y,upper,lower,lower1,dif; y := cdr u; u := simp!* car u; if null numr u then return (1 ./ 1) else if atom y then return u; if not atom cdr y then << lower := cadr y; lower1 := if numberp lower then lower - 1 else list('plus,lower,-1); upper := if not atom cddr y then caddr y else car y; dif := addsq(simp!* upper, negsq simp!* lower); if denr dif = 1 then if null numr dif then return subsq(u,list(!*a2k car y . upper)) else if fixp numr dif then dif := numr dif else dif := nil else dif := nil; if dif and dif <= 0 then return 1 ./ 1; if atom cddr y then upper := nil>>; v := !*a2k car y; return simp!-prod1(u,v,y,upper,lower,lower1,dif) end; symbolic procedure simp!-prod1(u,v,y,upper,lower,lower1,dif); begin scalar w,lst,x,z,flg; lst := prod!-split!-exp(u,v); w := car lst; lst := cdr lst; u := 1 ./ 1; a: if null w then go to b; x := simp!-sum1(cdar w,v,y,upper,lower,lower1,dif); u := multsq(u,simpexpt list(caar w, prepsq x)); w := cdr w; go to a; b: if null lst then return u; flg := nil; z := car lst; if !*trsum then << prin2!* "Product ";sqprint z;prin2!* " w.r.t "; xprinf(!*k2f v,nil,nil);terpri!* t >>; % z := reorder numr z ./ reorder denr z; w := prod!-sq(z,v); if w = 'failed then << if !*trsum then << prin2!* "PROD-SQ failed."; terpri!* t>>; w := prod!-unknown(z,v,y,lower,dif); flg := car w; w := cdr w>>; if !*trsum then << prin2!* "Result = "; sqprint w;terpri!* t >>; if flg then goto c; if upper then w := multsq(subsq(w,list(v . upper)), invsq subsq(w,list(v . lower1))) else if lower then w := multsq(w , invsq subsq(w, list(v . lower1))); c: u := multsq(u,w); lst := cdr lst; goto b end; put('prod,'simpfn,'simp!-prod); %********************************************************************* % Case of other functions %********************************************************************; symbolic procedure prod!-unknown(u,v,y,lower,dif); begin scalar z,w,uu; if null dif then << z := 'prod . (prepsq u . list car y); if w := opmtch z then return (nil . simp w) else if null cdr y then return (t . mksq(z,1)); load_package 'factor; % try to find factors uu := old_factorize prepf numr u; if length uu > 2 then << z := 'times . foreach uuu in cdr uu collect ('prod . ( prepsq multsq(if pairp uuu and eq(car uuu,'!*sq) then cadr uuu else simp uuu,1 ./ denr u)) . y); return (t . simp z) >>; z := 'prod . (prepsq u . y); % try to apply rules let prod_last_attempt_rules!*; w:= opmtch z; rule!-list (list prod_last_attempt_rules!*,nil); return (t . if w then simp w else mksq(z,1))>>; % return (t . if w := opmtch z then simp w else mksq(z,1))>>; z := 1 ./ 1; a: if dif < 0 then return (t . z); z := multsq(z,subsq(u,list(v . list('plus,lower,dif)))); dif := dif - 1; goto a end; prod_last_attempt_rules!* := algebraic << { prod(~f * ~g,~n,~anf,~ende) => prod(f,n,anf,ende) * prod(g,n,anf,ende) when g neq 1 and or(numberp prod(f,n,anf,ende), part(prod(f,n,anf,ende),0) neq prod, part(prod(g,n,anf,ende),0) neq prod), prod(~f / ~g,~n,~anf,~ende) => prod(f,n,anf,ende) / prod(g,n,anf,ende) when g neq 1 and or(numberp prod(f,n,anf,ende), % 1? part(prod(f,n,anf,ende),0) neq prod, part(prod(g,n,anf,ende),0) neq prod), prod(expt(~f,~k),~n,~anf,~ende) => (for ii:=1:k product prod(f,n,anf,ende)) when neq(part(prod(f,n,anf,ende),0),prod)} >>; %********************************************************************* % Product of rational function %********************************************************************; symbolic procedure prod!-sq(u,v); %ARGUMENT U : expression of s-q; % V : kernel; %value : expression of sq (result of product.); begin scalar gn,p1n,p2n,rn,qn,z,k,x,y; if null numr u then return 1 ./ 1; x := setkorder list v; qn := reorder numr u; rn := reorder denr u; if !*trsum then << prin2t " *** Product of A(n) = qn/rn with ***"; prin2!* "QN = ";xprinf(qn,nil,nil);terpri!* t; prin2!* "RN = ";xprinf(rn,nil,nil);terpri!* t>>; if nonpolyp(qn,v) or nonpolyp(rn,v) then go to fail; k := compress explode '!+j; z := integer!-root2(resultant(qn,nsubsf(rn,v,k),v),k); if !*trsum then << prin2 "Root of resultant(q(n),r(n+j)) are "; prin2t z >>; p2n := p1n := 1; while z do << k := car z; gn := gcdf!*(qn,nsubsf(rn,v,k)); qn := quotf!*(qn,gn); rn := quotf!*(rn,nsubsf(gn,v, -k)); if k > 0 then while (k := k - 1)>=0 do << p1n := multf(p1n,nsubsf(gn,v, -k)); if y := prod!-nsubsf(gn,v,-k) then p2n := multf(p2n,y)>> else if k < 0 then while k < 0 do << p2n := multf(p2n,nsubsf(gn,v, -k)); if y := prod!-nsubsf(gn,v,-k) then p1n := multf(p1n,y); k := k + 1>>; z := cdr z>>; if depend!-f(qn,v) or depend!-f(rn,v) then go to fail; u := multsq(p1n ./ p2n, simpexpt list(prepsq (qn ./ rn), v)); if !*trsum then << prin2t " *** Product of rational function calculated ***"; prin2!* " P(n) = "; sqprint u;terpri!* t; terpri!* t>>; setkorder x; return (reorder numr u ./ reorder denr u); return u; fail: if !*trsum then << prin2t " *** Product of rational function failed ***"; terpri!* t>>; setkorder x; return 'failed end; symbolic procedure prod!-nsubsf(u,kern,i); % ARGUMENT U : expression of sf; % KERN : kernel; % I : integer; % value : expression of sf; begin scalar x,z,n; x := setkorder list kern; u := reorder u; z := nil; a: if domainp u or not(mvar u eq kern) then goto b; z := addf(z,lc u); n := degr(u,kern) - degr(red u,kern); u := red u; if i = 0 then if n = 0 then nil else z := nil else z := multf(z,expt(i,n)); go to a; b: z := addf(z,u); setkorder x; return reorder z end; %********************************************************************* % integer (positive and negative) root isolation %********************************************************************; symbolic procedure integer!-root2(u,v); % Produce a list of all integer root of U; begin scalar x,root,n,w; x := setkorder list v; u := reorder u; if domainp u or not(mvar u eq v) then go to a; u := numr cancel(u ./ lc u); w := u; % get trailing term; while not domainp w and mvar w eq v and cdr w do w := cdr w; if (n := degr(w,v)) > 0 then << w := lc w; while n > 0 do << root := 0 . root; n := n - 1>>>>; n := dfactors lowcoef w; % factor tail coeff.; w := (v . 1) . 1; while n do << if not testdivide(u,v,car n) then << root := car n . root; u := quotf!*(u, (w . - car n))>> else if not testdivide(u,v,- car n) then << root := (- car n) . root; u := quotf!*(u, (w . car n))>> else n := cdr n>>; a: setkorder x; return root end; endmodule; end;