Artifact 99af573d0824375845632aa6f263f0377106d27ab3b63ad45a226bbb0a0aaaa9:
- File
r33/factor.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 74617) [annotate] [blame] [check-ins using] [more...]
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;