Artifact 243791274ddd1f632c1a5ea5c79e0cb04ae08dc8c950edd711b62c500536eb0c:
- Executable file
r38/packages/redlog/sfto.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: 21743) [annotate] [blame] [check-ins using] [more...]
% ---------------------------------------------------------------------- % $Id: sfto.red,v 1.14 2003/12/16 13:01:54 sturm Exp $ % ---------------------------------------------------------------------- % Copyright (c) 1995-1999 Andreas Dolzmann and Thomas Sturm % ---------------------------------------------------------------------- % $Log: sfto.red,v $ % Revision 1.14 2003/12/16 13:01:54 sturm % Fixed comment and removed unused scalar declaration in sfto_exteucd. % % Revision 1.13 2003/12/16 12:59:04 sturm % Added procedure sfto_exteucd for integer GCD's with cofactors. % AM frontend exteuc uses this now for integer arguments. % % Revision 1.12 2003/05/19 10:40:11 dolzmann % Changed comment. % % Revision 1.11 2003/01/13 10:00:52 dolzmann % Added predicated for linear terms and linear weak paramteric terms. % % Revision 1.10 2002/01/25 13:05:10 sturm % Added procedure sfto_dprpartksf(). % % Revision 1.9 2001/12/15 18:16:14 sturm % Added function sfto_exteucf() and AM interface. % % Revision 1.8 1999/03/22 15:26:43 dolzmann % Changed copyright information. % Added and reformatted comments. % % Revision 1.7 1999/01/17 15:33:13 dolzmann % Added procedure sfto_sqfpartz for computing the square-free part of an % integer. % % Revision 1.6 1999/01/10 12:09:16 dolzmann % Added procedures sfto_zdeqn, sfto_zdgtn, sfto_zdgen for the decomposition of % integers. % % Revision 1.5 1996/10/08 13:54:59 dolzmann % Renamed "degree parity decomposition" to "parity decomposition". % Adapted names of procedures and switches accordingly. % % Revision 1.4 1996/09/05 11:17:48 dolzmann % Added procedure sfto_monfp. % % Revision 1.3 1996/07/07 12:56:12 dolzmann % Fixed a bug in sfto_preducef and sfto_greducef. % % Revision 1.2 1996/05/13 13:54:26 dolzmann % Added procedure sfto_sqrtf. % % Revision 1.1 1996/04/30 12:06:46 sturm % Merged ioto, lto, and sfto into rltools. % % Revision 1.2 1996/04/30 09:13:33 sturm % Added procedure sfto_gcdf implementing the Davenport test. % % Revision 1.1 1996/03/22 12:19:17 sturm % Moved. % % Revision 1.5 1996/03/04 17:15:46 sturm % Added procedure sfto_decdegf. % Moved sfto_reorder from module ofsf to this module. % % Revision 1.4 1996/03/04 13:09:54 dolzmann % Moved procedures sfto_groebnerf, sfto_preducef, sfto_greducef and % loading of groebner packages from module ofsfgs to this module. % % Revision 1.3 1995/08/30 08:26:45 sturm % Fixed a bug in procedure sfto_dprpartf. % % Revision 1.2 1995/08/30 07:35:19 sturm % Added procedures sfto_dprpartf and sfto_dprpartf1. % % Revision 1.1 1995/05/29 14:47:23 sturm % Initial check-in. % % ---------------------------------------------------------------------- lisp << fluid '(sfto_rcsid!* sfto_copyright!*); sfto_rcsid!* := "$Id: sfto.red,v 1.14 2003/12/16 13:01:54 sturm Exp $"; sfto_copyright!* := "Copyright (c) 1995-1999 by A. Dolzmann and T. Sturm" >>; module sfto; % Standard form tools. load!-package 'groebner; load!-package 'groebnr2; fluid '(!*ezgcd !*gcd !*rldavgcd !*rational); switch sfto_yun,sfto_tobey,sfto_musser; !*sfto_yun := T; put('sqfpart,'polyfn,'sfto_sqfpartf); put('tsqsum,'psopfn,'sfto_tsqsum!$); put('sqfdec,'psopfn,'sfto_sqfdec!$); put('pdec,'psopfn,'sfto_pdec!$); put('sfto_yun,'simpfg, '((T (setq !*sfto_tobey nil) (setq !*sfto_musser nil)))); put('sfto_tobey,'simpfg, '((T (setq !*sfto_yun nil) (setq !*sfto_musser nil)))); put('sfto_musser,'simpfg, '((T (setq !*sfto_tobey nil) (setq !*sfto_yun nil)))); operator exteuc; procedure sfto_dcontentf(u); % Standard form tools domain content standard form. [u] is an SF. % Returns a domain element, which is the (non-negatice) content of % [u] as a multivariate polynomial over the current domain. sfto_dcontentf1(u,nil); procedure sfto_dcontentf1(u,g); % Standard form tools domain content standard form subroutine. [u] % is a term; [g] is a domain element. Returns the gcd of the % content of [u] and [g], which is a domain element. if g = 1 then g else if domainp u then sfto_gcdf(absf u,g) else sfto_dcontentf1(red u,sfto_dcontentf1(lc u,g)); procedure sfto_dprpartf(u); % Standard form tools domain primitive part standard form. [u] is % an SF. Returns an SF which is the primitive part of [u] as a % multivariate polynomial over the current domain. The sign of the % head coefficient in this primitive part is positive. sfto_dprpartf1(u,sfto_dcontentf u); procedure sfto_dprpartksf(u); % Standard form tools domain primitive part standard form keep % sign. [u] is an SF. Returns an SF which is the primitive part of % [u] as a multivariate polynomial over the current domain. The % sign of the head coefficient in this primitive part is that of % [u]. quotf(u,sfto_dcontentf u); procedure sfto_dprpartf1(u,c); % Standard form tools domain primitive part standard form % subroutine. [u] and [c] are SF's. Returns an SF which is the % primitive part of [u] as a multivariate polynomial over the % current domain. (if minusf w then negf w else w) where w = quotf(u,c); procedure sfto_sqfpartf(u); % Standard form tools square-free part. [u] is a non-zero SF. % Returns an SF which is the square-free part of [u] as a % multivariate polynomial. The (domain) content is dropped. begin scalar c,pp; if domainp u then return 1; c := sfto_ucontentf u; pp := quotf(u,c); return multf(sfto_sqfpartf(c),quotf(pp,sfto_gcdf!*(pp,diff(pp,mvar u)))) end; procedure sfto_ucontentf(u); % Standard form tools univariate content standard form. [u] is an % SF. Returns the content of [u] as a univariate polynomial in its % [mvar] over the polynomial ring in all other contained variables. if domainp u then u else sfto_ucontentf1(u,mvar u); procedure sfto_ucontentf1(u,v); % Standard form tools univariate content standard form subroutine. % [v] is a kernel; [u] is an SF with main variable [v]. Returns an % SF which is the content of [u] as an univariate polynomial in [v] % over the polynomial ring in all other contained variables. if domainp u or mvar u neq v then u else sfto_gcdf!*(lc u,sfto_ucontentf1(red u,v)); procedure sfto_uprpartf(u); % Standard form tools univariate primitive part. [u] is an SF. % Returns the primitive part of [u] as a univariate polynomial in % its [mvar] over the polynomial ring in all other contained % variables. quotf(u,sfto_ucontentf u); procedure sfto_tsqsumf(u); % Standard form tools trivial square sum standard form. [u] is an % SF. Returns one of [nil], ['stsq], or ['tsq]. ['stsq] means that % in the sparse distributive representation of [u] all exponents % are even and all coefficients are positive. ['tsq] means that all % exponents are even and all coefficients are positive except for % that there is no absolute summand. if domainp u then (if null u then 'tsq else if not minusf u then 'stsq) else evenp ldeg u and sfto_tsqsumf lc u and sfto_tsqsumf red u; procedure sfto_tsqsum!$(argl); sfto_tsqsumf(numr simp car argl); procedure sfto_sqfdecf(u); % Standard form tools multivariate square-free decomposition % standard form. [u] is an SF. Returns a (dense) list $((q_1 . % 1),(q_2 . 2),...,(q_n . n))$ such that $\prod q_i^i = u$ with the % $q_i$ square-free and pairwise relatively prime. The (integer) % content of u is dropped. Decomposition is performed by merging % univariate decompositions. The univariate decomposition method is % selected by turning on one of the switches [sfto_yun] (default), % [sfto_tobey], or [sfto_musser]. begin scalar c,pp; if domainp u then return {1 . 1}; c := sfto_ucontentf u; pp := quotf(u,c); return sfto_sqdmerge(sfto_sqfdecf(c),sfto_usqfdecf(pp)) end; procedure sfto_sqfdec!$(argl); % Standard form tools square free decomposition. [argl] is an % argument list. Returns an AM list of AM lists of the form % $(p_i,m_i)$, where the $p_i$ are polynomials represented as a % Lisp-prefix-form and the $m_i$ are integers. begin scalar w; return 'list . for each x in sfto_sqfdecf numr simp car argl join if (w := prepf car x) neq 1 then {{'list,w,cdr x}} end; procedure sfto_usqfdecf(u); if !*sfto_yun then sfto_yun!-usqfdecf u else if !*sfto_musser then sfto_musser!-usqfdecf u else if !*sfto_tobey then sfto_tobey!-usqfdecf u else rederr {"sfto_usqfdecf: select a decomposition method"}; procedure sfto_yun!-usqfdecf(p); % Standard form tools univariate square-free decomposition after % Yun. [p] is a an SF that is viewed a univariate Polynomial in its % [mvar] over the polynomial ring in all other variables; in this % sense, [p] must be primitive. Returns the square-free % decomposition of [p] as a (dense) list $((q_1 . 1),(q_2 . % 2),...,(q_n . n))$ such that $\prod q_i^i = u$ with the $q_i$ % square-free and pairwise relatively prime. begin scalar !*gcd,x,g,c,d,w,l; integer n; !*gcd := T; x := mvar p; g := sfto_gcdf(p,w := diff(p,x)); c := quotf(p,g); d := addf(quotf(w,g),negf(diff(c,x))); repeat << p := sfto_gcdf(c,d); l := (p . (n := n+1)) . l; c := quotf(c,p); d := addf(quotf(d,p),negf(diff(c,x))) >> until domainp c; return reversip l end; procedure sfto_musser!-usqfdecf(u); % Standard form tools univariate square-free decomposition after % Musser. [p] is a an SF that is viewed a univariate Polynomial in % its [mvar] over the polynomial ring in all other variables; in % this sense, [p] must be primitive. Returns the square-free % decomposition of [p] as a (dense) list $((q_1 . 1),(q_2 . % 2),...,(q_n . n))$ such that $\prod q_i^i = u$ with the $q_i$ % square-free and pairwise relatively prime. begin scalar !*gcd,v,u1,sqfp,sqfp1,l; integer n; !*gcd := T; v := mvar u; u1 := sfto_gcdf(u,diff(u,v)); sqfp := quotf(u,u1); while degr(u1,v)>0 do << sqfp1 := sfto_gcdf(sqfp,u1); l := (quotf(sqfp,sqfp1) . (n := n+1)) . l; u1 := quotf(u1,sqfp1); sqfp := sqfp1 >>; l := (sqfp . (n := n+1)) . l; return reversip l end; procedure sfto_tobey!-usqfdecf(u); % Standard form tools univariate square-free decomposition after % Tobey and Horowitz. [p] is a an SF that is viewed a univariate % Polynomial in its [mvar] over the polynomial ring in all other % variables; in this sense, [p] must be primitive. Returns the % square-free decomposition of [p] as a (dense) list $((q_1 . % 1),(q_2 . 2),...,(q_n . n))$ such that $\prod q_i^i = u$ with the % $q_i$ square-free and pairwise relatively prime. begin scalar !*gcd,v,h,q1,q2,l; integer n; !*gcd := T; v := mvar u; h := sfto_gcdf(u,diff(u,v)); q2 := quotf(u,h); while degr(u,v)>0 do << u := h; q1 := q2; h := sfto_gcdf(u,diff(u,v)); q2 := quotf(u,h); l := (quotf(q1,q2) . (n := n+1)) . l >>; return reversip l end; procedure sfto_sqdmerge(l1,l2); % Standard form tools square-free decomposition merge begin scalar l; l := l1; while l1 and l2 do << caar l1 := multf(caar l1,caar l2); l1 := cdr l1; l2 := cdr l2 >>; if l2 then l := nconc(l,l2); return l end; procedure sfto_pdecf(u); % Standard form tools multivariate parity decomposition. [u] is an % SF. Returns a consed pair $a . d$ such that $a$ is the product of % all square-free factors with an odd multiplicity in [u] and $d$ % is that of the even multiplicity square-free factors. The % (integer) content of u is dropped. Decomposition is performed by % merging univariate decompositions. The univariate decomposition % method is selected by turning on one of the switches [sfto_yun] % (default), [sfto_musser]. begin scalar c,dpdc,dpdpp; if domainp u then return 1 . 1; c := sfto_ucontentf u; dpdc := sfto_pdecf c; dpdpp := sfto_updecf quotf(u,c); return multf(car dpdc,car dpdpp) . multf(cdr dpdc,cdr dpdpp) end; procedure sfto_updecf(u); if !*sfto_yun then sfto_yun!-updecf u else if !*sfto_musser then sfto_musser!-updecf u else rederr {"sfto_updecf: select a decomposition method"}; procedure sfto_yun!-updecf(p); begin scalar !*gcd,x,g,c,d,w,l,od; !*gcd := T; l := 1 . 1; x := mvar p; g := sfto_gcdf(p,w := diff(p,x)); c := quotf(p,g); d := addf(quotf(w,g),negf(diff(c,x))); repeat << od := not od; p := sfto_gcdf(c,d); if od then car l := multf(car l,p) else cdr l := multf(cdr l,p); c := quotf(c,p); d := addf(quotf(d,p),negf(diff(c,x))) >> until domainp c; return l end; procedure sfto_musser!-updecf(u); begin scalar !*gcd,od,v,u1,sqfp,sqfp1,l; !*gcd := T; od := T; l := 1 . 1; v := mvar u; u1 := sfto_gcdf(u,diff(u,v)); sqfp := quotf(u,u1); while degr(u1,v)>0 do << sqfp1 := sfto_gcdf(sqfp,u1); if od then car l := multf(car l,quotf(sqfp,sqfp1)) else cdr l := multf(cdr l,quotf(sqfp,sqfp1)); u1 := quotf(u1,sqfp1); sqfp := sqfp1; od := not od >>; if od then car l := multf(car l,sqfp) else cdr l := multf(cdr l,sqfp); return l end; procedure sfto_pdec!$(argl); {'list,prepf car w,prepf cdr w} where w=sfto_pdecf numr simp car argl; procedure sfto_decdegf(u,k,n); % Standard form tools decrement degree standard form. [u] is an SF; % [k] is a variable; [n] is an integer. Returns an SF. Replace each % occurence of $[k]^d$ by $k^(d/n)$. reorder sfto_decdegf1(sfto_reorder(u,k),k,n); procedure sfto_decdegf1(u,k,n); % Standard form tools decrement degree standard form. [u] is an SF % with main variable [k]; [k] is a variable; [n] is an integer. % Returns an SF. Replace each occurence of $[k]^d$ by $k^(d/n)$. if degr(u,k)=0 then u else mvar u .** (ldeg u / n) .* lc u .+ sfto_decdegf1(red u,k,n); procedure sfto_reorder(u,v); % Standard form tools reorder. [u] is an SF; [v] is a kernel. % Returns the SF [u] reorderd wrt. [{v}]. begin scalar w; w := setkorder {v}; u := reorder u; setkorder w; return u end; procedure sfto_groebnerf(l); % Standard form tools Groebner calculation standard form. [l] is a % list of SF's. Returns a list of SF's. The returned list is the % reduced Groebner basis of [l] wrt. the current term order. begin scalar w; if null l then return nil; w := groebnereval {'list . for each sf in l collect prepf sf}; return for each x in cdr w collect numr simp x end; procedure sfto_preducef(f,gl); % Standard form tools polynomial reduction standard form. [f] is an % SF and [gl] a list of SF's. Returns the numerator of [f] reduced % wrt. [gl]. if null gl then f else if (null cdr gl) and (domainp car gl) then nil else numr simp preduceeval { prepf f,'list . for each sf in gl collect prepf sf}; procedure sfto_greducef(f,gl); % Standard form tools polynomial reduction standard form. [f] is an % SF and [gl] a list of SF's. Returns the SF [f] reduced wrt. a % Groebner basis of [gl]. if null gl then f else if (null cdr gl) and (domainp car gl) then nil else numr simp greduceeval { prepf f,'list . for each sf in gl collect prepf sf}; procedure sfto_gcdf!*(f,g); % Standard form tools greatest common divisor of standard forms. % [f] and [g] are SF's. Returns an SF, the GCD of [f] and [g]. % Compute the GCD of [f] and [g] via [gcdf!*] or [ezgcdf] according % to Davenport's criterion: If, in one polynomial, the number of % variables of a degree greater than 2 is greater than 1, then use % [ezgcd]. sfto_gcdf(f,g) where !*gcd=T; procedure sfto_gcdf(f,g); % Standard form tools greatest common divisor of standard forms. % [f] and [g] are SF's. Returns an SF, the GCD of [f] and [g]. % Compute the GCD of [f] and [g] via [gcdf!*] or [ezgcdf] according % to Davenport's criterion: If, in one polynomial, the number of % variables of a degree greater than 2 is greater than 1, then use % [ezgcd]. For computing the real gcd of [f] ang [g] this % procedures require, that [!*gcd] is set to [T]. if null !*rldavgcd then gcdf(f,g) else if sfto_davp(f,nil) or sfto_davp(g,nil) then gcdf(f,g) where !*ezgcd=nil else ezgcdf(f,g); procedure sfto_davp(f,badv); % Standard form tools Davenport predicate. [f] is an SF; [v] is a % kernel or [nil]. Returns Boolean. [T] means [gcdf] can be used. if domainp f then T else if ldeg f > 2 then if badv and mvar f neq badv then nil else sfto_davp(lc f,mvar f) and sfto_davp(red f,mvar f) else sfto_davp(lc f,badv) and sfto_davp(red f,badv); procedure sfto_sqrtf(f); % Standard form tools square root standard form. Returns [nil] or % an SF $g$, such that $g**2=[f]$. begin scalar a,c,w,sd,result; c := sfto_dcontentf(f); result := fix sqrt c; if result**2 neq c then return nil; sd := sfto_sqfdecf(f); w := sd; while sd do << a := car sd; sd := cdr sd; if not(evenp cdr a) and car a neq 1 then << sd := nil; a := 'break >> else result := multf(result,exptf(car a,cdr a / 2 )) >>; if a neq 'break and exptf(result,2) = f then return result end; procedure sfto_monfp(sf); % Standard form tools monomial predicate. [f] is an SF. Returns an % SF. Check if [sf] is of the form $a x_1 \dots x_n$ for a domain % element $a$ and kernels $x_i$. domainp sf or (null red sf and sfto_monfp lc sf); procedure sfto_sqfpartz(z); % Standard form tools square free part of an integer. [z] is an % integer with prime decomposition $p_1^{e_1}\cdots p_n^{e_n}$. % Returns $\prod \{p_i\}$. sfto_zdgen(z,0); procedure sfto_zdeqn(z,n); % Standard form tools z decomposition equal n. [z] is an integer % with prime decomposition $p_1^{e_1}\cdots p_n^{e_n}$; [n] is a % positive integer. Returns $\prod \{p_i:e_i=n\}$. for each x in zfactor z product if cdr x = n then car x else 1; procedure sfto_zdgtn(z,n); % Standard form tools z decomposition greater than n. [z] is an % integer with prime decomposition $p_1^{e_1}\cdots p_n^{e_n}$; [n] % is a positive integer. Returns $\prod \{p_i:e_i>n\}$. for each x in zfactor z product if cdr x > n then car x else 1; procedure sfto_zdgen(z,n); % Standard form tools z decomposition greater than or equal to n. % [z] is an integer with prime decomposition $p_1^{e_1}\cdots % p_n^{e_n}$; [n] is a positive integer. Returns $\prod % \{p_i:e_i\geq n\}$. for each x in zfactor z product if cdr x >= n then car x else 1; procedure sfto_exteucf(a,b); % Standard form tools extended Euclidean Algorithm for polynomials. % [a], [b] are SF's describing univariate polynomials both in the % same variable. Returns a list $(g,s,t)$ of SQ's such that % $\gcd([a],[b])=g=s[a]+t[b]$. If the GCD $g$ is a domain element, % then it is normalized to $1$. begin scalar w,s,tt,u,v,s1,t1,!*rational; on1 'rational; s := numr simp 1; tt := numr simp 0; u := numr simp 0; v := numr simp 1; while not null b do << w := qremf(a,b); a := b; b := cdr w; s1 := s; t1 := tt; s := u; tt := v; u := addf(s1,negf multf(car w,u)); v := addf(t1,negf multf(car w,v)) >>; if domainp a then << s := quotf(s,a); tt := quotf(tt,a); a := 1 >>; off1 'rational; return {resimp !*f2q a,resimp !*f2q s,resimp !*f2q tt} end; procedure exteuc(a,b); begin scalar af,bf,ka,kb; af := numr simp a; bf := numr simp b; if domainp af and domainp bf then return 'list . sfto_exteucd(a,b); ka := kernels af; if (ka and cdr ka) or (kb and cdr kb) then rederr "non-univariate input polynomial"; return 'list . for each x in sfto_exteucf(af,bf) collect prepsq x end; procedure sfto_exteucd(a,b); % Standard form tools extended Euclidean Algorithm for domain % elements (integers). [a], [b] are numbers. Returns a list % $(g,s,t)$ of numbers such that $g>=0$ and % $\gcd([a],[b])=g=s[a]+t[b]$. begin integer s,tt,u,v,s1,t1,q,r; s := 1; v := 1; while not eqn(b,0) do << q := a / b; r := remainder(a,b); a := b; b := r; s1 := s; t1 := tt; s := u; tt := v; u := s1 - q * u; v := t1 - q * v >>; if a < 0 then << s := -s; tt := -tt; a := -a >>; return {a,s,tt} end; procedure sfto_linp(f,vl); % Standard form tools linear predicate. [f] is a SF; [vl] is a list % of variables. Returns Bool. Returns [T] if [f] is linear in [vl]. sfto_linp1(f,vl,nil); procedure sfto_linp1(f,vl,oc); domainp f or (not(mvar f memq vl) and not(mvar f memq oc) and sfto_linp1(lc f,vl,oc) and sfto_linp1(red f,vl,oc)) or (mvar f memq vl and not(mvar f memq oc) and ldeg f = 1 and sfto_linp1(lc f,vl,mvar f . oc) and sfto_linp1(red f,vl,oc)); procedure sfto_linwpp(f,vl); % Standard form tools linear and weak parametric predicate. [f] is % a SF; [vl] is a list of variables. Returns Bool. Returns [T] if % [f] is linear and weak parametric in [vl]. domainp f or (not(mvar f memq vl) and null intersection(kernels lc f,vl) and sfto_linwpp(red f,vl)) or (mvar f memq vl and ldeg f = 1 and domainp lc f and sfto_linwpp(red f,vl)); endmodule; % [sfto] end; % of file