module heugcd;
%Authors: James Davenport & Julian Padget
% For full details of the algorithms see first Char et al. from the
% Proceedings of EUROSAM 84 (Springer LNCS #174), then Davenport and
% Padget in the Proceedings of EUROCAL 85 (Springer LNCS #204) and
% Davenport and Padget in the proceedings of Calcul Formel (Homage a
% Noel Gastinel) published by Masson-Wiley (France).
%exports heu!-gcd, heu!-gcd!-list;
%imports to be determined
%internal-functions
% univariatep, univariatep1, htc, kontent, kontent1,
% horner!-eval!-rat, horner!-eval!-rat!-and!-gcdl,
% horner!-eval!-rat!-and!-gcdl1, heu!-quotfl, heu!-quotfl1,
% heu!-quotf, xceiling, next!-even!-value, next!-odd!-value,
% heu!-gcdl, analyse!-polynomials, negshiftz, gen!-poly,
% gen!-poly!-forward, gen!-poly!-backward, gcdlist2, gcdf2
fluid '(!*heugcd reduction!-count);
global '(ee);
% ****************** Various polynomial utilities **********************
symbolic smacro procedure univariatep p; univariatep1(p,mvar p);
symbolic procedure univariatep1(p,v);
% checks that p is univariate in v;
if atom p then t
else if mvar p neq v then nil
else if atom lc p then univariatep1(red p,v)
else nil;
symbolic procedure htc p;
if atom p then p
else if null red p then lc p
else htc red p;
symbolic procedure kontent p;
% extract integer content of polynomial p
if domainp p then
if numberp p then p
else if null p then 1
else rederr "HEUGCD(kontent): unsupported domain element"
else if domainp red p then
if numberp red p then gcdn(lc p,red p)
else if null red p then lc p
else rederr "HEUGCD(kontent): unsupported domain element"
else kontent1(red red p,gcdn(lc p,lc red p));
symbolic procedure kontent1(p,a);
if a=1 then 1
else if domainp p then
if numberp p then gcdn(p,a)
else if null p then a
else rederr "HEUGCD(kontent1): unsupported domain element"
else kontent1(red p,gcdn(remainder(lc p,a),a));
symbolic procedure horner!-eval!-rat(p,v);
% evaluate the (sparse univariate) polynomial p at a rational v using
% Horner's scheme. Denominators are cleared by in fact calculating the
% following:
%
% for i:=min:max sum (a[i] * n**(i-min) * d**(max-min-i))
%
% note that if the polynomial does not end in a non-zero constant
% the routine it return the evaluation of p/(trailing exponent)
% s accumulates d**(max-min-i)
% ans accumulates the sum
% m is degree difference between current and previous term
% See specific routines below for further detail
if (numr v)=1 then horner!-eval!-integer(p,denr v,1,0)
else if (denr v)=1 then horner!-eval!-reciprocal(p,numr v,0,0)
else horner!-eval!-rational(p,numr v,denr v,0,1,0);
symbolic procedure horner!-eval!-rational(p,n,d,m,s,ans);
% general case of an arbitrary rational
if domainp p then
if p then ans*n**m+s*p else ans
else (lambda mp;
horner!-eval!-rational(red p,n,d,mp,s*d**mp,ans*n**m+s*lc p))
(ldeg p)-(if domainp red p then 0 else ldeg red p);
symbolic procedure horner!-eval!-integer(p,d,s,ans);
% simple sub case of an integer (n/1)
if domainp p then
if p then ans+s*p else ans
else horner!-eval!-integer(red p,d,
s*d**((ldeg p)-(if domainp red p then 0 else ldeg red p)),
ans+s*lc p);
symbolic procedure horner!-eval!-reciprocal(p,n,m,ans);
% simpler sub case of a straight reciprocal of an integer (1/n)
if domainp p then
if p then ans*n**m+p else ans
else horner!-eval!-reciprocal(red p,n,
(ldeg p)-(if domainp red p then 0 else ldeg red p),
ans*n**m+lc p);
symbolic procedure horner!-eval!-rat!-and!-gcdl(l,v);
% l is a list of polynomials to be evaluated at the point v
% and then take the GCD of these evaluations. We use an auxiliary
% routine with an accumulator variable to make the computation
% tail-recursive
if null cdr l then horner!-eval!-rat(car l,v)
else if null cddr l then
gcdn(horner!-eval!-rat(car l,v),horner!-eval!-rat(cadr l,v))
else horner!-eval!-rat!-and!-gcdl1(cddr l,v,
gcdn(horner!-eval!-rat(car l,v),horner!-eval!-rat(cadr l,v)));
symbolic procedure horner!-eval!-rat!-and!-gcdl1(l,v,a);
if a=1 then 1
else if null l then a
else horner!-eval!-rat!-and!-gcdl1(cdr l,v,
gcdn(horner!-eval!-rat(car l,v),a));
%*********** Polynomial division utilities and extensions *************
symbolic procedure heu!-quotfl(l,d);
% test division of each of a list of SF's (l) by the SF d
if null cdr l then heu!-quotf(car l,d)
else heu!-quotfl1(cdr l,d,heu!-quotf(car l,d));
symbolic procedure heu!-quotfl1(l,d,flag);
if null flag then nil
else if null cdr l then heu!-quotf(car l,d)
else heu!-quotfl1(cdr l,d,heu!-quotf(car l,d));
symbolic procedure heu!-quotf(p,q);
if domainp q then
if domainp p then
if null p then nil
else if null q then
rederr "HEUGCD(heu-quotf): division by zero"
else (lambda temp; if cdr temp=0 then car temp else nil)
divide(p,q)
else quotf(p,q)
else if domainp p then nil
else if ldeg p<ldeg q then nil
else if (cdr divide(lc p,lc q)) neq 0 then nil
else if p=q then 1
else (lambda qv;
if qv=0 then quotf(p,q)
else if remainder(horner!-eval!-rat(p,'(2 . 1)),qv)=0 then
quotf(p,q)
else nil)
horner!-eval!-rat(q,'(2 . 1));
%****************** Z-adic polynomial GCD routines ********************
symbolic smacro procedure xceiling(n,d); (n+d-1)/d;
symbolic smacro procedure force!-even x;
if evenp x then x else x+1;
symbolic smacro procedure force!-odd x;
if evenp x then x+1 else x;
symbolic smacro procedure next!-even!-value x;
if (denr x)=1 then force!-even fix(numr x * ee) ./ 1
else 1 ./ force!-even fix(denr x * ee);
symbolic smacro procedure next!-odd!-value x;
if (denr x)=1 then force!-odd fix(numr x * ee) ./ 1
else 1 ./ force!-odd fix(denr x * ee);
symbolic smacro procedure first!-value(inp,inq,lcp,lcq,tcp,tcq);
% Initial evaluation is based on Cauchy's inequality.
if lcp<tcp then
if lcq<tcq then
if (inp*tcq)<(inq*tcp) then 1 . (2+2*xceiling(inp,tcp))
else 1 . (2+2*xceiling(inq,tcq))
else if (inp*lcq)<(inq*tcp) then 1 . (2+2*xceiling(inp,tcp))
else (2+2*xceiling(inq,lcq)) . 1
else if lcq<tcq then
if (inp*tcq)<(inq*lcp) then (2+2*xceiling(inp,lcp)) . 1
else 1 . (2+2*xceiling(inq,tcq))
else if (inp*lcq)<(inq*lcp) then (2+2*xceiling(inp,lcp)) . 1
else (2+2*xceiling(inq,lcq)) . 1;
symbolic smacro procedure
second!-value(inp,inq,lcp,lcq,lgcd,tcp,tcq,tgcd);
% The second evaluation point is from a modified Mignotte bound.
(lambda (inp,inq,lcp,lcq,tcp,tcq);
if lcp<tcp then
if lcq<tcq then
if (inp*tcq)<(inq*tcp) then
1 . force!-even (2+xceiling(inp,tcp))
else
1 . force!-even (2+xceiling(inq,tcq))
else if (inp*lcq)<(inq*tcp) then
1 . force!-even (2+xceiling(inp,tcp))
else force!-even (2+xceiling(inq,lcq)) . 1
else if lcq<tcq then
if (inp*tcq)<(inq*lcp) then
force!-even (2+xceiling(inp,lcp)) . 1
else 1 . force!-even (2+xceiling(inq,tcq))
else if (inp*lcq)<(inq*lcp) then
force!-even (2+xceiling(inp,lcp)) . 1
else force!-even (2+xceiling(inq,lcq)) . 1)
(inp,inq,max(2,lcp/lgcd),max(2,lcq/lgcd),
max(2,tcp/tgcd),max(2,tcq/tgcd));
symbolic procedure heu!-gcd!-list l;
if null cdr l then car l
else if null cddr l then heu!-gcd(car l,cadr l)
else heu!-gcdl sort(l,function (lambda (p1,p2); ldeg p1<ldeg p2));
symbolic procedure heu!-gcdl l;
% Heuristic univariate polynomial GCD after Davenport & Padgets'
% extensions of Geddes' algorithm (EUROSAM 84) for a list of polynomials
begin scalar k,value,dval,d,xsx,inp,inq,lcp,lcq,lgcd,tcp,tcq,tgcd,tmp;
% check if first one is linear (input is sorted by leading degree)
if (ldeg car l=1) then return
(lambda pcarl; if heu!-quotfl(cdr l,pcarl) then pcarl else 1)
quotf(car l,kontent car l);
% general case - compute GCD of all of them at Cauchy's bound
tmp:=analyse!-polynomial car l;
if tmp then <<
inp:=car tmp; lcp:=lc car l; xsx:=cadr tmp; tcp:=caddr tmp;
tmp:=analyse!-polynomial cadr l;
if tmp then <<
inq:=car tmp; lcq:=lc cadr l;
xsx:=min(xsx,cadr tmp); tcq:=caddr tmp>>
else return nil>>
else return nil;
value:=first!-value(inp,inq,lcp,lcq,tcp,tcq);
% first check for trivial GCD
d:=gen!-poly(horner!-eval!-rat!-and!-gcdl(l,value),
value,mvar car l,xsx);
if heu!-quotfl(l,d) then return d;
% since that failed we pick a much higher evaluation point
% courtesy of a modified Mignotte inequality and just work on the
% first two
lgcd:=gcdn(lcp,lcq);
for each x in cddr l do lgcd:=gcdn(lc x,lgcd);
tgcd:=gcdn(tcp,tcq);
for each x in cddr l do tgcd:=gcdn(htc x,tgcd);
value:=second!-value(inp,inq,lcp,lcq,lgcd,tcp,tcq,tgcd);
loop:
d:=gen!-poly(horner!-eval!-rat!-and!-gcdl(l,value),
value,mvar car l,xsx);
if heu!-quotfl(l,d) then return d;
value:=next!-odd!-value value;
k:=k+1;
d:=gen!-poly(horner!-eval!-rat!-and!-gcdl(l,value),
value,mvar car l,xsx);
if heu!-quotfl(l,d) then return d;
value:=next!-even!-value value;
k:=k+1;
if k < 10 then goto loop;
print "(HEUGCD):heu-gcd-list fails";
return nil
end;
symbolic procedure heu!-gcd(p,q);
% Heuristic univariate polynomial GCD after Davenport & Padgets'
% extensions of Geddes' algorithm (EUROSAM 84)
% the method of choosing the evaluation point is quite complex (but not
% as general as it ought to be). It is
%
% min(infinity!-norm p/lc p, infinity!-norm p/htc p,
% infinity!-norm q/lc q, infinity!-norm q/htc q)
%
begin scalar k,value,d,dval,xsx,inp,inq,lcp,lcq,lgcd,tcp,tcq,tgcd,tmp;
% check if one of p and q is linear
if (ldeg q=1) or (ldeg p=1) then return
if univariatep p and univariatep q then
(lambda (pp,pq);
if (ldeg pq)=1 then
(lambda h; if null h then 1 else pq) heu!-quotf(pp,pq)
else
(lambda h; if null h then 1 else pp) heu!-quotf(pq,pp))
(quotf(p, kontent p), quotf(q, kontent q))
else nil;
% general case
if (ldeg p)>(ldeg q) then return heu!-gcd(q,p);
tmp:=analyse!-polynomial p;
if tmp then <<
inp:=car tmp; lcp:=lc p; xsx:=cadr tmp; tcp:=caddr tmp;
tmp:=analyse!-polynomial q;
if tmp then <<
inq:=car tmp; lcq:=lc q; xsx:=min(xsx,cadr tmp); tcq:=caddr tmp>>
else return nil>>
else return nil;
value:=first!-value(inp,inq,lcp,lcq,tcp,tcq);
% first check for trivial GCD
dval:=gcdn(horner!-eval!-rat(p,value),horner!-eval!-rat(q,value));
d:=gen!-poly(dval,value,mvar p,xsx);
if heu!-quotf(p,d) and heu!-quotf(q,d) then return d;
% if that failed we pick a much higher evaluation point
lgcd:=gcdn(lcp,lcq);
tgcd:=gcdn(lcp,lcq);
value:=second!-value(inp,inq,lcp,lcq,lgcd,tcp,tcq,tgcd);
k:=0;
loop:
dval:=gcdn(horner!-eval!-rat(p,value),horner!-eval!-rat(q,value));
d:=gen!-poly(dval,value,mvar p,xsx);
if heu!-quotf(p,d) and heu!-quotf(q,d) then return d;
value:=next!-odd!-value value;
k:=k+1;
dval:=gcdn(horner!-eval!-rat(p,value),horner!-eval!-rat(q,value));
d:=gen!-poly(dval,value,mvar p,xsx);
if heu!-quotf(p,d) and heu!-quotf(q,d) then return d;
value:=next!-even!-value value;
k:=k+1;
if k < 10 then goto loop;
print "(HEUGCD):heu-gcd fails";
return nil
end;
symbolic procedure analyse!-polynomial p;
% Determine the infinity norm of p and take note of the trailing
% coefficient, simultaneously check that p is univariate and take note
% of any trailing powers of the main variable. The result is a triple
% of (infinity-norm,excess powers,trailing coefficient)
analyse!-polynomial1(p,1,lc p,0,mvar p);
symbolic procedure analyse!-polynomial1 (p,inp,tcp,xsxp,mvarp);
if domainp p then
if p then list(max(inp,abs p),0,abs p)
else list(inp,xsxp,abs tcp)
else if ((mvar p) neq mvarp) then nil
else if domainp lc p then
analyse!-polynomial1(red p,max(inp,abs lc p),lc p,ldeg p,mvarp)
else nil;
%********** Reconstruction from the Z-adic representation *************
% given a number in [0,modulus), return the equivalent
% member of [-modulus/2,modulus/2)
% LAMBDA to ensure only one evaluation of arguments;
symbolic smacro procedure negshiftz(n,modulus);
(lambda (nn,mmodulus);
if nn>quotient(mmodulus,2) then nn-mmodulus else nn)
(n,modulus);
symbolic procedure gen!-poly(dval,value,var,xsx);
if (numr value)=1 then gen!-poly!-backward(dval,denr value,var,xsx)
else if (denr value)=1 then gen!-poly!-forward(dval,numr value,var,xsx)
else rederr "HEUGCD(gen-poly):point must be integral or reciprocal";
symbolic procedure gen!-poly!-forward(dval,value,var,xsx);
% generate a new polynomial in var from the value-adic representation
% provided by dval
begin scalar i,d,val,val1,kont;
kont:=0;
val:=dval;
i:=xsx;
if zerop i then <<
% an x**0 term is represented specially;
val1:=negshiftz(remainder(val,value),value);
if not zerop val1 then kont:=d:=val1;
val:=quotient(val-val1,value);
i:=1
>>;
while not zerop val do <<
val1:=negshiftz(remainder(val,value),value);
if not zerop val1 then <<
kont:=gcdn(val1,kont);
d:=var .** i .* val1 .+ d
>>;
val:=quotient(val-val1,value);
i:=1+i
>>;
return quotf(d,kont)
end;
symbolic procedure gen!-poly!-backward(dval,value,var,xsx);
% generate a new polynomial in var from the 1/value-adic representation
% provided by dval
begin scalar i,d,ans,val,val1,kont;
kont:=0;
val:=dval;
% because we are at the 1/value representation
% we need the implicit REVERSE that the two-loop strategy here
% provides;
while not zerop val do <<
val1:=negshiftz(remainder(val,value),value);
d:=val1 . d;
val:=quotient(val-val1,value)
>>;
i:=xsx;
if (zerop i and not zerop car d) then <<
% Handle x**0 term specially;
kont:=ans:=car d;
d:=cdr d;
i:=1
>>;
while d do <<
if not zerop car d then <<
kont:=gcdn(car d,kont);
ans:= var .** i .* car d .+ ans
>>;
d:=cdr d;
i:=i+1
>>;
return quotf(ans,kont)
end;
endmodule;
end;