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