module facform; % Factored form representation for standard form polys.
% Author: Anthony C. Hearn.
% Modifications by: Francis J. Wright.
% Copyright (c) 1990 The RAND Corporation. All rights reserved.
% INTEGER FACTORS?
% SHOULDN'T SYMMETRIC TESTS ETC BE RUN RECURSIVELY?
fluid '(!*exp !*ezgcd !*factor !*force!-prime !*gcd !*ifactor !*nopowers
!*kernreverse !*limitedfactors !*sqfree !*trfac current!-modulus
dmode!* m!-image!-variable ncmp!*);
switch limitedfactors,nopowers;
% switch sqfree;
put('sqfree,'simpfg,'((t (rmsubs) (setq !*exp nil))
(nil (rmsubs) (setq !*exp t))));
comment In this module, we consider the manipulation of factored forms.
These have the structure
<monomial> . <form-power-list>
where the monomial is a standard form (with numerator and
denominator satisfying the KERNLP test) and a form-power is a
dotted pair whose car is a standard form and cdr an integer>0.
We have thus represented the form as a product of a monomial
quotient and powers of non-monomial factors;
symbolic procedure fac!-merge(u,v);
% Returns the merge of the factored forms U and V.
multf(car u,car v) . append(cdr u,cdr v);
symbolic procedure factorize u;
% Factorize the polynomial u, returning the factors found.
(begin scalar x,y;
x := simp!* u;
y := denr x;
if not domainp y then typerr(u,"polynomial");
u := numr x;
if u = 1 then return
{'list, if !*nopowers then 1 else {'list,1,1}} % FJW
else if fixp u then !*ifactor := t; % Factor an integer.
if !*force!-prime and not primep !*force!-prime
then typerr(!*force!-prime,"prime");
u := if dmode!* and not(dmode!* memq '(!:rd!: !:cr!:))
then if get(dmode!*,'factorfn)
then begin scalar !*factor;
!*factor := t;
return fctrf u
end
else rerror(poly,14,
list("Factorization not supported over domain",
get(dmode!*,'dname)))
else fctrf u;
return facform2list(u,y)
end) where !*ifactor = !*ifactor;
symbolic procedure facform2list(x,y);
% x is a factored form.
% y is a possible numerical (domain) denominator.
begin scalar factor!-count,z;
if null car x and null cdr x then return list 'list
% car x is now expected to be a number.
else if null !*nopowers then z := facform2list2 x
else <<
z:= (0 . car x) . nil;
x := reversip!* cdr x; % This puts factors in better order.
factor!-count:=0;
for each fff in x do
for i:=1:cdr fff do
z := ((factor!-count:=factor!-count+1) .
mk!*sq(car fff ./ 1)) . z;
z := multiple!-result(z,nil);
if atom z then typerr(z,"factor form") % old style input.
else if numberp cadr z and cadr z<0 and cddr z
then z := car z .
(- cadr z) . mk!*sq negsq simp caddr z
. cdddr z;
% make numerical coefficient positive.
z := cdr z;
if car z = 1 then z := cdr z
else if not fixp car z then z := prepd car z . cdr z
else if !*ifactor
then z := append(pairlist2list reversip zfactor car z,
cdr z)>>;
if y neq 1 then z := list('recip,prepd y) . z;
return 'list . z
end;
symbolic procedure facform2list2 u;
begin scalar bool,x;
if !:minusp(x := car u) then <<bool := t; x := !:minus x>>;
u := cdr u;
if x neq 1
then if !*ifactor and fixp x
then u := append(reversip zfactor x,u)
else u := (x . 1) . u;
% Adjust for negative sign.
x := nil;
for each j in u do
if bool and not evenp cdr j
then <<bool := nil; x := (negf car j . cdr j) . x>>
else x := j . x;
% Convert terms to list form.
u := nil;
for each j in x do
if fixp car j then u := {'list,car j,cdr j} . u
else u := {'list,mk!*sq(car j ./ 1),cdr j} . u;
return if bool then '(list -1 1) . u else u
end;
symbolic procedure old_factorize u; factorize u where !*nopowers=t;
flag('(factorize old_factorize),'opfn);
symbolic procedure pairlist2list u;
for each x in u conc nlist(car x,cdr x);
symbolic procedure fctrf u;
% U is a standard form. Value is a factored form.
% The function FACTORF is an assumed entry point to a more complete
% factorization module. It returns a form power list.
(begin scalar !*ezgcd,!*gcd,denom,x,y;
if domainp u then return list u
else if ncmp!* and not noncomfp u then ncmp!* := nil;
!*gcd := t;
if null !*limitedfactors and null dmode!* then !*ezgcd := t;
if null !*mcd
then rerror(poly,15,"Factorization invalid with MCD off")
else if null !*exp
then <<!*exp := t; u := !*q2f resimp !*f2q u>>;
% Convert rationals to integers for factorization.
if dmode!* eq '!:rn!:
then <<dmode!* := nil; alglist!* := nil . nil;
u := simp prepf u;
denom := denr u; u := numr u>>;
% Check for homogeneous polynomials. This can't be done with
% current code though if non-commuting objects occur.
if null ncmp!*
then <<x := sf2ss u;
if homogp x
then <<if !*trfac
then prin2t
"This polynomial is homogeneous - variables scaled";
y := caaar x . listsum caaadr x;
x := fctrf1 ss2sf(car(x)
. (reverse subs0 cadr x . 1));
x := rconst(y,x);
return car x . sort!-factors cdr x>>>>;
u := fctrf1 u;
if denom
then <<alglist!* := nil . nil;
dmode!* := '!:rn!:; car u := quotf(car u,denom)>>;
return car u . sort!-factors cdr u
end) where !*exp = !*exp, ncmp!* = ncmp!*;
symbolic procedure fctrf1 u;
begin scalar x,y,z;
if domainp u then return list u; % Do this again just in case.
if flagp(dmode!*,'field) and ((z := lnc u) neq 1)
then u := multd(!:recip z,u)
else if dmode!* and (y := get(dmode!*,'unitsfn))
then <<x := apply2(y,1 . u,lnc u);
u := cdr x;
z := !:recip car x>>;
x := comfac u;
u := quotf(u,comfac!-to!-poly x);
y := fctrf1 cdr x; % factor the content.
if car x then y := car y . (!*k2f caar x . cdar x) . cdr y;
if z and (z neq 1) then y := multd(z,car y) . cdr y;
if domainp u then return multf(u,car y) . cdr y
else if minusf u
then <<u := negf u; y := negf car y . cdr y>>;
return fac!-merge(factor!-prim!-f u,y)
end;
symbolic procedure factorize!-form!-recursion u;
fctrf1 u;
symbolic procedure factor!-prim!-f u;
% U is a non-trivial form which is primitive in all its variables
% and has a positive leading numerical coefficient. Result is a
% form power list.
begin scalar v,w,x,y;
if ncmp!* then return list(1,u . 1);
if dmode!* and (x := get(dmode!*,'sqfrfactorfn))
then if !*factor then v := apply1(x,u) else v := list(1,u . 1)
else if flagp(dmode!*,'field) and ((w := lnc u) neq 1)
then v := w . sqfrf multd(!:recip w,u)
else if (w := get(dmode!*,'units)) and (w := assoc(y := lnc u,w))
then v := y . sqfrf multd(cdr w,u)
else v := 1 . sqfrf u;
if x and (x eq get(dmode!*,'factorfn))
then return v; % No point in re-factorizing.
w := list car v;
for each x in cdr v
do w := fac!-merge(factor!-prim!-sqfree!-f x,w);
return w
end;
symbolic procedure factor!-prim!-sqfree!-f u;
% U is of the form <square free poly> . <power>. Result is a factored
% form.
% Modified to work properly in rounded (real or complex),
% rational and complex modes. SLK.
begin scalar x,y,!*msg,r;
r := !*rounded;
% It's probable that lc numr u and denr u will always be 1 if
% u is univariate.
if r and univariatep numr u and lc numr u=1 and denr u=1
then return unifactor u
else if r or !*complex or !*rational then
<<if r then on rational;
u := numr resimp !*f2q car u . cdr u>>;
if null !*limitedfactors
then <<if null dmode!* then y := 'factorf
else <<x := get(dmode!*,'sqfrfactorfn);
y := get(dmode!*,'factorfn);
if x and not(x eq y) then y := 'factorf>>;
if y
then <<y := apply1(y,car u);
u := (exptf(car y,cdr u) . for each j in cdr y
collect(car j . cdr u));
go to ret>>>>;
u := factor!-prim!-sqfree!-f!-1(car u,cdr u);
ret: if r then
<<on rounded;
u := car u . for each j in cdr u collect
(numr resimp !*f2q car j . cdr j)>>;
return u
end;
symbolic procedure unifactor u;
if not eqcar(u := root_val list mk!*sq u,'list)
then errach {"unifactor1",u}
else 1 . for each j in cdr u collect
if not eqcar(j,'equal) then errach{"unifactor2",u}
else addsq(!*k2q cadr j,negsq simp caddr j);
symbolic procedure distribute!.multiplicity(factorlist,n);
% Factorlist is a simple list of factors of a square free primitive
% multivariate poly and n is their multiplicity in a square free
% decomposition of another polynomial. result is a list of form:
% ((f1 . n),(f2 . n),...) where fi are the factors.
for each w in factorlist collect (w . n);
symbolic procedure factorf u;
% NOTE: This is not an entry point to be used by novice programmers.
% Please used FCTRF instead.
% There is an inefficiency in this procedure relating to ordering.
% There is a final reordering done at every recursive level in order
% to make sure final result has the right order. However, this only
% need be done once at the top level, probably in fctrf. Since some
% programmers still use this function however, it's safer for it to
% return results in the correct order.
(begin scalar m!-image!-variable,new!-korder,old!-korder,sign,v,w;
if domainp u then return list u;
new!-korder:=kernord(u,nil);
if !*kernreverse then new!-korder:=reverse new!-korder;
old!-korder:=setkorder new!-korder;
u := reorder u; % Make var of lowest degree the main one.
if minusf u then <<sign := not sign; u := negf u>>;
v := comfac u; % Since new order may reveal more factors.
u := quotf1(u,cdr v);
if domainp u then errach list("Improper factors in factorf");
% The example on rounded; solve(df(e^x/(e^(2*x)+1)^1.5,x),{x});
% shows car v can be non-zero.
w := car v;
v := fctrf1 cdr v;
if w then v := car v . (!*k2f car w . cdr w) . cdr v;
m!-image!-variable := mvar u;
u :=
distribute!.multiplicity(factorize!-primitive!-polynomial u,1);
setkorder old!-korder;
if sign then u := (negf caar u . cdar u) . cdr u;
u := fac!-merge(v,1 . u);
return car u . for each w in cdr u collect (reorder car w . cdr w)
end) where current!-modulus = current!-modulus;
symbolic procedure factor!-prim!-sqfree!-f!-1(u,n);
(exptf(car x,n) . for each j in cdr x collect (j . n))
where x = prsqfrfacf u;
symbolic procedure sqfrf u;
% U is a non-trivial form which is primitive in all its variables
% and has an overall numerical coefficient which should be a unit.
% SQFRF performs square free factorization on U and returns a
% form power list.
% Modified to work properly in rounded (real or complex) modes. SLK.
begin integer n; scalar !*gcd,units,v,w,x,y,z,!*msg,r;
!*gcd := t;
if (r := !*rounded) then
<<on rational; u := numr resimp !*f2q u>>;
n := 1;
x := mvar u;
% With ezgcd off, some sqrts can take a long, long time.
v := gcdf(u,diff(u,x)) where !*ezgcd = t;
u := quotf(u,v);
% If domain is a field, or has non-trivial units, v can have a
% spurious numerical factor.
if flagp(dmode!*,'field) and ((y := lnc u) neq 1)
then <<u := multd(!:recip y,u); v := multd(y,v)>>
% The following check for units can result in the loss of such
% a unit.
% else if (units := get(dmode!*,'units))
% and (w := assoc(y:= lnc u,units))
% then <<u := multd(cdr w,u); v := multd(y,v)>>;
;
while degr(v,x)>0 do
<<w := gcdf(v,u);
if u neq w then z := (quotf(u,w) . n) . z;
% Don't add units to z.
v := quotf(v,w);
u := w;
n := n + 1>>;
if r then
<<on rounded;
u := numr resimp !*f2q u;
z := for each j in z
collect numr resimp !*f2q car j . cdr j>>;
if v neq 1 and assoc(v,units) then v := 1;
if v neq 1 then if n=1 then u := multf(v,u)
else if (w := rassoc(1,z)) then rplaca(w,multf(v,car w))
else if null z and ((w := rootxf(v,n)) neq 'failed)
then u := multf(w,u)
else if not domainp v then z := aconc(z,v . 1)
else errach {"sqfrf failure",u,n,z};
return (u . n) . z
end;
symbolic procedure square_free u;
'list . for each v in sqfrf !*q2f simp!* u
collect {'list,mk!*sq(car v . 1),cdr v};
flag('(square_free),'opfn);
symbolic procedure diff(u,v);
% A polynomial differentation routine which does not check
% indeterminate dependencies.
if domainp u then nil
else addf(addf(multpf(lpow u,diff(lc u,v)),
multf(lc u,diffp1(lpow u,v))),
diff(red u,v));
symbolic procedure diffp1(u,v);
if not( car u eq v) then nil
else if cdr u=1 then 1
else multd(cdr u,!*p2f(car u .** (cdr u-1)));
endmodule;
end;