module gcd; % Greatest common divisor routines.
% Author: Anthony C. Hearn.
% Copyright (c) 1987 The RAND Corporation. All rights reserved.
fluid '(!*exp !*ezgcd !*gcd !*heugcd !*mcd asymplis!* dmode!*);
switch ezgcd,heugcd;
% Note: The handling of non-commuting quantities in the following is
% dubious. The problem is that to do things properly, a left- and
% right-hand gcd and quotient would be necessary. For now, the code
% returns 1 if the quotient tests fail in gcdf1 for non-commuting
% arguments.
symbolic procedure comfac p;
% P is a non-atomic standard form.
% CAR of result is lowest common power of leading kernel in
% every term in P (or NIL). CDR is gcd of all coefficients of
% powers of leading kernel.
% If field elements are involved, lnc is normalized to 1.
% We need GCDF here since the same function is used by EZGCD.
begin scalar x,y;
if flagp(dmode!*,'field) and ((x := lnc p) neq 1)
then p := multd(!:recip x,p);
if null red p then return lt p;
x := lc p;
y := mvar p;
a: p := red p;
if degr(p,y)=0
then return nil
. if domainp p or not(noncomp y and noncomp mvar p)
then gcdf(x,p)
else 1
else if null red p then return lpow p . gcdf(x,lc p)
else x := gcdf(lc p,x);
go to a
end;
symbolic procedure degr(u,var);
if domainp u or not(mvar u eq var) then 0 else ldeg u;
put('gcd,'polyfn,'gcdf!*);
put('gcd,'number!-of!-args,2);
symbolic procedure gcdf!*(u,v);
begin scalar !*gcd; !*gcd := t; return gcdf(u,v) end;
symbolic procedure gcdf(u,v);
% U and V are standard forms.
% Value is the gcd of U and V, complete only if *GCD is true.
begin scalar !*exp,!*rounded;
% The next line was to prevent numerators moving to denominators
% as in weight x=1,y=2$ wtlevel 4$ wtest:=(z^4-z^3*y-z^3*x+z^2*y^2
% +2*z^2*x*y+z^2*x^2-3*z*x^2*y-z*x^3+x^4)/z^5; wtest where z=>a;
% However, the results are formally correct without it, and it
% causes other problems.
% if wtl!* then return 1;
!*exp := t;
u := if domainp u or domainp v or not !*ezgcd
% or dmode!* memq '(!:rn!: !:rd!:) % Should be generalized.
or dmode!* % I don't know what to do in this case.
or free!-powerp u or free!-powerp v
then gcdf1(u,v)
else ezgcdf(u,v);
return if minusf u then negf u else u
end;
symbolic procedure free!-powerp u;
not domainp u
and (not fixp ldeg u or free!-powerp lc u or free!-powerp red u);
symbolic procedure gcdf1(u,v);
begin scalar w;
if null u then return v
else if null v then return u
else if u=1 or v=1 then return 1
else if domainp u then return gcdfd(u,v)
else if domainp v then return gcdfd(v,u)
else if not num!-exponents u or not num!-exponents v then 1
else if quotf1(u,v) then return v
else if quotf1(v,u) then return u;
w := gcdf2(u,v);
if !*gcd and not(dmode!* memq '(!:rd!: !:cr!:))
and (null quotf1(u,w) or null quotf1(v,w))
then if noncomfp u or noncomfp v then return 1
else errach list("gcdf failed",prepf u,prepf v);
% This probably implies that integer overflow occurred.
return w
end;
symbolic procedure gcdf2(u,v);
% U and V are both non-trivial forms. Value is their GCD.
% We need to rebind asymplis!* to avoid setting higher powers to 0.
begin scalar asymplis!*,w,z;
if not num!-exponents u or not num!-exponents v then return 1;
if !*gcd and length(w := kernord(u,v))>1
then <<w := list setkorder w; % List used to make sure non-NIL
u := reorder u;
v := reorder v>>
else w := nil;
% Things can go wrong with noncom oprs. However, empirically we
% only need to make sure that both u and v do not have a leading
% noncom opr.
if mvar u eq mvar v
then begin scalar x,y;
x := comfac u;
y := comfac v;
z := gcdf1(cdr x,cdr y);
u := quotf1(u,comfac!-to!-poly x);
v := quotf1(v,comfac!-to!-poly y);
if !*gcd then z := multf(gcdk(u,v),z)
else if v and quotf1(u,v) then z := multf(v,z)
else if u and quotf1(v,u) then z := multf(u,z);
if car x and car y
then if pdeg car x>pdeg car y
then z := multpf(car y,z)
else z := multpf(car x,z)
end
else if noncomp mvar u and noncomp mvar v
then z := gcdfnc(u,v,mvar v)
else if ordop(mvar u,mvar v) then z := gcdf1(cdr comfac u,v)
else z := gcdf1(cdr comfac v,u);
if w then <<setkorder car w; z := reorder z>>;
return z
end;
symbolic procedure gcdfnc(x,p,y);
if domainp x or not noncomp mvar x then gcdf1(x,p)
else if null red x then gcdfnc(lc x,p,y)
else gcdf1(gcdfnc(lc x,p,y),gcdfnc(red x,p,y));
symbolic procedure num!-exponents u;
% check that all exponents are integers (this may not be true in
% rules).
domainp u or
fixp ldeg u and num!-exponents lc u and num!-exponents red u;
symbolic procedure gcdfd(u,v);
% U is a domain element, V a form.
% Value is gcd of U and V.
% if not atom u and flagp(car u,'field) then 1 else gcdfd1(u,v);
if flagp(dmode!*,'field) then 1 else gcdfd1(u,v);
symbolic procedure gcdfd1(u,v);
if null v then u
else if domainp v then gcddd(u,v)
else gcdfd1(gcdfd1(u,lc v),red v);
symbolic procedure gcddd(u,v);
%U and V are domain elements. If they are invertable, value is 1
%otherwise the gcd of U and V as a domain element;
if u=1 or v=1 then 1
% else if atom u and atom v then gcdn(u,v)
else if atom u then if atom v then gcdn(u,v)
else if fieldp v then 1
else dcombine(u,v,'gcd)
else if atom v
then if flagp(car u,'field) then 1 else dcombine(u,v,'gcd)
else if flagp(car u,'field) or flagp(car v,'field) then 1
else dcombine(u,v,'gcd);
symbolic procedure gcdk(u,v);
% U and V are primitive polynomials in the main variable VAR.
% Result is gcd of U and V.
begin scalar lclst,var,w,x;
if u=v then return u
else if domainp u or degr(v,(var := mvar u))=0 then return 1
else if ldeg u<ldeg v then <<w := u; u := v; v := w>>;
if quotf1(u,v) then return v
else if !*heugcd and (x := heu!-gcd(u,v)) then return x
% else if flagp(dmode!*,'field) then return 1
% otherwise problems arise.
else if ldeg v=1
or getd 'modular!-multicheck and modular!-multicheck(u,v,var)
or not !*mcd
then return 1;
a: w := remk(u,v);
if null w then return v
else if degr(w,var)=0 then return 1;
lclst := addlc(v,lclst);
if x := quotf1(w,lc w) then w := x
else for each y in lclst do
if atom y or not % prevent endless loop in !:gi!: dmode.
(domainp y and (x := get(car y,'units))
and y member (for each z in x collect car z))
then while (x := quotf1(w,y)) do w := x;
u := v; v := prim!-part w;
if degr(v,var)=0 then return 1 else go to a
end;
symbolic procedure addlc(u,v);
if u=1 then v
else (lambda x;
if x=1 or x=-1 or not atom x and flagp(car x,'field) then v
else x . v)
lc u;
symbolic procedure delallasc(u,v);
if null v then nil
else if u eq caar v then delallasc(u,cdr v)
else car v . delallasc(u,cdr v);
symbolic procedure kernord(u,v);
<<u := kernord!-split(u,v);
append(kernord!-sort car u,kernord!-sort cdr u)>>;
symbolic procedure kernord!-split(u,v);
% splits U and V into a set of powers of those kernels occurring in
% one form and not the other, and those occurring in both;
begin scalar x,y;
u := powers u;
v := powers v;
for each j in u do
if assoc(car j,v) then y := j . y else x := j . x;
for each j in v do
if assoc(car j,u) then y := j . y else x := j . x;
return reversip x . reversip y
end;
symbolic procedure kernord!-sort u;
% returns list of kernels ordered so that kernel with lowest maximum
% power in U (a list of powers) is first, and so on;
begin scalar x,y;
while u do
<<x := maxdeg(cdr u,car u);
u := delallasc(car x,u);
y := car x . y>>;
return y
end;
symbolic procedure maxdeg(u,v);
if null u then v
else if cdar u>cdr v then maxdeg(cdr u,car u)
else maxdeg(cdr u,v);
symbolic procedure powers form;
% returns a list of the maximum powers of each kernel in FORM.
% order tends to be opposite to original order.
powers0(form,nil);
symbolic procedure powers0(form,powlst);
if null form or domainp form then powlst
else begin scalar x;
if (x := atsoc(mvar form,powlst))
% then ldeg form>cdr x and rplacd(x,ldeg form)
then (if ldeg form>cdr x
then powlst := repasc(mvar form,ldeg form,powlst))
else powlst := (mvar form . ldeg form) . powlst;
return powers0(red form,powers0(lc form,powlst))
end;
put('lcm,'polyfn,'lcm!*);
put('lcm,'number!-of!-args,2);
symbolic procedure lcm!*(u,v);
begin scalar !*gcd; !*gcd := t; return lcm(u,v) end;
symbolic procedure lcm(u,v);
%U and V are standard forms. Value is lcm of U and V;
if null u or null v then nil
else if u=1 then v % ONEP
else if v=1 then u % ONEP
else multf(u,quotf(v,gcdf(u,v)));
symbolic procedure remk(u,v);
%modified pseudo-remainder algorithm
%U and V are polynomials, value is modified prem of U and V;
begin scalar f1,var,x; integer k,n;
f1 := lc v;
var := mvar v;
n := ldeg v;
while (k := degr(u,var)-n)>=0 do
<<x := negf multf(lc u,red v);
if k>0 then x := multpf(var .** k,x);
u := addf(multf(f1,red u),x)>>;
return u
end;
symbolic procedure prim!-part u;
%returns the primitive part of the polynomial U wrt leading var;
quotf1(u,comfac!-to!-poly comfac u);
symbolic procedure comfac!-to!-poly u;
if null car u then cdr u else list u;
endmodule;
end;