reduce-historical/r38/packages/factor/ezgcdf.red

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;