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