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

module modsr; % Modular Solve and Roots.

% Author: Herbert Melenk, ZIB Berlin, Jan 95.

create!-package('(modsr modsqrt modroots modsolve),'(solve));

fluid '(current!-modulus);

% Some routines from solve and factor(modpoly) are needed.

load!-package 'solve;
load!-package 'factor;

% Now a few things that MIGHT have been in the factorizer but were not
% It is quite possible that as a matter of style these few functions
% should be put into factor.red, even though they are not used there,
% since that way they live near their friends and are more generally
% useful???

symbolic procedure general!-evaluate!-mod!-p(a,v,n);
   % Evaluate polynomial A at the point V=N.
    if domainp a then a
    else if n=0 then general!-evaluate!-mod!-p(a,v,nil)
    else if v=nil then errorf "Variable=NIL in GENERAL-EVALUATE-MOD-P"
    else if mvar a=v
     then general!-horner!-rule!-mod!-p(lc a,ldeg a,red a,n,v)
    else adjoin!-term(lpow a,
      general!-evaluate!-mod!-p(lc a,v,n),
      general!-evaluate!-mod!-p(red a,v,n));

symbolic procedure general!-horner!-rule!-mod!-p(v,degg,a,n,var);
   % V is the running total, and it must be multiplied by n**deg and
   % added to the value of a at n.
    if domainp a or not(mvar a=var)
      then if null n or zerop n then a
	    else <<v:=general!-times!-mod!-p(v,
			 general!-expt!-mod!-p(n,degg));
                   general!-plus!-mod!-p(a,v)>>
    else begin scalar newdeg;
      newdeg:=ldeg a;
      return general!-horner!-rule!-mod!-p(
          if null n or zerop n then lc a
          else general!-plus!-mod!-p(lc a,
             general!-times!-mod!-p(v,
                general!-expt!-mod!-p(n,idifference(degg,newdeg)))),
       newdeg,red a,n,var)
    end;

symbolic procedure general!-expt!-mod!-p(a,n);
   % a**n.
    if n=0 then 1
    else if n=1 then a
    else begin scalar w,x;
     w:=divide(n,2);
     x:=general!-expt!-mod!-p(a,car w);
     x:=general!-times!-mod!-p(x,x);
     if not (cdr w = 0) then x:=general!-times!-mod!-p(x,a);
     return x
    end;

symbolic procedure general!-monic!-mod!-p a;
   % This procedure can only cope with polys that have a numeric
   % leading coeff.
   if a=nil then nil
   else if domainp a then 1
   else if lc a = 1 then a
   else if not domainp lc a then
       errorf "LC NOT NUMERIC IN GENERAL-MONIC-MOD-P"
   else general!-multiply!-by!-constant!-mod!-p(a,
     general!-modular!-reciprocal lc a);

symbolic procedure general!-quotient!-mod!-p(a,b);
   % Truncated quotient of a by b.
    if null b then errorf "B=0 IN GENERAL-QUOTIENT-MOD-P"
    else if domainp b then general!-multiply!-by!-constant!-mod!-p(a,
                             general!-modular!-reciprocal b)
    else if a=nil then nil
    else if domainp a then exact!-quotient!-flag:=nil
    else if mvar a=mvar b then general!-xquotient!-mod!-p(a,b,mvar b)
    else if ordop(mvar a,mvar b) then
       adjoin!-term(lpow a,
          general!-quotient!-mod!-p(lc a,b),
          general!-quotient!-mod!-p(red a,b))
    else exact!-quotient!-flag:=nil;


symbolic procedure general!-xquotient!-mod!-p(a,b,v);
   % Truncated quotient a/b given that b is nontrivial.
    if a=nil then nil
    else if (domainp a) or (not(mvar a=v)) or
      ilessp(ldeg a,ldeg b) then exact!-quotient!-flag:=nil
    else if ldeg a = ldeg b then begin scalar w;
      w:=general!-quotient!-mod!-p(lc a,lc b);
      if general!-difference!-mod!-p(a,general!-times!-mod!-p(w,b)) then
        exact!-quotient!-flag:=nil;
      return w
      end
    else begin scalar term;
      term:=mksp(mvar a,idifference(ldeg a,ldeg b)) .*
        general!-quotient!-mod!-p(lc a,lc b);
   % That is the leading term of the quotient.  Now subtract term*b from
   % a.
      a:=general!-plus!-mod!-p(red a,
		general!-times!-term!-mod!-p(general!-negate!-term term,
					     red b));
   % or a:=a-b*term given leading terms must cancel.
      return term .+ general!-xquotient!-mod!-p(a,b,v)
    end;

symbolic procedure general!-negate!-term term;
   % Negate a term.
    tpow term .* general!-minus!-mod!-p tc term;


symbolic procedure general!-remainder!-mod!-p(a,b);
   % Remainder when a is divided by b.
    if null b then errorf "B=0 IN GENERAL-REMAINDER-MOD-P"
    else if domainp b then nil
    else if domainp a then a
    else general!-xremainder!-mod!-p(a,b,mvar b);


symbolic procedure general!-xremainder!-mod!-p(a,b,v);
   % Remainder when the modular polynomial a is divided by b, given that
   % b is non degenerate.
   if (domainp a) or (not(mvar a=v)) or ilessp(ldeg a,ldeg b) then a
   else begin
    scalar q,w;
    q:=general!-quotient!-mod!-p(general!-minus!-mod!-p lc a,lc b);
   % compute -lc of quotient.
    w:=idifference(ldeg a,ldeg b); %ldeg of quotient;
    if w=0 then a:=general!-plus!-mod!-p(red a,
      general!-multiply!-by!-constant!-mod!-p(red b,q))
    else
      a:=general!-plus!-mod!-p(red a,general!-times!-term!-mod!-p(
            mksp(mvar b,w) .* q,red b));
   % The above lines of code use red a and red b because by construc-
   % tion the leading terms of the required % answers will cancel out.
     return general!-xremainder!-mod!-p(a,b,v)
   end;

symbolic procedure general!-multiply!-by!-constant!-mod!-p(a,n);
   % Multiply the polynomial a by the constant n.
   if null a then nil
   else if n=1 then a
   else if domainp a then !*n2f general!-modular!-times(a,n)
   else adjoin!-term(lpow a,
		     general!-multiply!-by!-constant!-mod!-p(lc a,n),
     general!-multiply!-by!-constant!-mod!-p(red a,n));

symbolic procedure general!-gcd!-mod!-p(a,b);
   % Return the monic gcd of the two modular univariate polynomials a
   % and b.  Set REDUCTION-COUNT to the number of steps taken in the
   % process.
 << reduction!-count := 0;
    if null a then monic!-mod!-p b
    else if null b then monic!-mod!-p a
    else if domainp a then 1
    else if domainp b then 1
    else if igreaterp(ldeg a,ldeg b) then
      general!-ordered!-gcd!-mod!-p(a,b)
    else general!-ordered!-gcd!-mod!-p(b,a) >>;

symbolic procedure general!-ordered!-gcd!-mod!-p(a,b);
   % As above, but deg a > deg b.
  begin scalar steps;
    steps := 0;
top:
    a := general!-reduce!-degree!-mod!-p(a,b);
    if null a then return general!-monic!-mod!-p b;
    steps := steps + 1;
    if domainp a then <<
        reduction!-count := reduction!-count+steps;
        return 1 >>
    else if ldeg a<ldeg b then begin
      scalar w;
      reduction!-count := reduction!-count + steps;
      steps := 0;
      w := a; a := b; b := w
      end;
    go to top
  end;

symbolic procedure general!-reduce!-degree!-mod!-p(a,b);
   % Compute A-Q*B where Q is a single term chosen so that the result
   % has lower degree than A did.
  begin
    scalar q,w;
    q:=general!-modular!-quotient(general!-modular!-minus lc a,lc b);
   % compute -lc of quotient;
    w:=idifference(ldeg a,ldeg b); %ldeg of quotient;
   % The next lines of code use red a and red b because by construction
   % the leading terms of the required answers will cancel out.
    if w=0 then return general!-plus!-mod!-p(red a,
      general!-multiply!-by!-constant!-mod!-p(red b,q))
    else return general!-plus!-mod!-p(red a,
		   general!-times!-term!-mod!-p(mksp(mvar b,w) .* q,
						red b))
   end;


%%%%%%%

symbolic procedure modp(a,p);
   <<a:=remainder(a,p); if a<0 then a+p else a>>;

symbolic procedure lowestdeg(f,x,n);
   if null f then n else
   if domainp f or mvar f neq x then 0 else
   lowestdeg(red f,x,ldeg f);

symbolic procedure reduce!-mod!-p!*(f,p);
  (general!-reduce!-mod!-p f) where current!-modulus = p;

symbolic procedure moduntag f;
  if eqcar(f,'!:mod!:) then cdr f else
  if atom f then f else
  moduntag car f . moduntag cdr f;

symbolic procedure safe!-modrecip u;
   % Return 1/u or nil.
   begin scalar q,!*msg,!*protfg;
      !*msg := nil; !*protfg := t;
      if eqcar(u,'!:mod!:) then u := cdr u;
      q := errorset({'general!-modular!-reciprocal, u},nil,nil);
      erfg!* := nil;
      return if errorp q then nil else car q
   end;

endmodule;

end;