reduce-historical/r37/packages/solve/modroots.red

module modroots; % Roots of a univariate polynomial mod m,
                 % m not necessarily prime.

% Author: Herbert Melenk, ZIB Berlin.

% Algebraic interface: m_roots(polynomial, modulus);

symbolic procedure modroots0(f,m);
 % f: univariate standard form with modular coeffients,
 % m: positive integer modulus.
 % Algorithm: compute roots modulo the biggest factor of m
 % and lift these for the remaining factors. During lifing
 % the number of factors may change in both directions.
 begin scalar ml;
   ml := sort(for each q in zfactor m join
                for i:=1:cdr q collect car q,'lessp);
   return sort(modroots1(f,ml),'lessp);
 end;

symbolic procedure modroots1(f,ml);
   if null cdr ml then modroots2(f,car ml,nil) else
  begin scalar f1,p,q,pq,r,s,x,y;
    p:=car ml; ml:=cdr ml;
    r := modroots1(f,ml);
    if null r then return nil;
    x:=mvar f; y:=gensym();
    q:=for each m in ml product m;
    pq:=p*q;
     % lift roots to p*q:
     %   if f(r)=0 mod q, solve f(n*q+r)=0 mod p.
    for each w in r do
    <<f1:=numr subf(f,{x . {'plus,{'times,y,q},w}});
      for each y in modroots2(reduce!-mod!-p!*(f1,p),p,t)
        do <<y:= modp(y*q+w,pq);
             if null reduce!-mod!-p!*(numr subf(f,{x . y}),pq)
               and not member(y,s) then s:=y.s>>;
    >>;
    return s;
  end;

symbolic procedure modroots2(f,p,rec);
  if domainp f and f then nil
  else if null f then
    if p=2 and rec then '(-1 0 1)
      else for i:=0:(p-1) collect i
  else if p=2 then modroots4(f,t,rec)
  else modroots3(f,p);

symbolic procedure x!*!*p!-w(x,p,w);
  % Make a form x^p - w mod p.
    general!-difference!-mod!-p(x .** p .*1 .+ nil,w);

symbolic procedure modroots3(f,current!-modulus);
 % Roots of a polynomial f mod p, p prime.
 % Algorithm:
 % H. Cohen: Computational Algebraic Number theory, 1.6.1
 begin scalar a,d,p,r,x;  integer n;

     % From now on, we compute with untagged modular coefficients
     % using the routines in "factor/modpoly".
   p := current!-modulus;
   f := general!-reduce!-mod!-p f;
   x := mvar f;

     % gcd(f, x^p - x)
   a := general!-gcd!-mod!-p(f , x!*!*p!-w(x,p,!*k2f x));
   d := ldeg a;
   n := lowestdeg(a,x,0);
   if n>0 then
   <<r:='(0); a:=general!-quotient!-mod!-p(a,x!*!*p!-w(x,n,nil))>>;
   return append(r,modroots31(a,x,p));
 end;

symbolic procedure modroots31(a,x,p);
  begin scalar a0,a1,a2,b,d,e,s,w;

  s2:
   if domainp a then return nil;
   if ldeg a = 1 then
      return {general!-modular!-quotient(
                 if red a then general!-modular!-minus red a else 0,
                 lc a)};
   if ldeg a = 2 then
   <<
     a2:=lc a; a:=red a;
     if not domainp a then
     <<a1:= lc a; a:=red a>> else a1:=0;
     a0:=if null a then 0 else a;
     d:=general!-modular!-difference(
           general!-modular!-times(a1,a1),
           general!-modular!-times(4,general!-modular!-times(a0,a2)));
     s:=legendre!-symbol(d,p);
     if s=-1 then return nil;
     e:= modsqrt(d,p);
     a2:=general!-modular!-reciprocal general!-modular!-plus(a2,a2);
     a1:=general!-modular!-minus a1;
     return
      {general!-modular!-times(general!-modular!-plus(a1,e),a2),
       general!-modular!-times(general!-modular!-difference(a1,e),a2)};
   >>;

   s3:
    e:=random(p);
      % compute gcd[x ^((p-1)/2) - 1, A(x - e)]
    w:=x!*!*p!-w(x,(p-1)/2,1);
    a1:=general!-reduce!-mod!-p numr subf(a,{x.{'difference,x,e}});
    b:=general!-gcd!-mod!-p(w,a1);
    if domainp b or ldeg b = ldeg a then go to s3;

   s4:
     % Compute both root groups and transform roots back to x - e;
    return
    for each w in union(modroots31(general!-quotient!-mod!-p(a1,b),x,p),
			modroots31(b,x,p))
       collect general!-modular!-difference(w,e)
 end;

symbolic procedure modroots4(f,w,rec);
  % roots of f mod 2: count terms.
  if domainp f then
  <<
   if f then w:=not w;
   append(
    if null f then '(0),
    if w then (if rec then '(-1 1) else '(1))
         )
  >>
  else modroots4(red f,not w,rec);

put('m_roots,'psopfn,
      function(lambda(u);
               'list . modroots0(numr simp car u,reval cadr u)));

endmodule;

end;