12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586158715881589159015911592159315941595159615971598159916001601160216031604160516061607160816091610161116121613161416151616161716181619162016211622162316241625162616271628162916301631163216331634163516361637163816391640164116421643164416451646164716481649165016511652165316541655165616571658165916601661166216631664166516661667166816691670167116721673167416751676167716781679168016811682168316841685168616871688168916901691169216931694169516961697169816991700170117021703170417051706170717081709171017111712171317141715171617171718171917201721172217231724172517261727172817291730173117321733173417351736173717381739174017411742174317441745174617471748174917501751175217531754175517561757175817591760176117621763176417651766176717681769177017711772177317741775177617771778177917801781178217831784178517861787178817891790179117921793179417951796179717981799180018011802180318041805180618071808180918101811181218131814181518161817181818191820182118221823182418251826182718281829183018311832183318341835183618371838183918401841184218431844184518461847184818491850185118521853185418551856185718581859186018611862186318641865186618671868186918701871187218731874187518761877187818791880188118821883188418851886188718881889189018911892189318941895189618971898189919001901190219031904190519061907190819091910191119121913191419151916191719181919192019211922192319241925192619271928192919301931193219331934193519361937193819391940194119421943194419451946194719481949195019511952195319541955195619571958195919601961196219631964196519661967196819691970197119721973197419751976197719781979198019811982198319841985198619871988198919901991199219931994199519961997199819992000200120022003200420052006200720082009201020112012201320142015201620172018201920202021202220232024202520262027202820292030203120322033203420352036203720382039204020412042204320442045204620472048204920502051205220532054205520562057205820592060206120622063206420652066206720682069207020712072207320742075207620772078 |
- 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;
|