module bigmodp; % Modular polynomial arithmetic where the modulus may
% be a bignum.
% Authors: A. C. Norman and P. M. A. Moore, 1981;
fluid '(current!-modulus modulus!/2);
symbolic procedure general!-plus!-mod!-p(a,b);
% form the sum of the two polynomials a and b
% working over the ground domain defined by the routines
% general!-modular!-plus, general!-modular!-times etc. the inputs to
% this routine are assumed to have coefficients already
% in the required domain;
if null a then b
else if null b then a
else if isdomain a then
if isdomain b then !*num2f general!-modular!-plus(a,b)
else (lt b) .+ general!-plus!-mod!-p(a,red b)
else if isdomain b then (lt a) .+ general!-plus!-mod!-p(red a,b)
else if lpow a = lpow b then
adjoin!-term(lpow a,
general!-plus!-mod!-p(lc a,lc b),
general!-plus!-mod!-p(red a,red b))
else if comes!-before(lpow a,lpow b) then
(lt a) .+ general!-plus!-mod!-p(red a,b)
else (lt b) .+ general!-plus!-mod!-p(a,red b);
symbolic procedure general!-times!-mod!-p(a,b);
if (null a) or (null b) then nil
else if isdomain a then gen!-mult!-by!-const!-mod!-p(b,a)
else if isdomain b then gen!-mult!-by!-const!-mod!-p(a,b)
else if mvar a=mvar b then general!-plus!-mod!-p(
general!-plus!-mod!-p(general!-times!-term!-mod!-p(lt a,b),
general!-times!-term!-mod!-p(lt b,red a)),
general!-times!-mod!-p(red a,red b))
else if ordop(mvar a,mvar b) then
adjoin!-term(lpow a,general!-times!-mod!-p(lc a,b),
general!-times!-mod!-p(red a,b))
else adjoin!-term(lpow b,
general!-times!-mod!-p(a,lc b),general!-times!-mod!-p(a,red b));
symbolic procedure general!-times!-term!-mod!-p(term,b);
%multiply the given polynomial by the given term;
if null b then nil
else if isdomain b then
adjoin!-term(tpow term,
gen!-mult!-by!-const!-mod!-p(tc term,b),nil)
else if tvar term=mvar b then
adjoin!-term(mksp(tvar term,iplus(tdeg term,ldeg b)),
general!-times!-mod!-p(tc term,lc b),
general!-times!-term!-mod!-p(term,red b))
else if ordop(tvar term,mvar b) then
adjoin!-term(tpow term,general!-times!-mod!-p(tc term,b),nil)
else adjoin!-term(lpow b,
general!-times!-term!-mod!-p(term,lc b),
general!-times!-term!-mod!-p(term,red b));
symbolic procedure gen!-mult!-by!-const!-mod!-p(a,n);
% multiply the polynomial a by the constant n;
if null a then nil
else if n=1 then a
else if isdomain a then !*num2f general!-modular!-times(a,n)
else adjoin!-term(lpow a,gen!-mult!-by!-const!-mod!-p(lc a,n),
gen!-mult!-by!-const!-mod!-p(red a,n));
symbolic procedure general!-difference!-mod!-p(a,b);
general!-plus!-mod!-p(a,general!-minus!-mod!-p b);
symbolic procedure general!-minus!-mod!-p a;
if null a then nil
else if isdomain a then general!-modular!-minus a
else (lpow a .* general!-minus!-mod!-p lc a) .+
general!-minus!-mod!-p red a;
symbolic procedure general!-reduce!-mod!-p a;
%converts a multivariate poly from normal into modular polynomial;
if null a then nil
else if isdomain a then !*num2f general!-modular!-number a
else adjoin!-term(lpow a,
general!-reduce!-mod!-p lc a,
general!-reduce!-mod!-p red a);
symbolic procedure general!-make!-modular!-symmetric a;
% input is a multivariate MODULAR poly A with nos in the range 0->(p-1).
% This folds it onto the symmetric range (-p/2)->(p/2);
if null a then nil
else if domainp a then
if a>modulus!/2 then !*num2f(a - current!-modulus)
else a
else adjoin!-term(lpow a,
general!-make!-modular!-symmetric lc a,
general!-make!-modular!-symmetric red a);
endmodule;
module degsets; % degree set processing.
% Authors: A. C. Norman and P. M. A. Moore, 1981.
fluid '(!*trallfac
!*trfac
bad!-case
best!-set!-pointer
dpoly
factor!-level
factor!-trace!-list
factored!-lc
irreducible
modular!-info
one!-complete!-deg!-analysis!-done
previous!-degree!-map
split!-list
valid!-image!-sets);
symbolic procedure check!-degree!-sets(n,multivariate!-case);
% MODULAR!-INFO (vector of size N) contains the modular factors now.
begin scalar degree!-sets,w,x!-is!-factor,degs;
w:=split!-list;
for i:=1:n do <<
if multivariate!-case then
x!-is!-factor:=not numberp get!-image!-content
getv(valid!-image!-sets,cdar w);
degs:=for each v in getv(modular!-info,cdar w) collect ldeg v;
degree!-sets:=
(if x!-is!-factor then 1 . degs else degs)
. degree!-sets;
w:=cdr w >>;
check!-degree!-sets!-1 degree!-sets;
best!-set!-pointer:=cdar split!-list;
if multivariate!-case and factored!-lc then <<
while null(w:=get!-f!-numvec
getv(valid!-image!-sets,best!-set!-pointer))
and (split!-list:=cdr split!-list) do
best!-set!-pointer:=cdar split!-list;
if null w then bad!-case:=t >>;
% make sure the set is ok for distributing the
% leading coefft where necessary;
end;
symbolic procedure check!-degree!-sets!-1 l;
% L is a list of degree sets. Try to discover if the entries
% in it are consistent, or if they imply that some of the
% modular splittings were 'false';
begin
scalar i,degree!-map,degree!-map1,dpoly,
plausible!-split!-found,target!-count;
factor!-trace <<
printc "Degree sets are:";
for each s in l do <<
prin2 " ";
for each n in s do <<
prin2 " "; prin2 n >>;
terpri() >> >>;
dpoly:=sum!-list car l;
target!-count:=length car l;
for each s in cdr l do
target!-count:=min(target!-count,length s);
% This used to be IMIN, but since it was the only use, it was
% eliminated.
if null previous!-degree!-map then <<
degree!-map:=mkvect dpoly;
% To begin with all degrees of factors may be possible;
for i:=0:dpoly do putv(degree!-map,i,t) >>
else <<
factor!-trace "Refine an existing degree map";
degree!-map:=previous!-degree!-map >>;
degree!-map1:=mkvect dpoly;
for each s in l do <<
% For each degree set S I will collect in DEGREE-MAP1 a
% bitmap showing what degree factors would be consistent
% with that set. By ANDing together all these maps
% (into DEGREE-MAP) I find what degrees for factors are
% consistent with the whole of the information I have;
for i:=0:dpoly do putv(degree!-map1,i,nil);
putv(degree!-map1,0,t);
putv(degree!-map1,dpoly,t);
for each d in s do for i:=dpoly#-d#-1 step -1 until 0 do
if getv(degree!-map1,i) then
putv(degree!-map1,i#+d,t);
for i:=0:dpoly do
putv(degree!-map,i,getv(degree!-map,i) and
getv(degree!-map1,i)) >>;
factor!-trace <<
printc "Possible degrees for factors are: ";
for i:=1:dpoly#-1 do
if getv(degree!-map,i) then << prin2 i; prin2 " " >>;
terpri() >>;
i:=dpoly#-1;
while i#>0 do if getv(degree!-map,i) then i:=-1
else i:=i#-1;
if i=0 then <<
factor!-trace
printc "Degree analysis proves polynomial irreducible";
return irreducible:=t >>;
for each s in l do if length s=target!-count then begin
% Sets with too many factors are not plausible anyway;
i:=s;
while i and getv(degree!-map,car i) do i:=cdr i;
% If I drop through with I null it was because the set was
% consistent, otherwise it represented a false split;
if null i then plausible!-split!-found:=t end;
previous!-degree!-map:=degree!-map;
if plausible!-split!-found or one!-complete!-deg!-analysis!-done
then return nil;
% PRINTC "Going to try getting some more images";
return bad!-case:=t
end;
symbolic procedure sum!-list l;
if null cdr l then car l
else car l #+ sum!-list cdr l;
endmodule;
module facmod; % Modular factorization: discover the factor count mod p.
% Authors: A. C. Norman and P. M. A. Moore, 1979.
fluid '(!*timings
current!-modulus
dpoly
dwork1
dwork2
known!-factors
linear!-factors
m!-image!-variable
modular!-info
null!-space!-basis
number!-needed
poly!-mod!-p
poly!-vector
safe!-flag
split!-list
work!-vector1
work!-vector2);
safe!-flag:=carcheck 0; % For speed of array access - important here;
symbolic procedure get!-factor!-count!-mod!-p
(n,poly!-mod!-p,p,x!-is!-factor);
% gets the factor count mod p from the nth image using the
% first half of Berlekamp's method;
begin scalar old!-m,f!-count,wtime;
old!-m:=set!-modulus p;
% PRIN2 "prime = ";% printc current!-modulus;
% PRIN2 "degree = ";% printc ldeg poly!-mod!-p;
trace!-time display!-time("Entered GET-FACTOR-COUNT after ",time());
wtime:=time();
f!-count:=modular!-factor!-count();
trace!-time display!-time("Factor count obtained in ",time()-wtime);
split!-list:=
((if x!-is!-factor then car f!-count#+1 else car f!-count) . n)
. split!-list;
putv(modular!-info,n,cdr f!-count);
set!-modulus old!-m
end;
symbolic procedure modular!-factor!-count();
begin
scalar poly!-vector,wvec1,wvec2,x!-to!-p,
n,wtime,w,lin!-f!-count,null!-space!-basis;
known!-factors:=nil;
dpoly:=ldeg poly!-mod!-p;
wvec1:=mkvect (2#*dpoly);
wvec2:=mkvect (2#*dpoly);
x!-to!-p:=mkvect dpoly;
poly!-vector:=mkvect dpoly;
for i:=0:dpoly do putv(poly!-vector,i,0);
poly!-to!-vector poly!-mod!-p;
w:=count!-linear!-factors!-mod!-p(wvec1,wvec2,x!-to!-p);
lin!-f!-count:=car w;
if dpoly#<4 then return
(if dpoly=0 then lin!-f!-count
else lin!-f!-count#+1) .
list(lin!-f!-count . cadr w,
dpoly . poly!-vector,
nil);
% When I use Berlekamp I certainly know that the polynomial
% involved has no linear factors;
wtime:=time();
null!-space!-basis:=use!-berlekamp(x!-to!-p,caddr w,wvec1);
trace!-time display!-time("Berlekamp done in ",time()-wtime);
n:=lin!-f!-count #+ length null!-space!-basis #+ 1;
% there is always 1 more factor than the number of
% null vectors we have picked up;
return n . list(
lin!-f!-count . cadr w,
dpoly . poly!-vector,
null!-space!-basis)
end;
%**********************************************************************;
% Extraction of linear factors is done specially;
symbolic procedure count!-linear!-factors!-mod!-p(wvec1,wvec2,x!-to!-p);
% Compute gcd(x**p-x,u). It will be the product of all the
% linear factors of u mod p;
begin scalar dx!-to!-p,lin!-f!-count,linear!-factors;
for i:=0:dpoly do putv(wvec2,i,getv(poly!-vector,i));
dx!-to!-p:=make!-x!-to!-p(current!-modulus,wvec1,x!-to!-p);
for i:=0:dx!-to!-p do putv(wvec1,i,getv(x!-to!-p,i));
if dx!-to!-p#<1 then <<
if dx!-to!-p#<0 then putv(wvec1,0,0);
putv(wvec1,1,modular!-minus 1);
dx!-to!-p:=1 >>
else <<
putv(wvec1,1,modular!-difference(getv(wvec1,1),1));
if dx!-to!-p=1 and getv(wvec1,1)=0 then
if getv(wvec1,0)=0 then dx!-to!-p:=-1
else dx!-to!-p:=0 >>;
if dx!-to!-p#<0 then
lin!-f!-count:=copy!-vector(wvec2,dpoly,wvec1)
else lin!-f!-count:=gcd!-in!-vector(wvec1,dx!-to!-p,
wvec2,dpoly);
linear!-factors:=mkvect lin!-f!-count;
for i:=0:lin!-f!-count do
putv(linear!-factors,i,getv(wvec1,i));
dpoly:=quotfail!-in!-vector(poly!-vector,dpoly,
linear!-factors,lin!-f!-count);
return list(lin!-f!-count,linear!-factors,dx!-to!-p)
end;
symbolic procedure make!-x!-to!-p(p,wvec1,x!-to!-p);
begin scalar dx!-to!-p,dw1;
if p#<dpoly then <<
for i:=0:p#-1 do putv(x!-to!-p,i,0);
putv(x!-to!-p,p,1);
return p >>;
dx!-to!-p:=make!-x!-to!-p(p/2,wvec1,x!-to!-p);
dw1:=times!-in!-vector(x!-to!-p,dx!-to!-p,x!-to!-p,dx!-to!-p,wvec1);
dw1:=remainder!-in!-vector(wvec1,dw1,
poly!-vector,dpoly);
if not(iremainder(p,2)=0) then <<
for i:=dw1 step -1 until 0 do
putv(wvec1,i#+1,getv(wvec1,i));
putv(wvec1,0,0);
dw1:=remainder!-in!-vector(wvec1,dw1#+1,
poly!-vector,dpoly) >>;
for i:=0:dw1 do putv(x!-to!-p,i,getv(wvec1,i));
return dw1
end;
symbolic procedure find!-linear!-factors!-mod!-p(p,n);
% P is a vector representing a polynomial of degree N which has
% only linear factors. Find all the factors and return a list of
% them;
begin
scalar root,var,w,vec1;
if n#<1 then return nil;
vec1:=mkvect 1;
putv(vec1,1,1);
root:=0;
while (n#>1) and not (root #> current!-modulus) do <<
w:=evaluate!-in!-vector(p,n,root);
if w=0 then << %a factor has been found!!;
if var=nil then
var:=mksp(m!-image!-variable,1) . 1;
w:=!*f2mod
adjoin!-term(car var,cdr var,!*n2f modular!-minus root);
known!-factors:=w . known!-factors;
putv(vec1,0,modular!-minus root);
n:=quotfail!-in!-vector(p,n,vec1,1) >>;
root:=root#+1 >>;
known!-factors:=
vector!-to!-poly(p,n,m!-image!-variable) . known!-factors
end;
%**********************************************************************;
% Berlekamp's algorithm part 1: find null space basis giving factor
% count;
symbolic procedure use!-berlekamp(x!-to!-p,dx!-to!-p,wvec1);
% Set up a basis for the set of remaining (nonlinear) factors
% using Berlekamp's algorithm;
begin
scalar berl!-m,berl!-m!-size,w,
dcurrent,current!-power,wtime;
berl!-m!-size:=dpoly#-1;
berl!-m:=mkvect berl!-m!-size;
for i:=0:berl!-m!-size do <<
w:=mkvect berl!-m!-size;
for j:=0:berl!-m!-size do putv(w,j,0); %initialize to zero;
putv(berl!-m,i,w) >>;
% Note that column zero of the matrix (as used in the
% standard version of Berlekamp's algorithm) is not in fact
% needed and is not used here;
% I want to set up a matrix that has entries
% x**p, x**(2*p), ... , x**((n-1)*p)
% as its columns,
% where n is the degree of poly-mod-p
% and all the entries are reduced mod poly-mod-p;
% Since I computed x**p I have taken out some linear factors,
% so reduce it further;
dx!-to!-p:=remainder!-in!-vector(x!-to!-p,dx!-to!-p,
poly!-vector,dpoly);
dcurrent:=0;
current!-power:=mkvect berl!-m!-size;
putv(current!-power,0,1);
for i:=1:berl!-m!-size do <<
if current!-modulus#>dpoly then
dcurrent:=times!-in!-vector(
current!-power,dcurrent,
x!-to!-p,dx!-to!-p,
wvec1)
else << % Multiply by shifting;
for i:=0:current!-modulus#-1 do
putv(wvec1,i,0);
for i:=0:dcurrent do
putv(wvec1,current!-modulus#+i,
getv(current!-power,i));
dcurrent:=dcurrent#+current!-modulus >>;
dcurrent:=remainder!-in!-vector(
wvec1,dcurrent,
poly!-vector,dpoly);
for j:=0:dcurrent do
putv(getv(berl!-m,j),i,putv(current!-power,j,
getv(wvec1,j)));
% also I need to subtract 1 from the diagonal of the matrix;
putv(getv(berl!-m,i),i,
modular!-difference(getv(getv(berl!-m,i),i),1)) >>;
wtime:=time();
% print!-m("Q matrix",berl!-m,berl!-m!-size);
w := find!-null!-space(berl!-m,berl!-m!-size);
trace!-time display!-time("Null space found in ",time()-wtime);
return w
end;
symbolic procedure find!-null!-space(berl!-m,berl!-m!-size);
% Diagonalize the matrix to find its rank and hence the number of
% factors the input polynomial had;
begin scalar null!-space!-basis;
% find a basis for the null-space of the matrix;
for i:=1:berl!-m!-size do
null!-space!-basis:=
clear!-column(i,null!-space!-basis,berl!-m,berl!-m!-size);
% print!-m("Null vectored",berl!-m,berl!-m!-size);
return
tidy!-up!-null!-vectors(null!-space!-basis,berl!-m,berl!-m!-size)
end;
symbolic procedure print!-m(m,berl!-m,berl!-m!-size);
<< printc m;
for i:=0:berl!-m!-size do <<
for j:=0:berl!-m!-size do <<
prin2 getv(getv(berl!-m,i),j);
ttab((4#*j)#+4) >>;
terpri() >> >>;
symbolic procedure clear!-column(i,
null!-space!-basis,berl!-m,berl!-m!-size);
% Process column I of the matrix so that (if possible) it
% just has a '1' in row I and zeros elsewhere;
begin
scalar ii,w;
% I want to bring a non-zero pivot to the position (i,i)
% and then add multiples of row i to all other rows to make
% all but the i'th element of column i zero. First look for
% a suitable pivot;
ii:=0;
search!-for!-pivot:
if getv(getv(berl!-m,ii),i)=0 or
((ii#<i) and not(getv(getv(berl!-m,ii),ii)=0)) then
if (ii:=ii#+1)#>berl!-m!-size then
return (i . null!-space!-basis)
else go to search!-for!-pivot;
% Here ii references a row containing a suitable pivot element for
% column i. Permute rows in the matrix so as to bring the pivot onto
% the diagonal;
w:=getv(berl!-m,ii);
putv(berl!-m,ii,getv(berl!-m,i));
putv(berl!-m,i,w);
% swop rows ii and i ;
w:=modular!-minus modular!-reciprocal getv(getv(berl!-m,i),i);
% w = -1/pivot, and is used in zeroing out the rest of column i;
for row:=0:berl!-m!-size do
if row neq i then begin
scalar r; %process one row;
r:=getv(getv(berl!-m,row),i);
if not(r=0) then <<
r:=modular!-times(r,w);
%that is now the multiple of row i that must be added to row ii;
for col:=i:berl!-m!-size do
putv(getv(berl!-m,row),col,
modular!-plus(getv(getv(berl!-m,row),col),
modular!-times(r,getv(getv(berl!-m,i),col)))) >>
end;
for col:=i:berl!-m!-size do
putv(getv(berl!-m,i),col,
modular!-times(getv(getv(berl!-m,i),col),w));
return null!-space!-basis
end;
symbolic procedure tidy!-up!-null!-vectors(null!-space!-basis,
berl!-m,berl!-m!-size);
begin
scalar row!-to!-use;
row!-to!-use:=berl!-m!-size#+1;
null!-space!-basis:=
for each null!-vector in null!-space!-basis collect
build!-null!-vector(null!-vector,
getv(berl!-m,row!-to!-use:=row!-to!-use#-1),berl!-m);
berl!-m:=nil; % Release the store for full matrix;
% prin2 "Null vectors: ";
% print null!-space!-basis;
return null!-space!-basis
end;
symbolic procedure build!-null!-vector(n,vec,berl!-m);
% At the end of the elimination process (the CLEAR-COLUMN loop)
% certain columns, indicated by the entries in NULL-SPACE-BASIS
% will be null vectors, save for the fact that they need a '1'
% inserted on the diagonal of the matrix. This procedure copies
% these null-vectors into some of the vectors that represented
% rows of the Berlekamp matrix;
begin
% putv(vec,0,0); % Not used later!!;
for i:=1:n#-1 do
putv(vec,i,getv(getv(berl!-m,i),n));
putv(vec,n,1);
% for i:=n#+1:berl!-m!-size do
% putv(vec,i,0);
return vec . n
end;
%**********************************************************************;
% Berlekamp's algorithm part 2: retrieving the factors mod p;
symbolic procedure get!-factors!-mod!-p(n,p);
% given the modular info (for the nth image) generated by the
% previous half of Berlekamp's method we can reconstruct the
% actual factors mod p;
begin scalar nth!-modular!-info,old!-m,wtime;
nth!-modular!-info:=getv(modular!-info,n);
old!-m:=set!-modulus p;
wtime:=time();
putv(modular!-info,n,
convert!-null!-vectors!-to!-factors nth!-modular!-info);
trace!-time display!-time("Factors constructed in ",time()-wtime);
set!-modulus old!-m
end;
symbolic procedure convert!-null!-vectors!-to!-factors m!-info;
% Using the null space found, complete the job
% of finding modular factors by taking gcd's of the
% modular input polynomial and variants on the
% null space generators;
begin
scalar number!-needed,factors,
work!-vector1,dwork1,work!-vector2,dwork2,wtime;
known!-factors:=nil;
wtime:=time();
find!-linear!-factors!-mod!-p(cdar m!-info,caar m!-info);
trace!-time display!-time("Linear factors found in ",time()-wtime);
dpoly:=caadr m!-info;
poly!-vector:=cdadr m!-info;
null!-space!-basis:=caddr m!-info;
if dpoly=0 then return known!-factors; % All factors were linear;
if null null!-space!-basis then
return known!-factors:=
vector!-to!-poly(poly!-vector,dpoly,m!-image!-variable) .
known!-factors;
number!-needed:=length null!-space!-basis;
% count showing how many more factors I need to find;
work!-vector1:=mkvect dpoly;
work!-vector2:=mkvect dpoly;
factors:=list (poly!-vector . dpoly);
try!-next!-null:
if null!-space!-basis=nil then
errorf "RAN OUT OF NULL VECTORS TOO EARLY";
wtime:=time();
factors:=try!-all!-constants(factors,
caar null!-space!-basis,cdar null!-space!-basis);
trace!-time display!-time("All constants tried in ",time()-wtime);
if number!-needed=0 then
return known!-factors:=append!-new!-factors(factors,
known!-factors);
null!-space!-basis:=cdr null!-space!-basis;
go to try!-next!-null
end;
symbolic procedure try!-all!-constants(list!-of!-polys,v,dv);
% use gcd's of v, v+1, v+2, ... to try to split up the
% polynomials in the given list;
begin
scalar a,b,aa,s;
% aa is a list of factors that can not be improved using this v,
% b is a list that might be;
aa:=nil; b:=list!-of!-polys;
s:=0;
try!-next!-constant:
putv(v,0,s); % Fix constant term of V to be S;
% wtime:=time();
a:=split!-further(b,v,dv);
% trace!-time display!-time("Polys split further in ",time()-wtime);
b:=cdr a; a:=car a;
aa:=nconc(a,aa);
% Keep aa up to date as a list of polynomials that this poly
% v can not help further with;
if b=nil then return aa; % no more progress possible here;
if number!-needed=0 then return nconc(b,aa);
% no more progress needed;
s:=s#+1;
if s#<current!-modulus then go to try!-next!-constant;
% Here I have run out of choices for the constant
% coefficient in v without splitting everything;
return nconc(b,aa)
end;
symbolic procedure split!-further(list!-of!-polys,v,dv);
% list-of-polys is a list of polynomials. try to split
% its members further by taking gcd's with the polynomial
% v. return (a . b) where the polys in a can not possibly
% be split using v+constant, but the polys in b might
% be;
if null list!-of!-polys then nil . nil
else begin
scalar a,b,gg,q;
a:=split!-further(cdr list!-of!-polys,v,dv);
b:=cdr a; a:=car a;
if number!-needed=0 then go to no!-split;
% if all required factors have been found there is no need to
% search further;
dwork1:=copy!-vector(v,dv,work!-vector1);
dwork2:=copy!-vector(caar list!-of!-polys,cdar list!-of!-polys,
work!-vector2);
dwork1:=gcd!-in!-vector(work!-vector1,dwork1,
work!-vector2,dwork2);
if dwork1=0 or dwork1=cdar list!-of!-polys then go to no!-split;
dwork2:=copy!-vector(caar list!-of!-polys,cdar list!-of!-polys,
work!-vector2);
dwork2:=quotfail!-in!-vector(work!-vector2,dwork2,
work!-vector1,dwork1);
% Here I have a splitting;
gg:=mkvect dwork1;
copy!-vector(work!-vector1,dwork1,gg);
a:=((gg . dwork1) . a);
copy!-vector(work!-vector2,dwork2,q:=mkvect dwork2);
b:=((q . dwork2) . b);
number!-needed:=number!-needed#-1;
return (a . b);
no!-split:
return (a . ((car list!-of!-polys) . b))
end;
symbolic procedure append!-new!-factors(a,b);
% Convert to REDUCE (rather than vector) form;
if null a then b
else
vector!-to!-poly(caar a,cdar a,m!-image!-variable) .
append!-new!-factors(cdr a,b);
carcheck safe!-flag; % Restore status quo;
endmodule;
module factrr; % Full factorization of polynomials.
% Author: P. M. A. Moore, 1979.
fluid '(!*all!-contents
!*exp
!*ezgcd
!*force!-prime
!*gcd
!*kernreverse
!*mcd
!*timings
!*trfac
base!-time
current!-modulus
dmode!*
factor!-count
factor!-level
factor!-trace!-list
gc!-base!-time
last!-displayed!-gc!-time
last!-displayed!-time
m!-image!-variable
modulus!/2
polynomial!-to!-factor
polyzero);
global '(!*ifactor);
symbolic procedure factoreval u;
% Factorize the polynomial in the car of u, returning the factors found.
% If cadr u exists then if it is a number, use it as a force prime.
% Otherwise, use cadr u as a fill object, and check to see if caddr u
% is now a force prime.
begin scalar p,w,!*force!-prime,x,z,factor!-count;
p := length u;
if p<1 or p>3
then rederr "FACTORIZE called with wrong number of arguments";
p := !*q2f simp!* car u;
if cdr u then
<<w := cadr u;
if fixp w then <<!*force!-prime := w; w := nil>>
else if cddr u and fixp caddr u
then !*force!-prime := caddr u;
if !*force!-prime and not primep !*force!-prime
then typerr(!*force!-prime,"prime")>>;
x := if dmode!*
then if z := get(dmode!*,'factorfn) then apply1(z,p)
else rederr
list("Factorization not supported over domain",
get(dmode!*,'dname))
else factorf1(p,!*force!-prime);
% Note that car x is expected to be a number.
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,w);
if numberp z then return z % 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.
return if cadr z = 1 then car z . cddr z
else if !*ifactor and numberp cadr z and fixp cadr z
then car z .
append(pairlist2list reversip zfactor cadr z,
cddr z)
else z
end;
put('factorize,'psopfn,'factoreval);
symbolic procedure pairlist2list u;
for each x in u conc nlist(car x,cdr x);
symbolic procedure factorf u;
% This is the entry to the factorizer that is to be used by programmers
% working at the symbolic level. U is to be a standard form. FACTORF
% hands back a list giving the factors of U. The format of said list is
% described below in the comments with FACTORIZE!-FORM. Entry to the
% factorizer at any level other than this is at the programmers own
% risk!! ;
factorf1(u,nil);
symbolic procedure factorf1(u,!*force!-prime);
% This entry to the factorizer allows one to force
% the code to use some particular prime for its
% modular factorization. It is not for casual
% use;
begin
scalar factor!-level,base!-time,last!-displayed!-time,
gc!-base!-time,last!-displayed!-gc!-time,current!-modulus,
modulus!/2,expsave,!*ezgcd,!*gcd;
if null !*mcd then rederr "Factorization invalid with MCD off";
expsave := !*exp;
!*exp := !*gcd := t; % This code will not work otherwise;
!*ezgcd := t;
if null expsave then u := !*q2f resimp !*f2q u;
set!-time();
factor!-level := 0;
u := factorize!-form u;
!*exp := expsave;
return u
end;
symbolic procedure factorize!-form p;
% input:
% p is a reduce standard form that is to be factorized
% over the integers
% result: (nc . l)
% where nc is numeric (may be just 1)
% and l is list of the form:
% ((p1 . x1) (p2 . x2) .. (pn . xn))
% where p<i> are standard forms and x<i> are integers,
% and p= product<i> p<i>**x<i>;
%
% method:
% (a) reorder polynomial to make the variable of lowest maximum
% degree the main one and the rest ordered similarly;
% (b) use contents and primitive parts to split p up as far as possible
% (c) use square-free decomposition to continue the process
% (c.1) detect & perform special processing on cyclotomic polynomials
% (d) use modular-based method to find factors over integers;
begin scalar new!-korder,old!-korder;
new!-korder:=kernord(p,polyzero);
if !*kernreverse then new!-korder:=reverse new!-korder;
old!-korder:=setkorder new!-korder;
p:=reorder p; % Make var of lowest degree the main one;
p:=factorize!-form1(p,new!-korder);
setkorder old!-korder;
p := (car p . for each w in cdr p collect
(reorder car w . cdr w));
return p
end;
symbolic procedure factorize!-form1(p,given!-korder);
% input:
% p is a reduce standard form that is to be factorized
% over the integers
% given-korder is a list of kernels in the order of importance
% (ie when finding leading terms etc. we use this list)
% See FACTORIZE-FORM above;
if domainp p then (p . nil)
else begin scalar m!-image!-variable,var!-list,
polynomial!-to!-factor,n;
if !*all!-contents then var!-list:=given!-korder
else <<
m!-image!-variable:=car given!-korder;
var!-list:=list m!-image!-variable >>;
return (lambda factor!-level;
<< factor!-trace <<
prin2!* "FACTOR : "; printsf p;
prin2!* "Chosen main variable is ";
printvar m!-image!-variable >>;
polynomial!-to!-factor:=p;
n:=numeric!-content p;
p:=quotf(p,n);
if poly!-minusp p then <<
p:=negf p;
n:=-n >>;
factor!-trace <<
prin2!* "Numeric content = ";
printsf n >>;
p:=factorize!-by!-contents(p,var!-list);
p:=n . sort!-factors p;
factor!-trace <<
terpri(); terpri();
printstr "Final result is:"; fac!-printfactors p >>;
p >>)
(factor!-level+1)
end;
symbolic procedure factorize!-form!-recursion p;
% this is essentially the same as FACTORIZE!-FORM except that
% we must be careful of stray minus signs due to a possible
% reordering in the recursive factoring;
begin scalar s,n,x,res,new!-korder,old!-korder;
new!-korder:=kernord(p,polyzero);
if !*kernreverse then new!-korder:=reverse new!-korder;
old!-korder:=setkorder new!-korder;
p:=reorder p; % Make var of lowest degree the main one;
x:=factorize!-form1(p,new!-korder);
setkorder old!-korder;
n := car x;
x := for each p in cdr x collect (reorder car p . cdr p);
if minusp n then << s:=-1; n:=-n >> else s:=1;
res:=for each ff in x collect
if poly!-minusp car ff then <<
s:=s*((-1)**cdr ff);
(negf car ff . cdr ff) >>
else ff;
if minusp s then errorf list(
"Stray minus sign in recursive factorisation:",x);
return (n . res)
end;
symbolic procedure sort!-factors l;
%sort factors as found into some sort of standard order. The order
%used here is more or less random, but will be self-consistent;
sort(l,function orderfactors);
% ***** Contents and primitive parts as applied to factorization *****
symbolic procedure factorize!-by!-contents(p,v);
%use contents wrt variables in list v to split the
%polynomial p. return a list of factors;
% specification is that on entry p *must* be positive;
if domainp p then
errorf list("FACTORIZE-BY-CONTENTS HANDED DOMAIN ELT:",p)
else if null v then square!.free!.factorize p
else begin scalar c,w,l,wtime;
w:=contents!-with!-respect!-to(p,car v);
% contents!-with!-respect!-to returns a pair (g . c) where
% if g=nil the content is just c, otherwise g is a power
% [ x ** n ] and g*c is the content;
if not null car w then <<
% here a power of v divides p;
l:=(!*k2f caar w . cdar w) . nil;
p:=quotfail(p,!*p2f car w);
if p=1 then return l
else if domainp p then
errorf "P SHOULD NOT BE CONSTANT HERE" >>;
c:=cdr w;
if c=1 then << %no progress here;
if null l then
factor!-trace << prin2!* "Polynomial is primitive wrt ";
prinvar car v; terpri!*(nil) >>
else factor!-trace << printstr "Content is: ";
fac!-printfactors(1 . l) >>;
return if !*all!-contents then
append(factorize!-by!-contents(p,cdr v),l)
else append(square!.free!.factorize p,l) >>;
p:=quotfail(p,c); %primitive part;
% p is now primitive, so if it is not a real polynomial it
% must be a unit. since input was +ve it had better be +1 !! ;
if p=-1 then
errorf "NEGATIVE PRIMITIVE PART IN FACTORIZE-BY-CONTENTS";
trace!-time printc "Factoring the content:";
wtime:=time();
l:=append(cdr1 factorize!-form!-recursion c,l);
trace!-time display!-time("Content factored in ",
time()-wtime);
factor!-trace <<
prin2!* "Content wrt "; prinvar car v; prin2!* " is: ";
printsf comfac!-to!-poly w;
printstr "Factors of content are: ";
fac!-printfactors(1 . l) >>;
if p=1 then return l
else if !*all!-contents then
return append(factorize!-by!-contents(p,cdr v),l)
else return append(square!.free!.factorize p,l)
end;
symbolic procedure cdr1 a;
if car a=1 then cdr a
else errorf list("NUMERIC CONTENT NOT EXTRACTED:",car a);
endmodule;
module facuni;
% Authors: A. C. Norman and P. M. A. Moore, 1979;
fluid '(!*force!-prime
!*trfac
alphalist
bad!-case
best!-factor!-count
best!-known!-factors
best!-modulus
best!-set!-pointer
chosen!-prime
factor!-level
factor!-trace!-list
forbidden!-primes
hensel!-growth!-size
input!-leading!-coefficient
input!-polynomial
irreducible
known!-factors
m!-image!-variable
modular!-info
no!-of!-best!-primes
no!-of!-random!-primes
non!-monic
null!-space!-basis
number!-of!-factors
one!-complete!-deg!-analysis!-done
poly!-mod!-p
previous!-degree!-map
reduction!-count
split!-list
target!-factor!-count
univariate!-factors
univariate!-input!-poly
valid!-primes);
symbolic procedure univariate!-factorize poly;
% input poly a primitive square-free univariate polynomial at least
% quadratic and with +ve lc. output is a list of the factors of poly
% over the integers ;
if testx!*!*n!+1 poly then
factorizex!*!*n!+1(m!-image!-variable,ldeg poly,1)
else if testx!*!*n!-1 poly then
factorizex!*!*n!-1(m!-image!-variable,ldeg poly,1)
else univariate!-factorize1 poly;
symbolic procedure univariate!-factorize1 poly;
begin scalar
valid!-primes,univariate!-input!-poly,best!-set!-pointer,
number!-of!-factors,irreducible,forbidden!-primes,
no!-of!-best!-primes,no!-of!-random!-primes,bad!-case,
target!-factor!-count,modular!-info,univariate!-factors,
hensel!-growth!-size,alphalist,previous!-degree!-map,
one!-complete!-deg!-analysis!-done,reduction!-count,
multivariate!-input!-poly;
%note that this code works by using a local database of
%fluid variables that are updated by the subroutines directly
%called here. this allows for the relativly complicated
%interaction between flow of data and control that occurs in
%the factorization algorithm;
factor!-trace <<
prin2!* "Univariate polynomial="; printsf poly;
printstr
"The polynomial is univariate, primitive and square-free";
printstr "so we can treat it slightly more specifically. We";
printstr "factorise mod several primes,then pick the best one";
printstr "to use in the Hensel construction." >>;
initialize!-univariate!-fluids poly;
% set up the fluids to start things off;
tryagain:
get!-some!-random!-primes();
choose!-the!-best!-prime();
if irreducible then <<
univariate!-factors:=list univariate!-input!-poly;
goto exit >>
else if bad!-case then <<
bad!-case:=nil; goto tryagain >>;
reconstruct!-factors!-over!-integers();
if irreducible then <<
univariate!-factors:=list univariate!-input!-poly;
goto exit >>;
exit:
factor!-trace <<
printstr "The univariate factors are:";
for each ff in univariate!-factors do printsf ff >>;
return univariate!-factors
end;
%**********************************************************************
% univariate factorization part 1. initialization and setting fluids;
symbolic procedure initialize!-univariate!-fluids u;
% Set up the fluids to be used in factoring primitive poly;
begin
if !*force!-prime then <<
no!-of!-random!-primes:=1;
no!-of!-best!-primes:=1 >>
else <<
no!-of!-random!-primes:=5;
% we generate this many modular images and calculate
% their factor counts;
no!-of!-best!-primes:=3;
% we find the modular factors of this many;
>>;
univariate!-input!-poly:=u;
target!-factor!-count:=ldeg u
end;
%**********************************************************************;
% univariate factorization part 2. creating modular images and picking
% the best one;
symbolic procedure get!-some!-random!-primes();
% here we create a number of random primes to reduce the input mod p;
begin scalar chosen!-prime,poly!-mod!-p,i;
valid!-primes:=mkvect no!-of!-random!-primes;
i:=0;
while i < no!-of!-random!-primes do <<
poly!-mod!-p:=
find!-a!-valid!-prime(lc univariate!-input!-poly,
univariate!-input!-poly,nil);
if not(poly!-mod!-p='not!-square!-free) then <<
i:=iadd1 i;
putv(valid!-primes,i,chosen!-prime . poly!-mod!-p);
forbidden!-primes:=chosen!-prime . forbidden!-primes
>>
>>
end;
symbolic procedure choose!-the!-best!-prime();
% given several random primes we now choose the best by factoring
% the poly mod its chosen prime and taking one with the
% lowest factor count as the best for hensel growth;
begin scalar split!-list,poly!-mod!-p,null!-space!-basis,
known!-factors,w,n;
modular!-info:=mkvect no!-of!-random!-primes;
for i:=1:no!-of!-random!-primes do <<
w:=getv(valid!-primes,i);
get!-factor!-count!-mod!-p(i,cdr w,car w,nil) >>;
split!-list:=sort(split!-list,function lessppair);
% this now contains a list of pairs (m . n) where
% m is the no: of factors in set no: n. the list
% is sorted with best split (smallest m) first;
if caar split!-list = 1 then <<
irreducible:=t; return nil >>;
w:=split!-list;
for i:=1:no!-of!-best!-primes do <<
n:=cdar w;
get!-factors!-mod!-p(n,car getv(valid!-primes,n));
w:=cdr w >>;
% pick the best few of these and find out their
% factors mod p;
split!-list:=delete(w,split!-list);
% throw away the other sets;
check!-degree!-sets(no!-of!-best!-primes,nil);
% the best set is pointed at by best!-set!-pointer;
one!-complete!-deg!-analysis!-done:=t;
factor!-trace <<
w:=getv(valid!-primes,best!-set!-pointer);
prin2!* "The chosen prime is "; printstr car w;
prin2!* "The polynomial mod "; prin2!* car w;
printstr ", made monic, is:";
printsf cdr w;
printstr "and the factors of this modular polynomial are:";
for each x in getv(modular!-info,best!-set!-pointer)
do printsf x;
>>
end;
%**********************************************************************;
% univariate factorization part 3. reconstruction of the
% chosen image over the integers;
symbolic procedure reconstruct!-factors!-over!-integers();
% the hensel construction from modular case to univariate
% over the integers;
begin scalar best!-modulus,best!-factor!-count,input!-polynomial,
input!-leading!-coefficient,best!-known!-factors,s;
s:=getv(valid!-primes,best!-set!-pointer);
best!-known!-factors:=getv(modular!-info,best!-set!-pointer);
input!-leading!-coefficient:=lc univariate!-input!-poly;
best!-modulus:=car s;
best!-factor!-count:=length best!-known!-factors;
input!-polynomial:=univariate!-input!-poly;
univariate!-factors:=reconstruct!.over!.integers();
if irreducible then return t;
number!-of!-factors:=length univariate!-factors;
if number!-of!-factors=1 then return irreducible:=t
end;
symbolic procedure reconstruct!.over!.integers();
begin scalar w,lclist,non!-monic;
set!-modulus best!-modulus;
for i:=1:best!-factor!-count do
lclist:=input!-leading!-coefficient . lclist;
if not (input!-leading!-coefficient=1) then <<
best!-known!-factors:=
for each ff in best!-known!-factors collect
multf(input!-leading!-coefficient,!*mod2f ff);
non!-monic:=t;
factor!-trace <<
printstr
"(a) Now the polynomial is not monic so we multiply each";
printstr
"of the modular factors, f(i), by the absolute value of";
prin2!* "the leading coefficient: ";
prin2!* input!-leading!-coefficient; printstr '!.;
printstr "To bring the polynomial into agreement with this, we";
prin2!* "multiply it by ";
if best!-factor!-count > 2 then
<< prin2!* input!-leading!-coefficient; prin2!* "**";
printstr isub1 best!-factor!-count >>
else printstr input!-leading!-coefficient >> >>;
w:=uhensel!.extend(input!-polynomial,
best!-known!-factors,lclist,best!-modulus);
if irreducible then return t;
if car w ='ok then return cdr w
else errorf w
end;
% Now some special treatment for cyclotomic polynomials;
symbolic procedure testx!*!*n!+1 u;
not domainp u and (
lc u=1 and
red u = 1);
symbolic procedure testx!*!*n!-1 u;
not domainp u and (
lc u=1 and
red u = -1);
symbolic procedure factorizex!*!*n!+1(var,degree,vorder);
% Deliver factors of (VAR**VORDER)**DEGREE+1 given that it is
% appropriate to treat VAR**VORDER as a kernel;
if evenp degree then factorizex!*!*n!+1(var,degree/2,2*vorder)
else begin
scalar w;
w := factorizex!*!*n!-1(var,degree,vorder);
w := negf car w . cdr w;
return for each p in w collect negate!-variable(var,2*vorder,p)
end;
symbolic procedure negate!-variable(var,vorder,p);
% VAR**(VORDER/2) -> -VAR**(VORDER/2) in the polynomial P;
if domainp p then p
else if mvar p=var then
if remainder(ldeg p,vorder)=0 then
lt p .+ negate!-variable(var,vorder,red p)
else (lpow p .* negf lc p) .+ negate!-variable(var,vorder,red p)
else (lpow p .* negate!-variable(var,vorder,lc p)) .+
negate!-variable(var,vorder,red p);
symbolic procedure integer!-factors n;
% Return integer factors of N, with attached multiplicities. Assumes
% that N is fairly small;
begin
scalar l,q,m,w;
% L is list of results generated so far, Q is current test divisor,
% and M is associated multiplicity;
if n=1 then return '((1 . 1));
q := 2; m := 0;
% Test divide by 2,3,5,7,9,11,13,...
top:
w := divide(n,q);
while cdr w=0 do << n := car w; w := divide(n,q); m := m+1 >>;
if not m=0 then l := (q . m) . l;
if q>car w then <<
if not n=1 then l := (n . 1) . l;
return reversewoc l >>;
% q := ilogor(1,iadd1 q);
q := iadd1 q;
if q #> 3 then q := iadd1 q;
m := 0;
go to top
end;
symbolic procedure factored!-divisors fl;
% FL is an association list of primes and exponents. Return a list
% of all subsets of this list, i.e. of numbers dividing the
% original integer. Exclude '1' from the list;
if null fl then nil
else begin
scalar l,w;
w := factored!-divisors cdr fl;
l := w;
for i := 1:cdar fl do <<
l := list (caar fl . i) . l;
for each p in w do
l := ((caar fl . i) . p) . l >>;
return l
end;
symbolic procedure factorizex!*!*n!-1(var,degree,vorder);
if evenp degree then append(factorizex!*!*n!+1(var,degree/2,vorder),
factorizex!*!*n!-1(var,degree/2,vorder))
else if degree=1 then list((mksp(var,vorder) .* 1) .+ (-1))
else begin
scalar facdeg;
facdeg := '((1 . 1)) . factored!-divisors integer!-factors degree;
return for each fl in facdeg
collect cyclotomic!-polynomial(var,fl,vorder)
end;
symbolic procedure cyclotomic!-polynomial(var,fl,vorder);
% Create Psi<degree>(var**order)
% where degree is given by the association list of primes and
% multiplicities FL;
if not cdar fl=1 then
cyclotomic!-polynomial(var,(caar fl . sub1 cdar fl) . cdr fl,
vorder*caar fl)
else if cdr fl=nil then
if caar fl=1 then (mksp(var,vorder) .* 1) .+ (-1)
else quotfail((mksp(var,vorder*caar fl) .* 1) .+ (-1),
(mksp(var,vorder) .* 1) .+ (-1))
else quotfail(cyclotomic!-polynomial(var,cdr fl,vorder*caar fl),
cyclotomic!-polynomial(var,cdr fl,vorder));
endmodule;
module imageset;
% Authors: A. C. Norman and P. M. A. Moore, 1979;
fluid '(!*force!-prime
!*force!-zero!-set
!*timings
!*trfac
bad!-case
chosen!-prime
current!-modulus
f!-numvec
factor!-level
factor!-trace!-list
factor!-x
factored!-lc
forbidden!-primes
forbidden!-sets
image!-content
image!-lc
image!-mod!-p
image!-poly
image!-set
image!-set!-modulus
kord!*
m!-image!-variable
modulus!/2
multivariate!-input!-poly
no!-of!-primes!-to!-try
othervars
polyzero
save!-zset
usable!-set!-found
vars!-to!-kill
zero!-set!-tried
zerovarset
zset);
%*******************************************************************;
%
% this section deals with the image sets used in
% factorising multivariate polynomials according
% to wang's theories.
% ref: math. comp. vol.32 no.144 oct 1978 pp 1217-1220
% 'an improved multivariate polynomial factoring algorithm'
%
%*******************************************************************;
%*******************************************************************;
% first we have routines for generating the sets
%*******************************************************************;
symbolic procedure generate!-an!-image!-set!-with!-prime
good!-set!-needed;
% given a multivariate poly (in a fluid) we generate an image set
% to make it univariate and also a random prime to use in the
% modular factorization. these numbers are random except that
% we will not allow anything in forbidden!-sets or forbidden!-primes;
begin scalar currently!-forbidden!-sets,u,wtime;
u:=multivariate!-input!-poly;
% a bit of a handful to type otherwise!!!! ;
image!-set:=nil;
currently!-forbidden!-sets:=forbidden!-sets;
tryanotherset:
if image!-set then
currently!-forbidden!-sets:=image!-set .
currently!-forbidden!-sets;
wtime:=time();
image!-set:=get!-new!-set currently!-forbidden!-sets;
% princ "Trying imageset= ";
% printc image!-set;
trace!-time <<
display!-time(" New image set found in ",time()-wtime);
wtime:=time() >>;
image!-lc:=make!-image!-lc!-list(lc u,image!-set);
% list of image lc's wrt different variables in IMAGE-SET;
% princ "Image set to try is:";% printc image!-set;
% prin2!* "L.C. of poly is:";% printsf lc u;
% printc "Image l.c.s with variables substituted on order:";
% for each imlc in image!-lc do printsf imlc;
trace!-time
display!-time(" Image of lc made in ",time()-wtime);
if (caar image!-lc)=0 then goto tryanotherset;
wtime:=time();
image!-poly:=make!-image(u,image!-set);
trace!-time <<
display!-time(" Image poly made in ",time()-wtime);
wtime:=time() >>;
image!-content:=get!.content image!-poly;
% note: the content contains the image variable if it
% is a factor of the image poly;
trace!-time
display!-time(" Content found in ",time()-wtime);
image!-poly:=quotfail(image!-poly,image!-content);
% make sure the image polynomial is primitive which includes
% making the leading coefft positive (-ve content if
% necessary).
% If the image polynomial was of the form k*v^2 where v is
% the image variable then GET.CONTENT will have taken out
% one v and the k leaving the polynomial v here.
% Divisibility by v here thus indicates that the image was
% not square free, and so we will not be able to find a
% sensible prime to use.
if not didntgo quotf(image!-poly,!*k2f m!-image!-variable) then
go to tryanotherset;
wtime:=time();
image!-mod!-p:=find!-a!-valid!-prime(image!-lc,image!-poly,
not numberp image!-content);
if image!-mod!-p='not!-square!-free then goto tryanotherset;
trace!-time <<
display!-time(" Prime and image mod p found in ",time()-wtime);
wtime:=time() >>;
if factored!-lc then
if f!-numvec:=unique!-f!-nos(factored!-lc,image!-content,
image!-set) then <<
usable!-set!-found:=t;
trace!-time
display!-time(" Nos for lc found in ",time()-wtime) >>
else <<
trace!-time display!-time(" Nos for lc failed in ",
time()-wtime);
if (not usable!-set!-found) and good!-set!-needed then
goto tryanotherset >>
end;
symbolic procedure get!-new!-set forbidden!-s;
% associate each variable in vars-to-kill with a random no. mod
% image-set-modulus. If the boolean tagged with a variable is true then
% a value of 1 or 0 is no good and so rejected, however all other
% variables can take these values so they are tried exhaustively before
% using truly random values. sets in forbidden!-s not allowed;
begin scalar old!.m,alist,n,nextzset,w;
if zero!-set!-tried then <<
if !*force!-zero!-set then
errorf "Zero set tried - possibly it was invalid";
image!-set!-modulus:=iadd1 image!-set!-modulus;
old!.m:=set!-modulus image!-set!-modulus;
alist:=for each v in vars!-to!-kill collect
<< n:=modular!-number next!-random!-number();
if n>modulus!/2 then n:=n-current!-modulus;
if cdr v then <<
while n=0
or n=1
or (n = (isub1 current!-modulus)) do
n:=modular!-number next!-random!-number();
if n>modulus!/2 then n:=n-current!-modulus >>;
car v . n >> >>
else <<
old!.m:=set!-modulus image!-set!-modulus;
nextzset:=car zset;
alist:=for each zv in zerovarset collect <<
w:=zv . car nextzset;
nextzset:=cdr nextzset;
w >>;
if othervars then alist:=
append(alist,for each v in othervars collect <<
n:=modular!-number next!-random!-number();
while n=0
or n=1
or (n = (isub1 current!-modulus)) do
n:=modular!-number next!-random!-number();
if n>modulus!/2 then n:=n-current!-modulus;
v . n >>);
if null(zset:=cdr zset) then
if null save!-zset then zero!-set!-tried:=t
else zset:=make!-next!-zset save!-zset;
alist:=for each v in cdr kord!* collect atsoc(v,alist);
% Puts the variables in alist in the right order;
>>;
set!-modulus old!.m;
return if member(alist,forbidden!-s) then
get!-new!-set forbidden!-s
else alist
end;
%**********************************************************************
% now given an image/univariate polynomial find a suitable random prime;
symbolic procedure find!-a!-valid!-prime(lc!-u,u,factor!-x);
% finds a suitable random prime for reducing a poly mod p.
% u is the image/univariate poly. we are not allowed to use
% any of the primes in forbidden!-primes (fluid).
% lc!-u is either numeric or (in the multivariate case) a list of
% images of the lc;
begin scalar currently!-forbidden!-primes,res,prime!-count,v,w;
if factor!-x then u:=multf(u,v:=!*k2f m!-image!-variable);
chosen!-prime:=nil;
currently!-forbidden!-primes:=forbidden!-primes;
prime!-count:=1;
tryanotherprime:
if chosen!-prime then
currently!-forbidden!-primes:=chosen!-prime .
currently!-forbidden!-primes;
chosen!-prime:=get!-new!-prime currently!-forbidden!-primes;
set!-modulus chosen!-prime;
if not atom lc!-u then <<
w:=lc!-u;
while w and
((domainp caar w and not(modular!-number caar w = 0))
or not (domainp caar w or
modular!-number l!-numeric!-c(caar w,cdar w)=0)) do
w:=cdr w;
if w then goto tryanotherprime >>
else if modular!-number lc!-u=0 then goto tryanotherprime;
res:=monic!-mod!-p reduce!-mod!-p u;
if not square!-free!-mod!-p res then
if multivariate!-input!-poly
and (prime!-count:=prime!-count+1)>no!-of!-primes!-to!-try
then <<no!-of!-primes!-to!-try := no!-of!-primes!-to!-try+1;
res:='not!-square!-free>>
else goto tryanotherprime;
if factor!-x and not(res='not!-square!-free) then
res:=quotfail!-mod!-p(res,!*f2mod v);
return res
end;
symbolic procedure get!-new!-prime forbidden!-p;
% get a small prime that is not in the list forbidden!-p;
% we pick one of the first 10 primes if we can;
if !*force!-prime then !*force!-prime
else begin scalar p,primes!-done;
for each pp in forbidden!-p do
if pp<32 then primes!-done:=pp.primes!-done;
tryagain:
if null(p:=random!-teeny!-prime primes!-done) then <<
p:=random!-small!-prime();
primes!-done:='all >>
else primes!-done:=p . primes!-done;
if member(p,forbidden!-p) then goto tryagain;
return p
end;
%***********************************************************************
% find the numbers associated with each factor of the leading
% coefficient of our multivariate polynomial. this will help
% to distribute the leading coefficient later.;
symbolic procedure unique!-f!-nos(v,cont!.u0,im!.set);
% given an image set (im!.set), this finds the numbers associated with
% each factor in v subject to wang's condition (2) on the image set.
% this is an implementation of his algorithm n. if the condition
% is met the result is a vector containing the images of each factor
% in v, otherwise the result is nil;
begin scalar d,k,q,r,lc!.image!.vec;
% v's integer factor is at the front: ;
k:=length cdr v;
% no. of non-trivial factors of v;
if not numberp cont!.u0 then cont!.u0:=lc cont!.u0;
putv(d:=mkvect k,0,abs(cont!.u0 * car v));
% d will contain the special numbers to be used in the
% loop below;
putv(lc!.image!.vec:=mkvect k,0,abs(cont!.u0 * car v));
% vector for result with 0th entry filled in;
v:=cdr v;
% throw away integer factor of v;
% k is no. of non-trivial factors (say f(i)) in v;
% d will contain the nos. associated with each f(i);
% v is now a list of the f(i) (and their multiplicities);
for i:=1:k do
<< q:=abs make!-image(caar v,im!.set);
putv(lc!.image!.vec,i,q);
v:=cdr v;
for j:=isub1 i step -1 until 0 do
<< r:=getv(d,j);
while not onep r do
<< r:=gcd(r,q); q:=q/r >>;
if onep q then <<lc!.image!.vec:=nil; j := -1>>
% if q=1 here then we have failed the condition so exit;
>>;
if null lc!.image!.vec then i := k+1 else putv(d,i,q);
% else q is the ith number we want;
>>;
return lc!.image!.vec
end;
symbolic procedure get!.content u;
% u is a univariate square free poly. gets the content of u (=integer);
% if lc u is negative then the minus sign is pulled out as well;
% nb. the content includes the variable if it is a factor of u;
begin scalar c;
c:=if poly!-minusp u then -(numeric!-content u)
else numeric!-content u;
if not didntgo quotf(u,!*k2f m!-image!-variable) then
c:=adjoin!-term(mksp(m!-image!-variable,1),c,polyzero);
return c
end;
%********************************************************************;
% finally we have the routines that use the numbers generated
% by unique.f.nos to determine the true leading coeffts in
% the multivariate factorization we are doing and which image
% factors will grow up to have which true leading coefft.
%********************************************************************;
symbolic procedure distribute!.lc(r,im!.factors,s,v);
% v is the factored lc of a poly, say u, whose image factors (r of
% them) are in the vector im.factors. s is a list containing the
% image information including the image set, the image poly etc.
% this uses wang's ideas for distributing the factors in v over
% those in im.factors. result is (delta . vector of the lc's of
% the full factors of u) , where delta is the remaining integer part
% of the lc that we have been unable to distribute. ;
(lambda factor!-level;
begin scalar k,delta,div!.count,q,uf,i,d,max!.mult,f,numvec,
dvec,wvec,dtwid,w;
delta:=get!-image!-content s;
% the content of the u image poly;
dist!.lc!.msg1(delta,im!.factors,r,s,v);
v:=cdr v;
% we are not interested in the numeric factors of v;
k:=length v;
% number of things to distribute;
numvec:=get!-f!-numvec s;
% nos. associated with factors in v;
dvec:=mkvect r;
wvec:=mkvect r;
for j:=1:r do <<
putv(dvec,j,1);
putv(wvec,j,delta*lc getv(im!.factors,j)) >>;
% result lc's will go into dvec which we initialize to 1's;
% wvec is a work vector that we use in the division process
% below;
v:=reverse v;
for j:=k step -1 until 1 do
<< % (for each factor in v, call it f(j) );
f:=caar v;
% f(j) itself;
max!.mult:=cdar v;
% multiplicity of f(j) in v (=lc u);
v:=cdr v;
d:=getv(numvec,j);
% number associated with f(j);
i:=1; % we trial divide d into lc of each image
% factor starting with 1st;
div!.count:=0;
% no. of d's that have been distributed;
factor!-trace <<
prin2!* "f("; prin2!* j; prin2!* ")= "; printsf f;
prin2!* "There are "; prin2!* max!.mult;
printstr " of these in the leading coefficient.";
prin2!* "The absolute value of the image of f("; prin2!* j;
prin2!* ")= "; printstr d >>;
while ilessp(div!.count,max!.mult)
and not igreaterp(i,r) do
<< q:=divide(getv(wvec,i),d);
% first trial division;
factor!-trace <<
prin2!* " Trial divide into ";
prin2!* getv(wvec,i); printstr " :" >>;
while (zerop cdr q) and ilessp(div!.count,max!.mult) do
<< putv(dvec,i,multf(getv(dvec,i),f));
% f(j) belongs in lc of ith factor;
factor!-trace <<
prin2!* " It goes so an f("; prin2!* j;
prin2!* ") belongs in ";
printsf getv(im!.factors,i);
printstr " Try again..." >>;
div!.count:=iadd1 div!.count;
% another d done;
putv(wvec,i,car q);
% save the quotient for next factor to distribute;
q:=divide(car q,d);
% try again;
>>;
i:=iadd1 i;
% as many d's as possible have gone into that
% factor so now try next factor;
factor!-trace
<<printstr " no good so try another factor ..." >>>>;
% at this point the whole of f(j) should have been
% distributed by dividing d the maximum no. of times
% (= max!.mult), otherwise we have an extraneous factor;
if ilessp(div!.count,max!.mult) then
<<bad!-case:=t; div!.count := max!.mult>>
>>;
if bad!-case then return;
dist!.lc!.msg2(dvec,im!.factors,r);
if onep delta then
<< for j:=1:r do <<
w:=lc getv(im!.factors,j) /
evaluate!-in!-order(getv(dvec,j),get!-image!-set s);
if w<0 then begin
scalar oldpoly;
delta:= -delta;
oldpoly:=getv(im!.factors,j);
putv(im!.factors,j,negf oldpoly);
% to keep the leading coefficients positive we negate the
% image factors when necessary;
multiply!-alphas(-1,oldpoly,getv(im!.factors,j));
% remember to fix the alphas as well;
end;
putv(dvec,j,multf(abs w,getv(dvec,j))) >>;
dist!.lc!.msg3(dvec,im!.factors,r);
return (delta . dvec)
>>;
% if delta=1 then we know the true lc's exactly so put in their
% integer contents and return with result.
% otherwise try spreading delta out over the factors: ;
dist!.lc!.msg4 delta;
for j:=1:r do
<< dtwid:=evaluate!-in!-order(getv(dvec,j),get!-image!-set s);
uf:=getv(im!.factors,j);
d:=gcddd(lc uf,dtwid);
putv(dvec,j,multf(lc uf/d,getv(dvec,j)));
putv(im!.factors,j,multf(dtwid/d,uf));
% have to fiddle the image factors by an integer multiple;
multiply!-alphas!-recip(dtwid/d,uf,getv(im!.factors,j));
% fix the alphas;
delta:=delta/(dtwid/d)
>>;
% now we've done all we can to distribute delta so we return with
% what's left: ;
if delta<=0 then
errorf list("FINAL DELTA IS -VE IN DISTRIBUTE!.LC",delta);
factor!-trace <<
printstr " Finally we have:";
for j:=1:r do <<
prinsf getv(im!.factors,j);
prin2!* " with l.c. ";
printsf getv(dvec,j) >> >>;
return (delta . dvec)
end) (factor!-level * 10);
symbolic procedure dist!.lc!.msg1(delta,im!.factors,r,s,v);
factor!-trace <<
terpri(); terpri();
printstr "We have a polynomial whose image factors (call";
printstr "them the IM-factors) are:";
prin2!* delta; printstr " (= numeric content, delta)";
printvec(" f(",r,")= ",im!.factors);
prin2!* " wrt the image set: ";
for each x in get!-image!-set s do <<
prinvar car x; prin2!* "="; prin2!* cdr x; prin2!* ";" >>;
terpri!*(nil);
printstr "We also have its true multivariate leading";
printstr "coefficient whose factors (call these the";
printstr "LC-factors) are:";
fac!-printfactors v;
printstr "We want to determine how these LC-factors are";
printstr "distributed over the leading coefficients of each";
printstr "IM-factor. This enables us to feed the resulting";
printstr "image factors into a multivariate Hensel";
printstr "construction.";
printstr "We distribute each LC-factor in turn by dividing";
printstr "its image into delta times the leading coefficient";
printstr "of each IM-factor until it finds one that it";
printstr "divides exactly. The image set is chosen such that";
printstr "this will only happen for the IM-factors to which";
printstr "this LC-factor belongs - (there may be more than";
printstr "one if the LC-factor occurs several times in the";
printstr "leading coefficient of the original polynomial).";
printstr "This choice also requires that we distribute the";
printstr "LC-factors in a specific order:"
>>;
symbolic procedure dist!.lc!.msg2(dvec,im!.factors,r);
factor!-trace <<
printstr "The leading coefficients are now correct to within an";
printstr "integer factor and are as follows:";
for j:=1:r do <<
prinsf getv(im!.factors,j);
prin2!* " with l.c. ";
printsf getv(dvec,j) >> >>;
symbolic procedure dist!.lc!.msg3(dvec,im!.factors,r);
factor!-trace <<
printstr "Since delta=1, we have no non-trivial content of the";
printstr
"image to deal with so we know the true leading coefficients";
printstr
"exactly. We fix the signs of the IM-factors to match those";
printstr "of their true leading coefficients:";
for j:=1:r do <<
prinsf getv(im!.factors,j);
prin2!* " with l.c. ";
printsf getv(dvec,j) >> >>;
symbolic procedure dist!.lc!.msg4 delta;
factor!-trace <<
prin2!* " Here delta is not 1 meaning that we have a content, ";
printstr delta;
printstr "of the image to distribute among the factors somehow.";
printstr "For each IM-factor we can divide its leading";
printstr "coefficient by the image of its determined leading";
printstr "coefficient and see if there is a non-trivial result.";
printstr "This will indicate a factor of delta belonging to this";
printstr "IM-factor's leading coefficient." >>;
endmodule;
module pfactor; % Factorization of polynomials modulo p.
% Author: A. C. Norman, 1978.
fluid '(!*backtrace
!*gcd
base!-time
current!-modulus
gc!-base!-time
last!-displayed!-gc!-time
last!-displayed!-time
m!-image!-variable
modular!-info
modulus!/2
user!-prime);
symbolic procedure pfactor(q,p);
% Q is a standard form. Factorize and return the factors mod p.
begin scalar base!-time,last!-displayed!-time,
gc!-base!-time,last!-displayed!-gc!-time,
user!-prime,current!-modulus,modulus!/2,r,x;
set!-time();
if not numberp p then typerr(p,"number")
else if not primep p then typerr(p,"prime");
user!-prime:=p;
set!-modulus p;
if domainp q or null reduce!-mod!-p lc q then
printc "*** Degenerate case in modular factorization";
if not (length variables!-in!-form q=1) then
rederr "Multivariate input to modular factorization";
r:=reduce!-mod!-p q;
% LNCOEFF := LC R;
x := lnc r;
r :=monic!-mod!-p r;
print!-time "About to call FACTOR-FORM-MOD-P";
r:=errorset(list('factor!-form!-mod!-p,mkquote r),t,!*backtrace);
print!-time "FACTOR-FORM-MOD-P returned";
if not errorp r then return x . car r;
printc "****** FACTORIZATION FAILED******";
return list(1,prepf q) % 1 needed by factorize.
end;
symbolic procedure factor!-form!-mod!-p p;
% input:
% p is a reduce standard form that is to be factorized
% mod prime;
% result:
% ((p1 . x1) (p2 . x2) .. (pn . xn))
% where p<i> are standard forms and x<i> are integers,
% and p= product<i> p<i>**x<i>;
sort!-factors factorize!-by!-square!-free!-mod!-p p;
symbolic procedure factorize!-by!-square!-free!-mod!-p p;
if p=1 then nil
else if domainp p then (p . 1) . nil
else
begin
scalar dp,v;
v:=(mksp(mvar p,1).* 1) .+ nil;
dp:=0;
while evaluate!-mod!-p(p,mvar v,0)=0 do <<
p:=quotfail!-mod!-p(p,v);
dp:=dp+1 >>;
if dp>0 then return ((v . dp) .
factorize!-by!-square!-free!-mod!-p p);
dp:=derivative!-mod!-p p;
if dp=nil then <<
%here p is a something to the power current!-modulus;
p:=divide!-exponents!-by!-p(p,current!-modulus);
p:=factorize!-by!-square!-free!-mod!-p p;
return multiply!-multiplicities(p,current!-modulus) >>;
dp:=gcd!-mod!-p(p,dp);
if dp=1 then return factorize!-pp!-mod!-p p;
%now p is not square-free;
p:=quotfail!-mod!-p(p,dp);
%factorize p and dp separately;
p:=factorize!-pp!-mod!-p p;
dp:=factorize!-by!-square!-free!-mod!-p dp;
% i feel that this scheme is slightly clumsy, but
% square-free decomposition mod p is not as straightforward
% as square free decomposition over the integers, and pfactor
% is probably not going to be slowed down too badly by
% this;
return mergefactors(p,dp)
end;
%**********************************************************************;
% code to factorize primitive square-free polynomials mod p;
symbolic procedure divide!-exponents!-by!-p(p,n);
if isdomain p then p
else (mksp(mvar p,exactquotient(ldeg p,n)) .* lc p) .+
divide!-exponents!-by!-p(red p,n);
symbolic procedure exactquotient(a,b);
begin
scalar w;
w:=divide(a,b);
if cdr w=0 then return car w;
error(50,list("Inexact division",list(a,b,w)))
end;
symbolic procedure multiply!-multiplicities(l,n);
if null l then nil
else (caar l . (n*cdar l)) .
multiply!-multiplicities(cdr l,n);
symbolic procedure mergefactors(a,b);
% a and b are lists of factors (with multiplicities),
% merge them so that no factor occurs more than once in
% the result;
if null a then b
else mergefactors(cdr a,addfactor(car a,b));
symbolic procedure addfactor(a,b);
%add factor a into list b;
if null b then list a
else if car a=caar b then
(car a . (cdr a + cdar b)) . cdr b
else car b . addfactor(a,cdr b);
symbolic procedure factorize!-pp!-mod!-p p;
%input a primitive square-free polynomial p,
% output a list of irreducible factors of p;
begin
scalar vars;
if p=1 then return nil
else if isdomain p then return (p . 1) . nil;
% now I am certain that p is not degenerate;
print!-time "primitive square-free case detected";
vars:=variables!-in!-form p;
if length vars=1 then return unifac!-mod!-p p;
errorf "SHAMBLED IN PFACTOR - MULTIVARIATE CASE RESURFACED"
end;
symbolic procedure unifac!-mod!-p p;
%input p a primitive square-free univariate polynomial
%output a list of the factors of p over z mod p;
begin
scalar modular!-info,m!-image!-variable;
if isdomain p then return nil
else if ldeg p=1 then return (p . 1) . nil;
modular!-info:=mkvect 1;
m!-image!-variable:=mvar p;
get!-factor!-count!-mod!-p(1,p,user!-prime,nil);
print!-time "Factor counts obtained";
get!-factors!-mod!-p(1,user!-prime);
print!-time "Actual factors extracted";
return for each z in getv(modular!-info,1) collect (z . 1)
end;
endmodule;
module vecpoly;
% Authors: A. C. Norman and P. M. A. Moore, 1979;
fluid '(current!-modulus safe!-flag);
%**********************************************************************;
% Routines for working with modular univariate polynomials
% stored as vectors. Used to avoid unwarranted storage management
% in the mod-p factorization process;
safe!-flag:=carcheck 0;
symbolic procedure copy!-vector(a,da,b);
% Copy A into B;
<< for i:=0:da do
putv(b,i,getv(a,i));
da >>;
symbolic procedure times!-in!-vector(a,da,b,db,c);
% Put the product of A and B into C and return its degree.
% C must not overlap with either A or B;
begin
scalar dc,ic,w;
if da#<0 or db#<0 then return minus!-one;
dc:=da#+db;
for i:=0:dc do putv(c,i,0);
for ia:=0:da do <<
w:=getv(a,ia);
for ib:=0:db do <<
ic:=ia#+ib;
putv(c,ic,modular!-plus(getv(c,ic),
modular!-times(w,getv(b,ib)))) >> >>;
return dc
end;
symbolic procedure quotfail!-in!-vector(a,da,b,db);
% Overwrite A with (A/B) and return degree of result.
% The quotient must be exact;
if da#<0 then da
else if db#<0 then errorf "Attempt to divide by zero"
else if da#<db then errorf "Bad degrees in QUOTFAIL-IN-VECTOR"
else begin
scalar dc;
dc:=da#-db; % Degree of result;
for i:=dc step -1 until 0 do begin
scalar q;
q:=modular!-quotient(getv(a,db#+i),getv(b,db));
for j:=0:db#-1 do
putv(a,i#+j,modular!-difference(getv(a,i#+j),
modular!-times(q,getv(b,j))));
putv(a,db#+i,q)
end;
for i:=0:db#-1 do if getv(a,i) neq 0 then
errorf "Quotient not exact in QUOTFAIL!-IN!-VECTOR";
for i:=0:dc do
putv(a,i,getv(a,db#+i));
return dc
end;
symbolic procedure remainder!-in!-vector(a,da,b,db);
% Overwrite the vector A with the remainder when A is
% divided by B, and return the degree of the result;
begin
scalar delta,db!-1,recip!-lc!-b,w;
if db=0 then return minus!-one
else if db=minus!-one then errorf "ATTEMPT TO DIVIDE BY ZERO";
recip!-lc!-b:=modular!-minus modular!-reciprocal getv(b,db);
db!-1:=db#-1; % Leading coeff of B treated specially, hence this;
while not((delta:=da#-db) #< 0) do <<
w:=modular!-times(recip!-lc!-b,getv(a,da));
for i:=0:db!-1 do
putv(a,i#+delta,modular!-plus(getv(a,i#+delta),
modular!-times(getv(b,i),w)));
da:=da#-1;
while not(da#<0) and getv(a,da)=0 do da:=da#-1 >>;
return da
end;
symbolic procedure evaluate!-in!-vector(a,da,n);
% Evaluate A at N;
begin
scalar r;
r:=getv(a,da);
for i:=da#-1 step -1 until 0 do
r:=modular!-plus(getv(a,i),
modular!-times(r,n));
return r
end;
symbolic procedure gcd!-in!-vector(a,da,b,db);
% Overwrite A with the gcd of A and B. On input A and B are
% vectors of coefficients, representing polynomials
% of degrees DA and DB. Return DG, the degree of the gcd;
begin
scalar w;
if da=0 or db=0 then << putv(a,0,1); return 0 >>
else if da#<0 or db#<0 then errorf "GCD WITH ZERO NOT ALLOWED";
top:
% Reduce the degree of A;
da:=remainder!-in!-vector(a,da,b,db);
if da=0 then << putv(a,0,1); return 0 >>
else if da=minus!-one then <<
w:=modular!-reciprocal getv(b,db);
for i:=0:db do putv(a,i,modular!-times(getv(b,i),w));
return db >>;
% Now reduce degree of B;
db:=remainder!-in!-vector(b,db,a,da);
if db=0 then << putv(a,0,1); return 0 >>
else if db=minus!-one then <<
w:=modular!-reciprocal getv(a,da);
if not (w=1) then
for i:=0:da do putv(a,i,modular!-times(getv(a,i),w));
return da >>;
go to top
end;
carcheck safe!-flag;
endmodule;
end;