module unihens; % Univariate case of Hensel code with quadratic growth.
% Author: P. M. A. Moore, 1979.
% Modifications by J.H. Davenport 1988 following a 1985 draft by Abbott,
% Bradford and Davenport.
fluid '(!*linear
!*overshoot
!*overview
!*trfac
alphalist
alphavec
coefftbd
current!-factor!-product
current!-modulus
delfvec
deltam
factor!-level
factor!-trace!-list
factors!-done
factorvec
facvec
fhatvec
hensel!-growth!-size
hensel!-poly
irreducible
m!-image!-variable
modfvec
multivariate!-input!-poly
non!-monic
number!-of!-factors
polyzero
prime!-base
reconstructing!-gcd);
global '(largest!-small!-modulus);
symbolic procedure uhensel!.extend(poly,best!-flist,lclist,p);
% Extend poly=product(factors in best!-flist) mod p even if poly is
% non-monic. Return a list (ok. list of factors) if factors can be
% extended to be correct over the integers, otherwise return a list
% (failed <reason> <reason>).
begin scalar w,k,old!-modulus,alphavec,modular!-flist,factorvec,
modfvec,coefftbd,fcount,fhatvec,deltam,mod!-symm!-flist,
current!-factor!-product,facvec,factors!-done,hensel!-poly;
prime!-base:=p;
old!-modulus:=set!-modulus p;
% timer:=readtime();
number!-of!-factors:=length best!-flist;
w:=expt(lc poly,number!-of!-factors -1);
if lc poly < 0 then errorf list("LC SHOULD NOT BE -VE",poly);
coefftbd:=max(110,p+1,lc poly*get!-coefft!-bound(poly,ldeg poly));
poly:=multf(poly,w);
modular!-flist:=for each ff in best!-flist collect
reduce!-mod!-p ff;
% Modular factors have been multiplied by a constant to
% fix the l.c.'s, so they may be out of range - this
% fixes that.
if not(w=1) then factor!-trace <<
prin2!* "Altered univariate polynomial: "; printsf poly >>;
% Make sure the leading coefft will not cause trouble
% in the hensel construction.
mod!-symm!-flist:=for each ff in modular!-flist collect
make!-modular!-symmetric ff;
if not !*overview then factor!-trace <<
prin2!* "The factors mod "; prin2!* p;
printstr " to start from are:";
fcount:=1;
for each ff in mod!-symm!-flist do <<
prin2!* " f("; prin2!* fcount; prin2!* ")=";
printsf ff; fcount:=iadd1 fcount >>;
terpri!*(nil) >>;
alphalist:=alphas(number!-of!-factors,modular!-flist,1);
% 'magic' polynomials associated with the image factors.
if not !*overview then factor!-trace <<
printstr
"The following modular polynomials are chosen such that:";
terpri();
prin2!* " a(1)*h(1) + ... + a(";
prin2!* number!-of!-factors;
prin2!* ")*h("; prin2!* number!-of!-factors;
prin2!* ") = 1 mod "; printstr p;
terpri();
printstr " where h(i)=(product of all f(j) [see below])/f(i)";
printstr " and degree of a(i) < degree of f(i).";
fcount:=1;
for each a in modular!-flist do <<
prin2!* " a("; prin2!* fcount; prin2!* ")=";
printsf cdr get!-alpha a;
prin2!* " f("; prin2!* fcount; prin2!* ")=";
printsf a;
fcount:=iadd1 fcount >>
>>;
k:=0;
factorvec:=mkvect number!-of!-factors;
modfvec:=mkvect number!-of!-factors;
alphavec:=mkvect number!-of!-factors;
for each modsymmf in mod!-symm!-flist do
<< putv(factorvec,k:=k+1,force!-lc(modsymmf,car lclist));
lclist:=cdr lclist
>>;
k:=0;
for each modfactor in modular!-flist do
<< putv(modfvec,k:=k+1,modfactor);
putv(alphavec,k,cdr get!-alpha modfactor);
>>;
% best!-fvec is now a vector of factors of poly correct
% mod p with true l.c.s forced in.
fhatvec:=mkvect number!-of!-factors;
w:=hensel!-mod!-p(poly,modfvec,factorvec,coefftbd,nil,p);
if car w='overshot then w := uhensel!.extend1(poly,w)
else w := uhensel!.extend2 w;
set!-modulus old!-modulus;
if irreducible then <<
factor!-trace
printstr "Two factors and overshooting means irreducible";
return t >>;
factor!-trace begin scalar k;
k:=0;
printstr "Univariate factors, possibly with adjusted leading";
printstr "coefficients, are:";
for each ww in cdr w do <<
prin2!* " f("; prin2!* (k:=k #+ 1);
prin2!* ")="; printsf ww >>
end;
return if non!-monic then
(car w . primitive!.parts(cdr w,m!-image!-variable,t))
else w
end;
symbolic procedure uhensel!.extend1(poly,w);
begin scalar oklist,badlist,m,r,ff,om,pol;
m:=cadr w; % the modulus.
r:=getv(factorvec,0); % the number of factors.
if r=2 then return (irreducible:=t);
if factors!-done then <<
poly:=hensel!-poly;
for each ww in factors!-done do
poly:=multf(poly,ww) >>;
pol:=poly;
om:=set!-modulus hensel!-growth!-size;
alphalist:=nil;
for i:=r step -1 until 1 do
alphalist:=
(reduce!-mod!-p getv(factorvec,i) . getv(alphavec,i))
. alphalist;
set!-modulus om;
% bring alphalist up to date.
for i:=1:r do <<
ff:=getv(factorvec,i);
if not didntgo(w:=quotf(pol,ff)) then
<< oklist:=ff . oklist; pol:=w>>
else badlist:=(i . ff) . badlist >>;
if null badlist then w:='ok . oklist
else <<
if not !*overview then factor!-trace <<
printstr "Overshot factors are:";
for each f in badlist do <<
prin2!* " f("; prin2!* car f; prin2!* ")=";
printsf cdr f >>
>>;
w:=try!.combining(badlist,pol,m,nil);
if car w='one! bad! factor then begin scalar x;
w:=append(oklist,cdr w);
x:=1;
for each v in w do x:=multf(x,v);
w:='ok . (quotfail(pol,x) . w)
end
else w:='ok . append(oklist,w) >>;
if (not !*linear) and multivariate!-input!-poly then <<
poly:=1;
number!-of!-factors:=0;
for each facc in cdr w do <<
poly:=multf(poly,facc);
number!-of!-factors:=1 #+ number!-of!-factors >>;
% make sure poly is the product of the factors we have,
% we recalculate it this way because we may have the wrong
% lc in old value of poly.
reset!-quadratic!-step!-fluids(poly,cdr w,
number!-of!-factors);
if m=deltam then errorf list("Coefft bound < prime ?",
coefftbd,m);
m:=deltam*deltam;
while m<largest!-small!-modulus do <<
quadratic!-step(m,number!-of!-factors);
m:=m*deltam >>;
hensel!-growth!-size:=deltam;
om:=set!-modulus hensel!-growth!-size;
alphalist:=nil;
for i:=number!-of!-factors step -1 until 1 do
alphalist:=
(reduce!-mod!-p getv(factorvec,i) . getv(alphavec,i))
. alphalist;
set!-modulus om >>;
return w
end;
symbolic procedure uhensel!.extend2 w;
begin scalar r,faclist,om;
r:=getv(factorvec,0); % no of factors.
om:=set!-modulus hensel!-growth!-size;
alphalist:=nil;
for i:=r step -1 until 1 do
alphalist:=(reduce!-mod!-p getv(factorvec,i) . getv(alphavec,i))
. alphalist;
set!-modulus om;
% bring alphalist up to date.
for i:=r step -1 until 1 do
faclist:=getv(factorvec,i) . faclist;
return (car w . faclist)
end;
symbolic procedure get!-coefft!-bound(poly,ddeg);
% This uses Mignotte's bound which is minimal I believe.
% NB. poly had better be univariate as bound only valid for this.
binomial!-coefft(ddeg/2,ddeg/4) * root!-squares(poly,0);
symbolic procedure binomial!-coefft(n,r);
if n<r then nil
else if n=r then 1
else if r=1 then n
else begin scalar n!-c!-r,b;
n!-c!-r:=1;
b:=min(r,n-r);
for i:=1:b do
n!-c!-r:=(n!-c!-r * (n - i + 1)) / i;
return n!-c!-r
end;
symbolic procedure find!-alphas!-in!-a!-ring(n,mflist,fhatlist,gamma);
% Find the alphas (as below) given that the modulus may not be prime
% but is a prime power.
begin scalar gg,m,ppow,i,gg!-mod!-p,modflist,wvec,alpha,alphazeros,w;
if null prime!-base then errorf
list("Prime base not set for finding alphas",
current!-modulus,n,mflist);
m:=set!-modulus prime!-base;
modflist:= if m=prime!-base then mflist
else for each fthing in mflist collect
reduce!-mod!-p !*mod2f fthing;
alphalist:=alphas(n,modflist,gamma);
if m=prime!-base then <<
set!-modulus m;
return alphalist >>;
i:=0;
alphazeros:=mkvect n;
wvec:=mkvect n;
for each modfthing in modflist do <<
putv(modfvec,i:=iadd1 i,modfthing);
putv(alphavec,i,!*f2mod(alpha:=cdr get!-alpha modfthing));
putv(alphazeros,i,alpha);
putv(wvec,i,alpha);
putv(fhatvec,i,car fhatlist);
fhatlist:=cdr fhatlist >>;
gg:=gamma;
ppow:=prime!-base;
while ppow<m do <<
set!-modulus m;
gg:=!*f2mod quotfail(!*mod2f difference!-mod!-p(gg,
form!-sum!-and!-product!-mod!-m(wvec,fhatvec,n)),prime!-base);
set!-modulus prime!-base;
gg!-mod!-p:=reduce!-mod!-p !*mod2f gg;
for k:=1:n do <<
putv(wvec,k,w:=remainder!-mod!-p(
times!-mod!-p(getv(alphazeros,k),gg!-mod!-p),
getv(modfvec,k)));
putv(alphavec,k,addf(getv(alphavec,k),multf(!*mod2f w,ppow)))>>;
ppow:=ppow*prime!-base >>;
set!-modulus m;
i:=0;
return (for each fthing in mflist collect
(fthing . !*f2mod getv(alphavec,i:=iadd1 i)))
end;
symbolic procedure alphas(n,flist,gamma);
% Finds alpha,beta,delta,... wrt factors f(i) in flist s.t.
% alpha*g(1) + beta*g(2) + delta*g(3) + ... = gamma mod p,
% where g(i)=product(all the f(j) except f(i) itself).
% (cf. xgcd!-mod!-p below). n is number of factors in flist.
if n=1 then list(car flist . gamma)
else begin scalar k,w,f1,f2,i,gamma1,gamma2;
k:=n/2;
f1:=1; f2:=1;
i:=1;
for each f in flist do
<< if i>k then f2:=times!-mod!-p(f,f2)
else f1:=times!-mod!-p(f,f1);
i:=i+1 >>;
w:=xgcd!-mod!-p(f1,f2,1,polyzero,polyzero,1);
if atom w then
return 'factors! not! coprime;
gamma1:=remainder!-mod!-p(times!-mod!-p(cdr w,gamma),f1);
gamma2:=remainder!-mod!-p(times!-mod!-p(car w,gamma),f2);
i:=1; f1:=nil; f2:=nil;
for each f in flist do
<< if i>k then f2:=f . f2
else f1:=f . f1;
i:=i+1 >>;
return append(
alphas(k,f1,gamma1),
alphas(n-k,f2,gamma2))
end;
symbolic procedure xgcd!-mod!-p(a,b,x1,y1,x2,y2);
% Finds alpha and beta s.t. alpha*a+beta*b=1.
% Returns alpha . beta or nil if a and b are not coprime.
if null b then nil
else if domainp b then begin
b:=modular!-reciprocal b;
x2:=multiply!-by!-constant!-mod!-p(x2,b);
y2:=multiply!-by!-constant!-mod!-p(y2,b);
return x2 . y2 end
else begin scalar q;
q:=quotient!-mod!-p(a,b); % Truncated quotient here.
return xgcd!-mod!-p(b,difference!-mod!-p(a,times!-mod!-p(b,q)),
x2,y2,
difference!-mod!-p(x1,times!-mod!-p(x2,q)),
difference!-mod!-p(y1,times!-mod!-p(y2,q)))
end;
symbolic procedure hensel!-mod!-p(poly,mvec,fvec,cbd,vset,p);
% Hensel construction building up in powers of p.
% Given that poly=product(factors in factorvec) mod p, find the full
% factors over the integers. Mvec contains the univariate factors mod p
% while fvec contains our best knowledge of the factors to date.
% Fvec includes leading coeffts (and in multivariate case possibly other
% coeffts) of the factors. return a list whose first element is a flag
% with one of the following values:
% ok construction worked, the cdr of the result is a list of
% the correct factors.
% failed inputs must have been incorrect
% overshot factors are correct mod some power of p (say p**m),
% but are not correct over the integers.
% result is (overshot,p**m,list of factors so far).
begin scalar w,u0,delfvec,old!.mod,res,m;
u0:=initialize!-hensel(number!-of!-factors,p,poly,mvec,fvec,cbd);
% u0 contains the product (over integers) of factors mod p.
m := p;
old!.mod := set!-modulus nil;
if number!-of!-factors=1
then <<putv(fvec,1,current!-factor!-product:=poly);
% Added JHD 28.9.87
return hensel!-exit(m,old!.mod,p,vset,w)>>;
% only one factor to grow! but need to go this deep to
% construct the alphas and set things up for the
% multivariate growth which may follow.
hensel!-msg1(p,u0);
old!.mod:=set!-modulus p;
res:=addf(hensel!-poly,negf u0);
% calculate the residue. from now on this is always
% kept in res.
m:=p;
% measure of how far we have built up factors - at this
% stage we know the constant terms mod p in the factors.
a: if polyzerop res then return hensel!-exit(m,old!.mod,p,vset,w);
if (m/2)>coefftbd then
<<
% we started with a false split of the image so some
% of the factors we have built up must amalgamate in
% the complete factorization.
if !*overshoot then <<
prin2 if null vset then "Univariate " else "Multivariate ";
prin2t "coefft bound overshoot" >>;
if not !*overview then
factor!-trace printstr "We have overshot the coefficient bound";
return hensel!-exit(m,old!.mod,p,vset,'overshot)>>;
res:=quotfail(res,deltam);
% next term in residue.
if not !*overview then factor!-trace <<
prin2!* "Residue divided by "; prin2!* m; prin2!* " is ";
printsf res >>;
if (not !*linear) and null vset
and m<=largest!-small!-modulus and m>p then
quadratic!-step(m,number!-of!-factors);
w:=reduce!-mod!-p res;
if not !*overview then factor!-trace <<
prin2!* "Next term in residue to kill is:";
prinsf w; prin2!* " which is of size ";
printsf (deltam*m);
>>;
solve!-for!-corrections(w,fhatvec,modfvec,delfvec,vset);
% delfvec is vector of next correction terms to factors.
make!-vec!-modular!-symmetric(delfvec,number!-of!-factors);
if not !*overview then factor!-trace <<
printstr "Correction terms are:";
w:=1;
for i:=1:number!-of!-factors do <<
prin2!* " To f("; prin2!* w; prin2!* "): ";
printsf multf(m,getv(delfvec,i));
w:=iadd1 w >>;
>>;
w:=terms!-done(factorvec,delfvec,m);
res:=addf(res,negf w);
% subtract out the terms generated by these corrections
% from the residue.
current!-factor!-product:=
addf(current!-factor!-product,multf(m,w));
% add in the correction terms to give new factor product.
for i:=1:number!-of!-factors do
putv(factorvec,i,
addf(getv(factorvec,i),multf(getv(delfvec,i),m)));
% add the corrections into the factors.
if not !*overview then factor!-trace <<
printstr " giving new factors as:";
w:=1;
for i:=1:number!-of!-factors do <<
prin2!* " f("; prin2!* w; prin2!* ")=";
printsf getv(factorvec,i); w:=iadd1 w >>
>>;
m:=m*deltam;
if not polyzerop res and null vset and
not reconstructing!-gcd then
begin scalar j,u,fac;
j:=0;
while (j:=j #+ 1)<=number!-of!-factors do
% IF NULL GETV(DELFVEC,J) AND
% - Try dividing out every time for now.
if not didntgo
(u:=quotf(hensel!-poly,fac:=getv(factorvec,j))) then <<
hensel!-poly:=u;
res:=adjust!-growth(fac,j,m);
j:=number!-of!-factors >>
end;
go to a
end;
symbolic procedure hensel!-exit(m,old!.mod,p,vset,w);
begin
if factors!-done then <<
if not(w='overshot) then m:=p*p;
set!-hensel!-fluids!-back p >>;
if (not (w='overshot)) and null vset
and (not !*linear) and multivariate!-input!-poly then
while m<largest!-small!-modulus do <<
if not(m=deltam) then quadratic!-step(m,number!-of!-factors);
m:=m*deltam >>;
% set up the alphas etc so that multivariate growth can
% use a Hensel growth size of about word size.
set!-modulus old!.mod;
% reset the old modulus.
hensel!-growth!-size:=deltam;
putv(factorvec,0,number!-of!-factors);
return
if w='overshot then list('overshot,m,factorvec)
else 'ok . factorvec
end;
symbolic procedure hensel!-msg1(p,u0);
begin scalar w;
factor!-trace <<
printstr
"We are now ready to use the Hensel construction to grow";
prin2!* "in powers of "; printstr current!-modulus;
if not !*overview then <<prin2!* "Polynomial to factor (=U): ";
printsf hensel!-poly>>;
prin2!* "Initial factors mod "; prin2!* p;
printstr " with some correct coefficients:";
w:=1;
for i:=1:number!-of!-factors do <<
prin2!* " f("; prin2!* w; prin2!* ")=";
printsf getv(factorvec,i); w:=iadd1 w >>;
if not !*overview then << prin2!* "Coefficient bound = ";
prin2!* coefftbd;
terpri!*(nil);
prin2!* "The product of factors over the integers is ";
printsf u0;
printstr "In each step below, the residue is U - (product of the";
printstr
"factors as far as we know them). The correction to each";
printstr "factor, f(i), is (a(i)*v) mod f0(i) where f0(i) is";
prin2!* "f(i) mod "; prin2!* p;
printstr "(ie. the f(i) used in calculating the a(i))"
>>>>
end;
symbolic procedure initialize!-hensel(r,p,poly,mvec,fvec,cbd);
% Set up the vectors and initialize the fluids.
begin scalar u0;
delfvec:=mkvect r;
facvec:=mkvect r;
hensel!-poly:=poly;
modfvec:=mvec;
factorvec:=fvec;
coefftbd:=cbd;
factors!-done:=nil;
deltam:=p;
u0:=1;
for i:=1:r do u0:=multf(getv(factorvec,i),u0);
current!-factor!-product:=u0;
return u0
end;
% symbolic procedure reset!-quadratic!-step!-fluids(poly,faclist,n);
% begin scalar i,om,modf;
% current!-factor!-product:=poly;
% om:=set!-modulus hensel!-growth!-size;
% i:=0;
% for each fac in faclist do <<
% putv(factorvec,i:=iadd1 i,fac);
% putv(modfvec,i,modf:=reduce!-mod!-p fac);
% putv(alphavec,i,cdr get!-alpha modf) >>;
% for i:=1:n do <<
% prin2 "F("; % prin2 i; % prin2 ") = ";
% printsf getv(factorvec,i);
% prin2 "F("; % prin2 i; % prin2 ") MOD P = ";
% printsf getv(modfvec,i);
% prin2 "A("; % prin2 i; % prin2 ") = ";
% printsf getv(alphavec,i) >>;
% set!-modulus om
% end;
symbolic procedure reset!-quadratic!-step!-fluids(poly,faclist,n);
begin scalar i,om,facpairlist,cfp!-mod!-p,fhatlist;
current!-factor!-product:=poly;
om:=set!-modulus hensel!-growth!-size;
cfp!-mod!-p:=reduce!-mod!-p current!-factor!-product;
i:=0;
facpairlist:=for each fac in faclist collect <<
i:= i #+ 1;
(fac . reduce!-mod!-p fac) >>;
fhatlist:=for each facc in facpairlist collect
quotfail!-mod!-p(cfp!-mod!-p,cdr facc);
if factors!-done then alphalist:=
find!-alphas!-in!-a!-ring(i,
for each facpr in facpairlist collect cdr facpr,
fhatlist,1);
% a bug has surfaced such that the alphas get out of step.
% In this case so recalculate them to stop the error for now.
i:=0;
for each facpair in facpairlist do <<
putv(factorvec,i:=iadd1 i,car facpair);
putv(modfvec,i,cdr facpair);
putv(alphavec,i,cdr get!-alpha cdr facpair) >>;
% for i:=1:n do <<
% prin2 "f("; % prin2 i; % prin2 ") = ";
% printsf getv(factorvec,i);
% prin2 "f("; % prin2 i; % prin2 ") mod p = ";
% printsf getv(modfvec,i);
% prin2 "a("; % prin2 i; % prin2 ") = ";
% printsf getv(alphavec,i) >>;
set!-modulus om
end;
symbolic procedure quadratic!-step(m,r);
% Code for adjusting the hensel variables to take quadratic steps in
% the growing process.
begin scalar w,s,cfp!-mod!-p;
set!-modulus m;
cfp!-mod!-p:=reduce!-mod!-p current!-factor!-product;
for i:=1:r do putv(facvec,i,reduce!-mod!-p getv(factorvec,i));
for i:=1:r do putv(fhatvec,i,
quotfail!-mod!-p(cfp!-mod!-p,getv(facvec,i)));
w:=form!-sum!-and!-product!-mod!-m(alphavec,fhatvec,r);
w:=!*mod2f plus!-mod!-p(1,minus!-mod!-p w);
s:=quotfail(w,deltam);
set!-modulus deltam;
s:=!*f2mod s;
% Boxes S up to look like a poly mod deltam.
for i:=1:r do <<
w:=remainder!-mod!-p(times!-mod!-p(s,getv(alphavec,i)),
getv(modfvec,i));
putv(alphavec,i,
addf(!*mod2f getv(alphavec,i),multf(!*mod2f w,deltam))) >>;
s:=modfvec;
modfvec:=facvec;
facvec:=s;
deltam:=m;
% this is our new growth rate.
set!-modulus deltam;
for i:=1:r do <<
putv(facvec,i,"RUBBISH");
% we will want to overwrite facvec next time so we
% had better point it to the old (no longer needed)
% modvec. Also mark it as containing rubbish for safety.
putv(alphavec,i,!*f2mod getv(alphavec,i)) >>;
% Make sure the alphas are boxed up as being mod new deltam.
if not !*overview then factor!-trace <<
printstr "The new modular polynomials are chosen such that:";
terpri();
prin2!* " a(1)*h(1) + ... + a(";
prin2!* r;
prin2!* ")*h("; prin2!* r;
prin2!* ") = 1 mod "; printstr m;
terpri();
printstr " where h(i)=(product of all f(j) [see below])/f(i)";
printstr " and degree of a(i) < degree of f(i).";
for i:=1:r do <<
prin2!* " a("; prin2!* i; prin2!* ")=";
printsf getv(alphavec,i);
prin2!* " f("; prin2!* i; prin2!* ")=";
printsf getv(modfvec,i) >>
>>
end;
symbolic procedure terms!-done(fvec,delfvec,m);
begin scalar flist,delflist;
for i:=1:number!-of!-factors do <<
flist:=getv(fvec,i) . flist;
delflist:=getv(delfvec,i) . delflist >>;
return terms!.done(number!-of!-factors,flist,delflist,
number!-of!-factors,m)
end;
symbolic procedure terms!.done(n,flist,delflist,r,m);
if n=1 then (car flist) . (car delflist)
else begin scalar k,i,f1,f2,delf1,delf2;
k:=n/2; i:=1;
for each f in flist do
<< if i>k then f2:=(f . f2)
else f1:=(f . f1);
i:=i+1 >>;
i:=1;
for each delf in delflist do
<< if i>k then delf2:=(delf . delf2)
else delf1:=(delf . delf1);
i:=i+1 >>;
f1:=terms!.done(k,f1,delf1,r,m);
delf1:=cdr f1; f1:=car f1;
f2:=terms!.done(n-k,f2,delf2,r,m);
delf2:=cdr f2; f2:=car f2;
delf1:=
addf(addf(
multf(f1,delf2),
multf(f2,delf1)),
multf(multf(delf1,m),delf2));
if n=r then return delf1;
return (multf(f1,f2) . delf1)
end;
symbolic procedure try!.combining(l,poly,m,sofar);
try!.combining1(l,poly,m,sofar,2);
% The following code is not optimal. If we find a combination of k
% factors, we will re-rty all the previous combinations of k factors
% already tried.
% This is not frightfully serious, since in practice most such
% combinations will have had something in common with the set found, and
% so won't re-appear. This is definitely better than the previous
% version, which re-tried all combinations. JHD 14.1.88.
symbolic procedure try!.combining1(l,poly,m,sofar,k);
% l is a list of factors, f(i), s.t. (product of the f(i) mod m) = poly
% but no f(i) divides poly over the integers. We find the combinations
% of the f(i) that yield the true factors of poly over the integers.
% Sofar is a list of these factors found so far.
% start combining them K at a time
if poly=1 then
if null l then sofar
else errorf(list("TOO MANY BAD FACTORS:",l))
else begin scalar n,res,ff,v,w,w1,combined!.factors,ll,lcfinv,oldmod;
n:=length l;
if n=1 then
if ldeg car l > (ldeg poly)/2 then
return ('one! bad! factor . sofar)
else errorf(list("ONE BAD FACTOR DOES NOT FIT:",l));
if n=2 or n=3 then <<
w:=lc cdar l; % The LC of all the factors is the same.
while not (w=lc poly) do poly:=quotfail(poly,w);
% poly's LC may be a higher power of w than we want
% and we must return a result with the same
% LC as each of the combined factors.
if not !*overview then factor!-trace <<
printstr "We combine:";
for each lf in l do printsf cdr lf;
prin2!* " mod "; prin2!* m;
printstr " to give correct factor:";
printsf poly >>;
combine!.alphas(l,t);
return (poly . sofar) >>;
ll:=for each ff in l collect (cdr ff . car ff);
% K := 2; % K is now an argument
oldmod := set!-general!-modulus m;
lcfinv := general!-modular!-reciprocal lc cdar l;
set!-general!-modulus oldmod;
loop1:
if k > n/2 then go to exit;
w:=koutof(k,if 2*k=n then cdr l else l,nil);
% We needn't try a combination and its complement
while w and
(v:=factor!-trialdiv(poly,car w,m,ll,lcfinv))='didntgo do
<< w:=cdr w;
while w and
((car w = '!*lazyadjoin) or (car w = '!*lazykoutof)) do
if car w= '!*lazyadjoin then
w:=lazy!-adjoin(cadr w,caddr w,cadr cddr w)
else w:=koutof(cadr w,caddr w,cadr cddr w)
>>;
if not(v='didntgo) then <<
ff:=car v; v:=cdr v;
if not !*overview then factor!-trace <<
printstr "We combine:";
for each a in car w do printsf a;
prin2!* " mod "; prin2!* m;
printstr " to give correct factor:";
printsf ff >>;
for each a in car w do <<
w1:=l;
while not (a = cdar w1) do w1:=cdr w1;
combined!.factors:=car w1 . combined!.factors;
l:=delete(car w1,l) >>;
combine!.alphas(combined!.factors,t);
% Now combine the rest, starting with k-tuples
res:=try!.combining1(l,v,m,ff . sofar,k);
go to exit>>;
k := k + 1;
go to loop1;
exit:
if res then return res
else <<
w:=lc cdar l; % The LC of all the factors is the same.
while not (w=lc poly) do poly:=quotfail(poly,w);
% poly's LC may be a higher power of w than we want
% and we must return a result with the same
% LC as each of the combined factors.
if not !*overview then factor!-trace <<
printstr "We combine:";
for each ff in l do printsf cdr ff;
prin2!* " mod "; prin2!* m;
printstr " to give correct factor:";
printsf poly >>;
combine!.alphas(l,t);
return (poly . sofar) >>
end;
symbolic procedure koutof(k,l,sofar);
% Produces all permutations of length k from list l accumulating them
% in sofar as we go. We use lazy evaluation in that this results in
% a permutation dotted with:
% ( '!*lazy . (argument for eval) )
% except when k=1 when the permutations are explicitly given.
if k=1 then append(
for each f in l collect list cdr f,sofar)
else if k>length l then sofar
else <<
while eqcar(l,'!*lazyadjoin) or eqcar(l,'!*lazykoutof) do
if car l='!*lazyadjoin then
l := lazy!-adjoin(cadr l,caddr l,cadr cddr l)
else l := koutof(cadr l,caddr l,cadr cddr l);
if k=length l then
(for each ll in l collect cdr ll ) . sofar
else koutof(k,cdr l,
list('!*lazyadjoin,cdar l,
list('!*lazykoutof,(k-1),cdr l,nil),
sofar)) >>;
symbolic procedure lazy!-adjoin(item,l,tail);
% Dots item with each element in l using lazy evaluation on l.
% If l is null tail results.
<< while eqcar(l,'!*lazyadjoin) or eqcar(l,'!*lazykoutof) do
if car l ='!*lazyadjoin then
l:=lazy!-adjoin(cadr l,caddr l,cadr cddr l)
else l:=koutof(cadr l,caddr l,cadr cddr l);
if null l then tail
else (item . car l) .
if null cdr l then tail
else list('!*lazyadjoin,item,cdr l,tail) >>;
symbolic procedure factor!-trialdiv(poly,flist,m,llist,lcfinv);
% Combines the factors in FLIST mod M and test divides the result
% into POLY (over integers) to see if it goes. If it doesn't
% then DIDNTGO is returned, else the pair (D . Q) is
% returned where Q is the quotient obtained and D is the product
% of the factors mod M;
% Abbott,J.A., Bradford,R.J. & Davenport,J.H.,
% A Remark on Factorisation.
% SIGSAM Bulletin 19(1985) 2, pp. 31-33, 37.
if polyzerop poly then errorf "Test dividing into zero?"
else begin scalar d,q,tcpoly,tcoeff,x,oldmod,w,poly1,try1;
factor!-trace <<
prin2!* "We combine factors ";
for each ff in flist do <<
w:=assoc(ff,llist);
prin2!* "f(";
prin2!* cdr w;
prin2!* "), " >> ;
prin2!* "and try dividing : " >>;
x := mvar poly;
tcpoly :=trailing!.coefft(poly,x);
tcoeff := trailing!.coefft(car flist,x);
oldmod := set!-general!-modulus m;
for each fac in cdr flist do
tcoeff := general!-modular!-times(
general!-modular!-times(tcoeff,lcfinv),
trailing!.coefft(fac,x));
if not zerop remainder(tcpoly,
w:=general!-make!-modular!-symmetric tcoeff) then <<
factor!-trace printstr " it didn't go (tc test)";
set!-general!-modulus oldmod;
% if not(w = trailing!.coefft(car COMBINE(FLIST,M,LLIST,lcfinv),x))
% then <<
% printstr "incompatibility: we have";
% prin2!* w;
% printstr "which should be the trailing coefficient of:";
% prin2!* car combine(flist,m,llist,lcfinv) >>;
return 'didntgo >>;
% it has passed the tc test - now try evaluating at 1;
poly1 := eval!-at!-1 poly;
try1 := eval!-at!-1 car flist;
for each fac in cdr flist do
try1 := general!-modular!-times(
general!-modular!-times(try1,lcfinv),eval!-at!-1 fac);
if (zerop try1 and not zerop poly1) or
not zerop remainder(poly1,general!-make!-modular!-symmetric try1)
then <<
factor!-trace printstr " it didn't go (test at 1)";
set!-general!-modulus oldmod;
return 'didntgo >>;
% it has passed both tests - work out longhand;
set!-general!-modulus oldmod;
d:=combine(flist,m,llist,lcfinv);
if didntgo(q:=quotf(poly,car d)) then <<
factor!-trace printstr " it didn't go (division fail)";
return 'didntgo >>
else <<
factor!-trace printstr " it worked !";
return (car d . quotf(q,cdr d)) >>
end;
symbolic procedure eval!-at!-1 f;
% f a univariate standard form over Z with f(0) neq 0;
% return the integer f(1);
if atom f then f else (lc f) + eval!-at!-1(red f);
symbolic procedure combine(flist,m,l,lcfinv);
% Multiply factors in flist mod m.
% L is a list of the factors for use in FACTOR!-TRACE.
begin scalar om,res,lcf,lcfprod;
lcf := lc car flist; % all leading coeffts should be the same.
lcfprod := 1;
% This is one of only two places in the entire factorizer where
% it is ever necessary to use a modulus larger than word-size.
if m>largest!-small!-modulus then <<
om:=set!-general!-modulus m;
% lcfinv := general!-modular!-reciprocal lcf; Done once and for all
res:=general!-reduce!-mod!-p car flist;
for each ff in cdr flist do <<
if not(lcf=lc ff) then errorf "BAD LC IN FLIST";
res:=general!-times!-mod!-p(
general!-times!-mod!-p(lcfinv,
general!-reduce!-mod!-p ff),res);
lcfprod := lcfprod*lcf >>;
res:=general!-make!-modular!-symmetric res;
set!-modulus om;
return (res . lcfprod) >>
else <<
om:=set!-modulus m;
lcfinv := modular!-reciprocal lcf;
res:=reduce!-mod!-p car flist;
for each ff in cdr flist do <<
if not(lcf=lc ff) then errorf "BAD LC IN FLIST";
res:=times!-mod!-p(times!-mod!-p(lcfinv,reduce!-mod!-p ff),res);
lcfprod := lcfprod*lcf >>;
res:=make!-modular!-symmetric res;
set!-modulus om;
return (res . lcfprod) >>
end;
symbolic procedure combine!.alphas(flist,fixlcs);
% Combine the alphas associated with each of these factors to
% give the one alpha for their combination.
begin scalar f1,a1,ff,aa,oldm,lcfac,lcfinv,saveflist;
oldm:=set!-modulus hensel!-growth!-size;
flist:=for each fac in flist collect <<
saveflist:= (reduce!-mod!-p cdr fac) . saveflist;
(car fac) . car saveflist >>;
if fixlcs then <<
lcfinv:=modular!-reciprocal lc cdar flist;
lcfac:=modular!-expt(lc cdar flist,sub1 length flist)
>>
else << lcfinv:=1; lcfac:=1 >>;
% If FIXLCS is set then we have combined n factors
% (each with the same l.c.) to give one and we only need one
% l.c. in the result, we have divided the combination by
% lc**(n-1) and we must be sure to do the same for the
% alphas.
ff:=cdar flist;
aa:=cdr get!-alpha ff;
flist:=cdr flist;
while flist do <<
f1:=cdar flist;
a1:=cdr get!-alpha f1;
flist:=cdr flist;
aa:=plus!-mod!-p(times!-mod!-p(aa,f1),times!-mod!-p(a1,ff));
ff:=times!-mod!-p(ff,f1)
>>;
for each a in alphalist do
if not member(car a,saveflist) then
flist:=(car a . times!-mod!-p(cdr a,lcfac)) . flist;
alphalist:=(quotient!-mod!-p(ff, lcfac) . aa) . flist;
set!-modulus oldm
end;
% The following code is for dividing out factors in the middle
% of the Hensel construction and adjusting all the associated
% variables that go with it.
symbolic procedure adjust!-growth(facdone,k,m);
% One factor (at least) divides out so we can reconfigure the
% problem for Hensel constrn giving a smaller growth and hopefully
% reducing the coefficient bound considerably.
begin scalar w,u,bound!-scale,modflist,factorlist,fhatlist,
modfdone,b;
factorlist:=vec2list!-without!-k(factorvec,k);
modflist:=vec2list!-without!-k(modfvec,k);
fhatlist:=vec2list!-without!-k(fhatvec,k);
w:=number!-of!-factors;
modfdone:=getv(modfvec,k);
top:
factors!-done:=facdone . factors!-done;
if (number!-of!-factors:=number!-of!-factors #- 1)=1 then <<
factors!-done:=hensel!-poly . factors!-done;
number!-of!-factors:=0;
hensel!-poly:=1;
if not !*overview then factor!-trace <<
printstr " All factors found:";
for each fd in factors!-done do printsf fd >>;
return polyzero >>;
fhatlist:=for each fhat in fhatlist collect
quotfail!-mod!-p(if null fhat then polyzero else fhat,modfdone);
u:=comfac facdone; % Take contents and prim. parts.
if car u then
errorf(list("Factor divisible by main variable: ",facdone,car u));
facdone:=quotfail(facdone,cdr u);
bound!-scale:=cdr u;
if not((b:=lc facdone)=1) then begin scalar b!-inv,old!-m;
hensel!-poly:=quotfail(hensel!-poly,b**number!-of!-factors);
b!-inv:=modular!-reciprocal modular!-number b;
modflist:=for each modf in modflist collect
times!-mod!-p(b!-inv,modf);
% This is one of only two places in the entire factorizer where
% it is ever necessary to use a modulus larger than word-size.
if m>largest!-small!-modulus then <<
old!-m:=set!-general!-modulus m;
factorlist:=for each facc in factorlist collect
adjoin!-term(lpow facc,quotfail(lc facc,b),
general!-make!-modular!-symmetric(
general!-times!-mod!-p(
general!-modular!-reciprocal general!-modular!-number b,
general!-reduce!-mod!-p red facc))) >>
else <<
old!-m:=set!-modulus m;
factorlist:=for each facc in factorlist collect
adjoin!-term(lpow facc,quotfail(lc facc,b),
make!-modular!-symmetric(
times!-mod!-p(modular!-reciprocal modular!-number b,
reduce!-mod!-p red facc))) >>;
% We must be careful not to destroy the information
% that we have about the leading coefft.
set!-modulus old!-m;
fhatlist:=for each fhat in fhatlist collect
times!-mod!-p(
modular!-expt(b!-inv,number!-of!-factors #- 1),fhat)
end;
try!-another!-factor:
if (w:=w #- 1)>0 then
if not didntgo
(u:=quotf(hensel!-poly,facdone:=car factorlist)) then <<
hensel!-poly:=u;
factorlist:=cdr factorlist;
modfdone:=car modflist;
modflist:=cdr modflist;
fhatlist:=cdr fhatlist;
goto top >>
else <<
factorlist:=append(cdr factorlist,list car factorlist);
modflist:=append(cdr modflist,list car modflist);
fhatlist:=append(cdr fhatlist,list car fhatlist);
goto try!-another!-factor >>;
set!-fluids!-for!-newhensel(factorlist,fhatlist,modflist);
bound!-scale:=
bound!-scale * get!-coefft!-bound(
quotfail(hensel!-poly,bound!-scale**(number!-of!-factors #- 1)),
ldeg hensel!-poly);
% We expect the new coefficient bound to be smaller, but on
% dividing out a factor our polynomial's height may have grown
% more than enough to compensate in the bound formula for
% the drop in degree. Anyway, the bound we computed last time
% will still be valid, so let's stick with the smaller.
if bound!-scale < coefftbd then coefftbd := bound!-scale;
w:=quotfail(addf(hensel!-poly,negf current!-factor!-product),
m/deltam);
if not !*overview then factor!-trace <<
printstr " Factors found to be correct:";
for each fd in factors!-done do
printsf fd;
printstr "Remaining factors are:";
printvec(" f(",number!-of!-factors,") = ",factorvec);
prin2!* "New coefficient bound is "; printstr coefftbd;
prin2!* " and the residue is now "; printsf w >>;
return w
end;
symbolic procedure vec2list!-without!-k(v,k);
% Turn a vector into a list leaving out Kth element.
begin scalar w;
for i:=1:number!-of!-factors do
if not(i=k) then w:=getv(v,i) . w;
return w
end;
symbolic procedure set!-fluids!-for!-newhensel(flist,fhatlist,modflist);
<< current!-factor!-product:=1;
alphalist:=
find!-alphas!-in!-a!-ring(number!-of!-factors,modflist,fhatlist,1);
for i:=number!-of!-factors step -1 until 1 do <<
putv(factorvec,i,car flist);
putv(modfvec,i,car modflist);
putv(fhatvec,i,car fhatlist);
putv(alphavec,i,cdr get!-alpha car modflist);
current!-factor!-product:=multf(car flist,current!-factor!-product);
flist:=cdr flist;
modflist:=cdr modflist;
fhatlist:=cdr fhatlist >>
>>;
symbolic procedure set!-hensel!-fluids!-back p;
% After the Hensel growth we must be careful to set back any fluids
% that have been changed when we divided out a factor in the middle
% of growing. Since calculating the alphas involves modular division
% we cannot do it mod DELTAM which is generally a non-trivial power of
% P (prime). So we calculate them mod P and if necessary we can do a
% few quadratic growth steps later.
begin scalar n,fd,modflist,fullf,modf;
set!-modulus p;
deltam:=p;
n:=number!-of!-factors #+ length (fd:=factors!-done);
current!-factor!-product:=hensel!-poly;
for i:=(number!-of!-factors #+ 1):n do <<
putv(factorvec,i,fullf:=car fd);
putv(modfvec,i,modf:=reduce!-mod!-p fullf);
current!-factor!-product:=multf(fullf,current!-factor!-product);
modflist:=modf . modflist;
fd:=cdr fd >>;
for i:=1:number!-of!-factors do <<
modf:=reduce!-mod!-p !*mod2f getv(modfvec,i);
% need to 'unbox' a modpoly before reducing it mod p as we
% know that the input modpoly is wrt a larger modulus
% (otherwise this would be a stupid thing to do anyway!)
% and so we are just pretending it is a full poly.
modflist:=modf . modflist;
putv(modfvec,i,modf) >>;
alphalist:=alphas(n,modflist,1);
for i:=1:n do putv(alphavec,i,cdr get!-alpha getv(modfvec,i));
number!-of!-factors:=n
end;
endmodule;
end;