123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911 |
- module ezgcdf; % Polynomial GCD algorithms.
- % Author: A. C. Norman, 1981.
- fluid '(!*exp
- !*gcd
- !*heugcd
- !*overview
- !*trfac
- alphalist
- bad!-case
- best!-known!-factors
- current!-modulus
- dmode!*
- factor!-level
- factor!-trace!-list
- full!-gcd
- hensel!-growth!-size
- image!-factors
- image!-set
- irreducible
- kord!*
- m!-image!-variable
- multivariate!-factors
- multivariate!-input!-poly
- non!-monic
- no!-of!-primes!-to!-try
- number!-of!-factors
- prime!-base
- reconstructing!-gcd
- reduced!-degree!-lclst
- reduction!-count
- target!-factor!-count
- true!-leading!-coeffts
- unlucky!-case);
- global '(erfg!*);
- symbolic procedure ezgcdf(u,v);
- % Entry point for REDUCE call in GCDF. We have to make sure that
- % the kernel order is correct if an error occurs in ezgcdf1.
- begin scalar erfgx,kordx,x;
- erfgx := erfg!*;
- kordx := kord!*;
- x := errorset2{'ezgcdf1,mkquote u,mkquote v};
- if null errorp x then return first x;
- % If ezgcdf fails, erfg!* can be set to t,
- % (e.g., in invlap(c/(p^3/8-9p^2/4+27/2*p-27)^2,p,t)), and
- % the kernel order not properly reset.
- erfg!* := erfgx;
- setkorder kordx;
- return gcdf1(u,v)
- end;
-
- symbolic procedure ezgcdf1(u,v);
- % Entry point for REDUCE call in GCDF.
- poly!-abs gcdlist list(u,v) where factor!-level=0;
-
- %symbolic procedure simpezgcd u;
- % calculate the gcd of the polynomials given as arguments;
- % begin
- % scalar factor!-level,w;
- % factor!-level:=0;
- % u := for each p in u collect <<
- % w := simp!* p;
- % if (denr w neq 1) then
- % rederr "EZGCD requires polynomial arguments";
- % numr w >>;
- % return (poly!-abs gcdlist u) ./ 1
- % end;
- %put('ezgcd,'simpfn,'simpezgcd);
- symbolic procedure simpnprimitive p;
- % Remove any simple numeric factors from the expression P;
- begin
- scalar np,dp;
- if atom p or not atom cdr p then
- rerror(ezgcd,2,"NPRIMITIVE requires just one argument");
- p := simp!* car p;
- if polyzerop(numr p) then return nil ./ 1;
- np := quotfail(numr p,numeric!-content numr p);
- dp := quotfail(denr p,numeric!-content denr p);
- return (np ./ dp)
- end;
-
- put('nprimitive,'simpfn,'simpnprimitive);
-
-
- symbolic procedure poly!-gcd(u,v);
- % U and V are standard forms.
- % Value is the gcd of U and V.
- begin scalar !*exp,z;
- if polyzerop u then return poly!-abs v
- else if polyzerop v then return poly!-abs u
- else if u=1 or v=1 then return 1;
- !*exp := t;
- % The case of one argument exactly dividing the other is
- % detected specially here because it is perhaps a fairly
- % common circumstance.
- if quotf1(u,v) then z := v
- else if quotf1(v,u) then z := u
- else if !*gcd then z := gcdlist list(u,v)
- else z := 1;
- return poly!-abs z
- end;
-
- % moved('gcdf,'poly!-gcd);
-
- symbolic procedure ezgcd!-comfac p;
- %P is a standard form
- %CAR of result is lowest common power of leading kernel in
- %every term in P (or NIL). CDR is gcd of all coefficients of
- %powers of leading kernel;
- if domainp p then nil . poly!-abs p
- else if null red p then lpow p . poly!-abs lc p
- else begin
- scalar power,coeflist,var;
- % POWER will be the first part of the answer returned,
- % COEFLIST will collect a list of all coefs in the polynomial
- % P viewed as a poly in its main variable,
- % VAR is the main variable concerned;
- var := mvar p;
- while mvar p=var and not domainp red p do <<
- coeflist := lc p . coeflist;
- p:=red p >>;
- if mvar p=var then <<
- coeflist := lc p . coeflist;
- if null red p then power := lpow p
- else coeflist := red p . coeflist >>
- else coeflist := p . coeflist;
- return power . gcdlist coeflist
- end;
-
- symbolic procedure gcd!-with!-number(n,a);
- % n is a number, a is a polynomial - return their gcd, given that
- % n is non-zero;
- if n=1 or not atom n or flagp(dmode!*,'field) then 1
- else if domainp a
- then if a=nil then abs n
- else if not atom a then 1
- else gcddd(n,a)
- else gcd!-with!-number(gcd!-with!-number(n,lc a),red a);
- % moved('gcdfd,'gcd!-with!-number);
-
-
- symbolic procedure contents!-with!-respect!-to(p,v);
- if domainp p then nil . poly!-abs p
- else if mvar p=v then ezgcd!-comfac p
- else begin
- scalar w,y;
- y := updkorder v;
- w := ezgcd!-comfac reorder p;
- setkorder y;
- return w
- end;
-
- symbolic procedure numeric!-content form;
- % Find numeric content of non-zero polynomial.
- if domainp form then absf form
- else if null red form then numeric!-content lc form
- else begin
- scalar g1;
- g1 := numeric!-content lc form;
- if not (g1=1) then g1 := gcddd(g1,numeric!-content red form);
- return g1
- end;
-
- symbolic procedure gcdlist l;
- % Return the GCD of all the polynomials in the list L.
- %
- % First find all variables mentioned in the polynomials in L,
- % and remove monomial content from them all. If in the process
- % a constant poly is found, take special action. If then there
- % is some variable that is mentioned in all the polys in L, and
- % which occurs only linearly in one of them establish that as
- % main variable and proceed to GCDLIST3 (which will take
- % a special case exit). Otherwise, if there are any variables that
- % do not occur in all the polys in L they can not occur in the GCD,
- % so take coefficients with respect to them to get a longer list of
- % smaller polynomials - restart. Finally we have a set of polys
- % all involving exactly the same set of variables;
- if null l then nil
- else if null cdr l then poly!-abs car l
- else if domainp car l then gcdld(cdr l,car l)
- else begin
- scalar l1,gcont,x;
- % Copy L to L1, but on the way detect any domain elements
- % and deal with them specially;
- while not null l do <<
- if null car l then l := cdr l
- else if domainp car l then <<
- l1 := list list gcdld(cdr l,gcdld(mapcarcar l1,car l));
- l := nil >>
- else <<
- l1 := (car l . powers1 car l) . l1;
- l := cdr l >> >>;
- if null l1 then return nil
- else if null cdr l1 then return poly!-abs caar l1;
- % Now L1 is a list where each polynomial is paired with information
- % about the powers of variables in it;
- gcont := nil; % Compute monomial content on things in L;
- x := nil; % First time round flag;
- l := for each p in l1 collect begin
- scalar gcont1,gcont2,w;
- % Set GCONT1 to least power information, and W to power
- % difference;
- w := for each y in cdr p
- collect << gcont1 := (car y . cddr y) . gcont1;
- car y . (cadr y-cddr y) >>;
- % Now get the monomial content as a standard form (in GCONT2);
- gcont2 := numeric!-content car p;
- if null x then << gcont := gcont1; x := gcont2 >>
- else << gcont := vintersection(gcont,gcont1);
- % Accumulate monomial gcd;
- x := gcddd(x,gcont2) >>;
- for each q in gcont1 do if not(cdr q=0) then
- gcont2 := multf(gcont2,!*p2f mksp(car q,cdr q));
- return quotfail1(car p,gcont2,"Term content division failed")
- . w
- end;
- % Here X is the numeric part of the final GCD.
- for each q in gcont do x := multf(x,!*p2f mksp(car q,cdr q));
- % trace!-time <<
- % prin2!* "Term gcd = ";
- % printsf x >>;
- return poly!-abs multf(x,gcdlist1 l)
- end;
-
-
- symbolic procedure gcdlist1 l;
- % Items in L are monomial-primitive, and paired with power information.
- % Find out what variables are common to all polynomials in L and
- % remove all others;
- begin
- scalar unionv,intersectionv,vord,x,l1,reduction!-count;
- unionv := intersectionv := cdar l;
- for each p in cdr l do <<
- unionv := vunion(unionv,cdr p);
- intersectionv := vintersection(intersectionv,cdr p) >>;
- if null intersectionv then return 1;
- for each v in intersectionv do
- unionv := vdelete(v,unionv);
- % Now UNIONV is list of those variables mentioned that
- % are not common to all polynomials;
- intersectionv := sort(intersectionv,function lesspcdr);
- if cdar intersectionv=1 then <<
- % I have found something that is linear in one of its variables;
- vord := mapcarcar append(intersectionv,unionv);
- l1 := setkorder vord;
- % trace!-time <<
- % prin2 "Selecting "; prin2 caar intersectionv;
- % prin2t " as main because some poly is linear in it" >>;
- x := gcdlist3(for each p in l collect reorder car p,nil,vord);
- setkorder l1;
- return reorder x >>
- else if null unionv then return gcdlist2(l,intersectionv);
- % trace!-time <<
- % prin2 "The variables "; prin2 unionv; prin2t " can be removed" >>;
- vord := setkorder mapcarcar append(unionv,intersectionv);
- l1 := nil;
- for each p in l do
- l1:=split!-wrt!-variables(reorder car p,mapcarcar unionv,l1);
- setkorder vord;
- return gcdlist1(for each p in l1 collect
- (reorder p . total!-degree!-in!-powers(p,nil)))
- end;
-
-
- symbolic procedure gcdlist2(l,vars);
- % Here all the variables in VARS are used in every polynomial
- % in L. Select a good variable ordering;
- begin
- scalar x,x1,gg,lmodp,onestep,vord,oldmod,image!-set,gcdpow,
- unlucky!-case;
- % In the univariate case I do not need to think very hard about
- % the selection of a main variable!! ;
- if null cdr vars
- then return if !*heugcd and (x := heu!-gcd!-list(mapcarcar l))
- then x
- else gcdlist3(mapcarcar l,nil,list caar vars);
- oldmod := set!-modulus nil;
- % If some variable appears at most to degree two in some pair of the
- % polynomials then that will do as a main variable. Note that this is
- % not so useful if the two polynomials happen to be duplicates of each
- % other, but still... ;
- vars := mapcarcar sort(vars,function greaterpcdr);
- % Vars is now arranged with the variable that appears to highest
- % degree anywhere in L first, and the rest in descending order;
- l := for each p in l collect car p .
- sort(cdr p,function lesspcdr);
- l := sort(l,function lesspcdadr);
- % Each list of degree information in L is sorted with lowest degree
- % vars first, and the polynomial with the lowest degree variable
- % of all will come first;
- x := intersection(deg2vars(cdar l),deg2vars(cdadr l));
- if not null x then <<
- % trace!-time << prin2 "Two inputs are at worst quadratic in ";
- % prin2t car x >>;
- go to x!-to!-top >>; % Here I have found two polys with a common
- % variable that they are quadratic in;
- % Now generate modular images of the gcd to guess its degree wrt
- % all possible variables;
-
- % If either (a) modular gcd=1 or (b) modular gcd can be computed with
- % just 1 reduction step, use that information to choose a main variable;
- try!-again: % Modular images may be degenerate;
- set!-modulus random!-prime();
- unlucky!-case := nil;
- image!-set := for each v in vars
- collect (v . modular!-number next!-random!-number());
- % trace!-time <<
- % prin2 "Select variable ordering using P=";
- % prin2 current!-modulus;
- % prin2 " and substitutions from ";
- % prin2t image!-set >>;
- x1 := vars;
- try!-vars:
- if null x1 then go to images!-tried;
- lmodp := for each p in l collect make!-image!-mod!-p(car p,car x1);
- if unlucky!-case then go to try!-again;
- lmodp := sort(lmodp,function lesspdeg);
- gg := gcdlist!-mod!-p(car lmodp,cdr lmodp);
- if domainp gg or (reduction!-count<2 and (onestep:=t)) then <<
- % trace!-time << prin2 "Select "; prin2t car x1 >>;
- x := list car x1; go to x!-to!-top >>;
- gcdpow := (car x1 . ldeg gg) . gcdpow;
- x1 := cdr x1;
- go to try!-vars;
- images!-tried:
- % In default of anything better to do, use image variable such that
- % degree of gcd wrt it is as large as possible;
- vord := mapcarcar sort(gcdpow,function greaterpcdr);
- % trace!-time << prin2 "Select order by degrees: ";
- % prin2t gcdpow >>;
- go to order!-chosen;
-
- x!-to!-top:
- for each v in x do vars := delete(v,vars);
- vord := append(x,vars);
- order!-chosen:
- % trace!-time << prin2 "Selected Var order = "; prin2t vord >>;
- set!-modulus oldmod;
- vars := setkorder vord;
- x := gcdlist3(for each p in l collect reorder car p,onestep,vord);
- setkorder vars;
- return reorder x
- end;
-
- symbolic procedure gcdlist!-mod!-p(gg,l);
- if null l then gg
- else if gg=1 then 1
- else gcdlist!-mod!-p(gcd!-mod!-p(gg,car l),cdr l);
-
- symbolic procedure deg2vars l;
- if null l then nil
- else if cdar l>2 then nil
- else caar l . deg2vars cdr l;
-
- symbolic procedure vdelete(a,b);
- if null b then nil
- else if car a=caar b then cdr b
- else car b . vdelete(a,cdr b);
-
- symbolic procedure vintersection(a,b);
- begin
- scalar c;
- return if null a then nil
- else if null (c:=assoc(caar a,b)) then vintersection(cdr a,b)
- else if cdar a>cdr c then
- if cdr c=0 then vintersection(cdr a,b)
- else c . vintersection(cdr a,b)
- else if cdar a=0 then vintersection(cdr a,b)
- else car a . vintersection(cdr a,b)
- end;
-
-
- symbolic procedure vunion(a,b);
- begin
- scalar c;
- return if null a then b
- else if null (c:=assoc(caar a,b)) then car a . vunion(cdr a,b)
- else if cdar a>cdr c then car a . vunion(cdr a,delete(c,b))
- else c . vunion(cdr a,delete(c,b))
- end;
-
-
- symbolic procedure mapcarcar l;
- for each x in l collect car x;
-
-
- symbolic procedure gcdld(l,n);
- % GCD of the domain element N and all the polys in L;
- if n=1 or n=-1 then 1
- else if l=nil then abs n
- else if car l=nil then gcdld(cdr l,n)
- else gcdld(cdr l,gcd!-with!-number(n,car l));
-
- symbolic procedure split!-wrt!-variables(p,vl,l);
- % Push all the coeffs in P wrt variables in VL onto the list L
- % Stop if 1 is found as a coeff;
- if p=nil then l
- else if not null l and car l=1 then l
- else if domainp p then abs p . l
- else if member(mvar p,vl) then
- split!-wrt!-variables(red p,vl,split!-wrt!-variables(lc p,vl,l))
- else p . l;
-
- symbolic procedure gcdlist3(l,onestep,vlist);
- % GCD of the nontrivial polys in the list L given that they all
- % involve all the variables that any of them mention,
- % and they are all monomial-primitive.
- % ONESTEP is true if it is predicted that only one PRS step
- % will be needed to compute the gcd - if so try that PRS step.
- begin
- scalar unlucky!-case,image!-set,gg,gcont,l1,w,w1,w2,
- reduced!-degree!-lclst,p1,p2;
- % Make all the polys primitive;
- l1:=for each p in l collect p . ezgcd!-comfac p;
- l:=for each c in l1 collect
- quotfail1(car c,comfac!-to!-poly cdr c,
- "Content divison in GCDLIST3 failed");
- % All polys in L are now primitive.
- % Because all polys were monomial-primitive, there should
- % be no power of V to go in the result.
- gcont:=gcdlist for each c in l1 collect cddr c;
- if domainp gcont then if not(gcont=1)
- then errorf "GCONT has numeric part";
- % GCD of contents complete now;
- % Now I will remove duplicates from the list;
- % trace!-time <<
- % prin2t "GCDLIST3 on the polynomials";
- % for each p in l do print p >>;
- l := sort(for each p in l collect poly!-abs p,function ordp);
- w := nil;
- while l do <<
- w := car l . w;
- repeat l := cdr l until null l or not(car w = car l)>>;
- l := reversip w;
- w := nil;
- % trace!-time <<
- % prin2t "Made positive, with duplicates removed...";
- % for each p in l do print p >>;
- if null cdr l then return multf(gcont,car l);
- % That left just one poly;
- if domainp (gg:=car (l:=sort(l,function degree!-order))) then
- return gcont;
- % Primitive part of one poly is a constant (must be +/-1);
- if ldeg gg=1 then <<
- % True gcd is either GG or 1;
- if division!-test(gg,l) then return multf(poly!-abs gg,gcont)
- else return gcont >>;
- % All polys are now primitive and nontrivial. Use a modular
- % method to extract GCD;
- if onestep then <<
- % Try to take gcd in just one pseudoremainder step, because some
- % previous modular test suggests it may be possible;
- p1 := poly!-abs car l; p2 := poly!-abs cadr l;
- % Because polynomials are primitive and they have been normalised
- % wrt sign the only way that just one PRS step could lead to zero
- % would be if the two polys are identical. In which case that
- % should be the GCD. Note that because I got to gcdlist3 at all
- % both polys should use (all of) the same set of variables, and
- % in particular should have the same main variable.
- if p1=p2 then <<
- if division!-test(p1,cddr l) then return multf(p1,gcont)>>
- else <<
- % trace!-time prin2t "Just one pseudoremainder step needed?";
- gg := poly!-gcd(lc p1,lc p2);
- w1 := multf(red p1, quotfail1(lc p2, gg,
- "Division failure when just one pseudoremainder step needed"));
- w2 := multf(red p2,negf quotfail1(lc p1, gg,
- "Division failure when just one pseudoremainder step needed"));
- w := ldeg p1 - ldeg p2;
- if w > 0 then w2 := multf(w2, (mksp(mvar p2, w) .* 1) .+ nil)
- else if w < 0
- then w1 := multf(w1, (mksp(mvar p1, -w) .* 1) .+ nil);
- gg := ezgcd!-pp addf(w1, w2);
- % trace!-time printsf gg;
- if division!-test(gg,l) then return multf(gg,gcont) >>>>;
- return gcdlist31(l,vlist,gcont,gg,l1)
- end;
- symbolic procedure gcdlist31(l,vlist,gcont,gg,l1);
- begin scalar cofactor,lcg,old!-modulus,prime,w,w1,zeros!-list;
- old!-modulus:=set!-modulus nil; %Remember modulus;
- lcg:=for each poly in l collect lc poly;
- % trace!-time << prin2t "L.C.S OF L ARE:";
- % for each lcpoly in lcg do printsf lcpoly >>;
- lcg:=gcdlist lcg;
- % trace!-time << prin2!* "LCG (=GCD OF THESE) = ";
- % printsf lcg >>;
- try!-again:
- unlucky!-case:=nil;
- image!-set:=nil;
- set!-modulus(prime:=random!-prime());
- % Produce random univariate modular images of all the
- % polynomials;
- w:=l;
- if not zeros!-list then <<
- image!-set:=
- zeros!-list:=try!-max!-zeros!-for!-image!-set(w,vlist);
- % trace!-time << prin2t image!-set;
- % prin2 " Zeros-list = ";
- % prin2t zeros!-list >>
- >>;
- % trace!-time prin2t list("IMAGE SET",image!-set);
- gg:=make!-image!-mod!-p(car w,car vlist);
- % trace!-time prin2t list("IMAGE SET",image!-set," GG",gg);
- if unlucky!-case then <<
- % trace!-time << prin2t "Unlucky case, try again";
- % print image!-set >>;
- go to try!-again >>;
- l1:=list(car w . gg);
- make!-images:
- if null (w:=cdr w) then go to images!-created!-successfully;
- l1:=(car w . make!-image!-mod!-p(car w,car vlist)) . l1;
- if unlucky!-case then <<
- % trace!-time << prin2t "UNLUCKY AGAIN...";
- % prin2t l1;
- % print image!-set >>;
- go to try!-again >>;
- gg:=gcd!-mod!-p(gg,cdar l1);
- if domainp gg then <<
- set!-modulus old!-modulus;
- % trace!-time print "Primitive parts are coprime";
- return gcont >>;
- go to make!-images;
- images!-created!-successfully:
- l1:=reversip l1; % Put back in order with smallest first;
- % If degree of gcd seems to be same as that of smallest item
- % in input list, that item should be the gcd;
- if ldeg gg=ldeg car l then <<
- gg:=poly!-abs car l;
- % trace!-time <<
- % prin2!* "Probable GCD = ";
- % printsf gg >>;
- go to result >>
- else if (ldeg car l=add1 ldeg gg) and
- (ldeg car l=ldeg cadr l) then <<
- % Here it seems that I have just one pseudoremainder step to
- % perform, so I might as well do it;
- % trace!-time <<
- % prin2t "Just one pseudoremainder step needed"
- % >>;
- gg := poly!-gcd(lc car l,lc cadr l);
- gg := ezgcd!-pp addf(multf(red car l,
- quotfail1(lc cadr l,gg,
- "Division failure when just one pseudoremainder step needed")),
- multf(red cadr l,negf quotfail1(lc car l,gg,
- "Divison failure when just one pseudoremainder step needed")));
- % trace!-time printsf gg;
- go to result >>;
- w:=l1;
- find!-good!-cofactor:
- if null w then go to special!-case; % No good cofactor available;
- if domainp gcd!-mod!-p(gg,cofactor:=quotient!-mod!-p(cdar w,gg))
- then go to good!-cofactor!-found;
- w:=cdr w;
- go to find!-good!-cofactor;
- good!-cofactor!-found:
- cofactor:=monic!-mod!-p cofactor;
- % trace!-time prin2t "*** Good cofactor found";
- w:=caar w;
- % trace!-time << prin2!* "W= ";
- % printsf w;
- % prin2!* "GG= ";
- % printsf gg;
- % prin2!* "COFACTOR= ";
- % printsf cofactor >>;
- image!-set:=sort(image!-set,function ordopcar);
- % trace!-time << prin2 "IMAGE-SET = ";
- % prin2t image!-set;
- % prin2 "PRIME= "; prin2t prime;
- % prin2t "L (=POLYLIST) IS:";
- % for each ll in l do printsf ll >>;
- gg:=reconstruct!-gcd(w,gg,cofactor,prime,image!-set,lcg);
- if gg='nogood then go to try!-again;
- go to result;
- special!-case: % Here I have to do the first step of a PRS method;
- % trace!-time << prin2t "*** SPECIAL CASE IN GCD ***";
- % prin2t l;
- % prin2t "----->";
- % prin2t gg >>;
- reduced!-degree!-lclst:=nil;
- try!-reduced!-degree!-again:
- % trace!-time << prin2t "L1 =";
- % for each ell in l1 do print ell >>;
- w1:=reduced!-degree(caadr l1,caar l1);
- w:=car w1; w1:=cdr w1;
- if not domainp w and
- (domainp w1 or ldeg w neq ldeg w1) then go to try!-again;
- % trace!-time << prin2 "REDUCED!-DEGREE = "; printsf w;
- % prin2 " and its image = "; printsf w1 >>;
- % reduce the degree of the 2nd poly using the 1st. Result is
- % a pair : (new poly . image new poly);
- if domainp w and not null w then <<
- set!-modulus old!-modulus; return gcont >>;
- % we're done as they're coprime;
- if w and ldeg w = ldeg gg then <<
- gg:=w; go to result >>;
- % possible gcd;
- if null w then <<
- % the first poly divided the second one;
- l1:=(car l1 . cddr l1); % discard second poly;
- if null cdr l1 then <<
- gg := poly!-abs caar l1;
- go to result >>;
- go to try!-reduced!-degree!-again >>;
- % haven't made progress yet so repeat with new polys;
- if ldeg w<=ldeg gg then <<
- gg := poly!-abs w;
- go to result >>
- else if domainp gcd!-mod!-p(gg,cofactor:=quotient!-mod!-p(w1,gg))
- then <<
- w := list list w;
- go to good!-cofactor!-found >>;
- l1:= if ldeg w <= ldeg caar l1 then
- ((w . w1) . (car l1 . cddr l1))
- else (car l1 . ((w . w1) . cddr l1));
- % replace first two polys by the reduced poly and the first
- % poly ordering according to degree;
- go to try!-reduced!-degree!-again;
- % need to repeat as we still haven't found a good cofactor;
- result: % Here GG holds a tentative gcd for the primitive parts of
- % all input polys, and GCONT holds a proper one for the content;
- if division!-test(gg,l) then <<
- set!-modulus old!-modulus;
- return multf(gg,gcont) >>;
- % trace!-time prin2t list("Trial division by ",gg," failed");
- go to try!-again
- end;
- symbolic procedure make!-a!-list!-of!-variables l;
- begin scalar vlist;
- for each ll in l do vlist:=variables!.in!.form(ll,vlist);
- return make!-order!-consistent(vlist,kord!*)
- end;
-
- symbolic procedure make!-order!-consistent(l,m);
- % L is a subset of M. Make its order consistent with that
- % of M;
- if null l then nil
- else if null m then errorf("Variable missing from KORD*")
- else if car m member l then car m .
- make!-order!-consistent(delete(car m,l),cdr m)
- else make!-order!-consistent(l,cdr m);
-
- symbolic procedure try!-max!-zeros!-for!-image!-set(l,vlist);
- if null vlist then error(50,"VLIST NOT SET IN TRY-MAX-ZEROS-...")
- else begin scalar z;
- z:=for each v in cdr vlist collect
- if domainp lc car l or null quotf(lc car l,!*k2f v) then
- (v . 0) else (v . modular!-number next!-random!-number());
- for each ff in cdr l do
- z:=for each w in z collect
- if zerop cdr w then
- if domainp lc ff or null quotf(lc ff,!*k2f car w) then w
- else (car w . modular!-number next!-random!-number())
- else w;
- return z
- end;
-
- symbolic procedure
- reconstruct!-gcd(full!-poly,gg,cofactor,p,imset,lcg);
- if null addf(full!-poly,negf multf(gg,cofactor)) then gg
- else (lambda factor!-level;
- begin scalar number!-of!-factors,image!-factors,
- true!-leading!-coeffts,multivariate!-input!-poly,
- no!-of!-primes!-to!-try,
- irreducible,non!-monic,bad!-case,target!-factor!-count,
- multivariate!-factors,hensel!-growth!-size,alphalist,
- best!-known!-factors,prime!-base,
- m!-image!-variable, reconstructing!-gcd,full!-gcd;
- if not(current!-modulus=p) then
- errorf("GCDLIST HAS NOT RESTORED THE MODULUS");
- % *WARNING* GCDLIST does not restore the modulus so
- % I had better reset it here! ;
- if poly!-minusp lcg then error(50,list("Negative GCD: ",lcg));
- full!-poly:=poly!-abs full!-poly;
- initialise!-hensel!-fluids(full!-poly,gg,cofactor,p,lcg);
- % trace!-time << prin2t "TRUE LEADING COEFFTS ARE:";
- % for i:=1:2 do <<
- % printsf getv(image!-factors,i);
- % prin2!* " WITH L.C.:";
- % printsf getv(true!-leading!-coeffts,i) >> >>;
- if determine!-more!-coeffts()='done then
- return full!-gcd;
- if null alphalist then alphalist:=alphas(2,
- list(getv(image!-factors,1),getv(image!-factors,2)),1);
- if alphalist='factors! not! coprime then
- errorf list("image factors not coprime?",image!-factors);
- if not !*overview then factor!-trace <<
- printstr
- "The following modular polynomials are chosen such that:";
- terpri();
- prin2!* " a(2)*f(1) + a(1)*f(2) = 1 mod ";
- printstr hensel!-growth!-size;
- terpri();
- printstr " where degree of a(1) < degree of f(1),";
- printstr " and degree of a(2) < degree of f(2),";
- printstr " and";
- for i:=1:2 do <<
- prin2!* " a("; prin2!* i; prin2!* ")=";
- printsf cdr get!-alpha getv(image!-factors,i);
- prin2!* "and f("; prin2!* i; prin2!* ")=";
- printsf getv(image!-factors,i);
- terpri!* t >>
- >>;
- reconstruct!-multivariate!-factors(
- for each v in imset collect (car v . modular!-number cdr v));
- if irreducible or bad!-case then return 'nogood
- else return full!-gcd
- end) (factor!-level+1) ;
-
- symbolic procedure initialise!-hensel!-fluids(fpoly,fac1,fac2,p,lcf1);
- % ... ;
- begin scalar lc1!-image,lc2!-image;
- reconstructing!-gcd:=t;
- multivariate!-input!-poly:=multf(fpoly,lcf1);
- no!-of!-primes!-to!-try := 5;
- prime!-base:=hensel!-growth!-size:=p;
- number!-of!-factors:=2;
- lc1!-image:=make!-numeric!-image!-mod!-p lcf1;
- lc2!-image:=make!-numeric!-image!-mod!-p lc fpoly;
- % Neither of the above leading coefficients will vanish;
- fac1:=times!-mod!-p(lc1!-image,fac1);
- fac2:=times!-mod!-p(lc2!-image,fac2);
- image!-factors:=mkvect 2;
- true!-leading!-coeffts:=mkvect 2;
- putv(image!-factors,1,fac1);
- putv(image!-factors,2,fac2);
- putv(true!-leading!-coeffts,1,lcf1);
- putv(true!-leading!-coeffts,2,lc fpoly);
- % If the GCD is going to be monic, we know the lc
- % of both cofactors exactly;
- non!-monic:=not(lcf1=1);
- m!-image!-variable:=mvar fpoly
- end;
-
- symbolic procedure division!-test(gg,l);
- % Predicate to test if GG divides all the polynomials in the list L;
- if null l then t
- else if null quotf(car l,gg) then nil
- else division!-test(gg,cdr l);
-
-
- symbolic procedure degree!-order(a,b);
- % Order standard forms using their degrees wrt main vars;
- if domainp a then t
- else if domainp b then nil
- else ldeg a<ldeg b;
-
- symbolic procedure make!-image!-mod!-p(p,v);
- % Form univariate image, set UNLUCKY!-CASE if leading coefficient
- % gets destroyed;
- begin
- scalar lp;
- lp := degree!-in!-variable(p,v);
- p := make!-univariate!-image!-mod!-p(p,v);
- if not(degree!-in!-variable(p,v)=lp) then unlucky!-case := t;
- return p
- end;
-
-
- symbolic procedure make!-univariate!-image!-mod!-p(p,v);
- % Make a modular image of P, keeping only the variable V;
- if domainp p then
- if p=nil then nil
- else !*n2f modular!-number p
- else if mvar p=v then
- adjoin!-term(lpow p,
- make!-univariate!-image!-mod!-p(lc p,v),
- make!-univariate!-image!-mod!-p(red p,v))
- else plus!-mod!-p(
- times!-mod!-p(image!-of!-power(mvar p,ldeg p),
- make!-univariate!-image!-mod!-p(lc p,v)),
- make!-univariate!-image!-mod!-p(red p,v));
-
- symbolic procedure image!-of!-power(v,n);
- begin
- scalar w;
- w := assoc(v,image!-set);
- if null w then <<
- w := modular!-number next!-random!-number();
- image!-set := (v . w) . image!-set >>
- else w := cdr w;
- return modular!-expt(w,n)
- end;
-
- symbolic procedure make!-numeric!-image!-mod!-p p;
- % Make a modular image of P;
- if domainp p then
- if p=nil then 0
- else modular!-number p
- else modular!-plus(
- modular!-times(image!-of!-power(mvar p,ldeg p),
- make!-numeric!-image!-mod!-p lc p),
- make!-numeric!-image!-mod!-p red p);
-
-
- symbolic procedure total!-degree!-in!-powers(form,powlst);
- % Returns a list where each variable mentioned in FORM is paired
- % with the maximum degree it has. POWLST collects the list, and should
- % normally be NIL on initial entry;
- if null form or domainp form then powlst
- else begin scalar x;
- if (x := atsoc(mvar form,powlst))
- then ldeg form>cdr x and rplacd(x,ldeg form)
- else powlst := (mvar form . ldeg form) . powlst;
- return total!-degree!-in!-powers(red form,
- total!-degree!-in!-powers(lc form,powlst))
- end;
-
-
- symbolic procedure powers1 form;
- % For each variable V in FORM collect (V . (MAX . MIN)) where
- % MAX and MIN are limits to the degrees V has in FORM;
- powers2(form,powers3(form,nil),nil);
-
- symbolic procedure powers3(form,l);
- % Start of POWERS1 by collecting power information for
- % the leading monomial in FORM;
- if domainp form then l
- else powers3(lc form,(mvar form . (ldeg form . ldeg form)) . l);
-
- symbolic procedure powers2(form,powlst,thismonomial);
- if domainp form then
- if null form then powlst else powers4(thismonomial,powlst)
- else powers2(lc form,
- powers2(red form,powlst,thismonomial),
- lpow form . thismonomial);
-
- symbolic procedure powers4(new,old);
- % Merge information from new monomial into old information,
- % updating MAX and MIN details;
- if null new then for each v in old collect (car v . (cadr v . 0))
- else if null old then for each v in new collect (car v . (cdr v . 0))
- else if caar new=caar old then <<
- % variables match - do MAX and MIN on degree information;
- if cdar new>cadar old then rplaca(cdar old,cdar new);
- if cdar new<cddar old then rplacd(cdar old,cdar new);
- rplacd(old,powers4(cdr new,cdr old)) >>
- else if ordop(caar new,caar old) then <<
- rplacd(cdar old,0); % Some variable not mentioned in new monomial;
- rplacd(old,powers4(new,cdr old)) >>
- else (caar new . (cdar new . 0)) . powers4(cdr new,old);
-
-
- symbolic procedure ezgcd!-pp u;
- %returns the primitive part of the polynomial U wrt leading var;
- quotf1(u,comfac!-to!-poly ezgcd!-comfac u);
-
- symbolic procedure ezgcd!-sqfrf p;
- %P is a primitive standard form;
- %value is a list of square free factors;
- begin
- scalar pdash,p1,d,v;
- pdash := diff(p,v := mvar p);
- d := poly!-gcd(p,pdash); % p2*p3**2*p4**3*... ;
- if domainp d then return list p;
- p := quotfail1(p,d,"GCD division in FACTOR-SQFRF failed");
- p1 := poly!-gcd(p,
- addf(quotfail1(pdash,d,"GCD division in FACTOR-SQFRF failed"),
- negf diff(p,v)));
- return p1 . ezgcd!-sqfrf d
- end;
- symbolic procedure reduced!-degree(u,v);
- %U and V are primitive polynomials in the main variable VAR;
- %result is pair: (reduced poly of U by V . its image) where by
- % reduced I mean using V to kill the leading term of U;
- begin scalar var,w,x;
- % trace!-time << prin2t "ARGS FOR REDUCED!-DEGREE ARE:";
- % printsf u; printsf v >>;
- if u=v or quotf1(u,v) then return (nil . nil)
- else if ldeg v=1 then return (1 . 1);
- % trace!-time prin2t "CASE NON-TRIVIAL SO TAKE A REDUCED!-DEGREE:";
- var := mvar u;
- if ldeg u=ldeg v then x := negf lc u
- else x:=(mksp(var,ldeg u - ldeg v) .* negf lc u) .+ nil;
- w:=addf(multf(lc v,u),multf(x,v));
- % trace!-time printsf w;
- if degr(w,var)=0 then return (1 . 1);
- % trace!-time << prin2 "REDUCED!-DEGREE-LCLST = ";
- % print reduced!-degree!-lclst >>;
- reduced!-degree!-lclst := addlc(v,reduced!-degree!-lclst);
- % trace!-time << prin2 "REDUCED!-DEGREE-LCLST = ";
- % print reduced!-degree!-lclst >>;
- if x := quotf1(w,lc w) then w := x
- else for each y in reduced!-degree!-lclst do
- while (x := quotf1(w,y)) do w := x;
- u := v; v := ezgcd!-pp w;
- % trace!-time << prin2t "U AND V ARE NOW:";
- % printsf u; printsf v >>;
- if degr(v,var)=0 then return (1 . 1)
- else return (v . make!-univariate!-image!-mod!-p(v,var))
- end;
-
-
- % moved('comfac,'ezgcd!-comfac);
- % moved('pp,'ezgcd!-pp);
-
- endmodule;
- end;
|