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

module modpoly;   % Routines for performing arithmetic on multivariate
                  % polynomials with coefficients that are modular
                  % numbers as defined by modular!-plus etc.

% Authors: A. C. Norman and P. M. A. Moore, 1979.

fluid '(current!-modulus
        exact!-quotient!-flag
        m!-image!-variable
        modulus!/2
        reduction!-count);

% Note that the datastructure used is the same as that used in
% REDUCE except that it is assumed that domain elements are atomic.

symbolic smacro procedure comes!-before(p1,p2);
   % Similar to the REDUCE function ORDPP, but does not cater for non-
   % commutative terms and assumes that exponents are small integers.
    (car p1=car p2 and igreaterp(cdr p1,cdr p2)) or
       (not(car p1=car p2) and ordop(car p1,car p2));

symbolic procedure plus!-mod!-p(a,b);
   % form the sum of the two polynomials a and b working over the
   % ground domain defined by the routines % modular!-plus,
   % modular!-times etc. the inputs to this % routine are assumed to
   % have coefficients already in the required domain.
   if null a then b
   else if null b then a
   else if domainp a then
      if domainp b then !*n2f modular!-plus(a,b)
      else (lt b) .+ plus!-mod!-p(a,red b)
   else if domainp b then (lt a) .+ plus!-mod!-p(red a,b)
   else if lpow a = lpow b then
      adjoin!-term(lpow a,
         plus!-mod!-p(lc a,lc b),plus!-mod!-p(red a,red b))
   else if comes!-before(lpow a,lpow b) then
         (lt a) .+ plus!-mod!-p(red a,b)
   else (lt b) .+ plus!-mod!-p(a,red b);

symbolic procedure times!-mod!-p(a,b);
   if (null a) or (null b) then nil
   else if domainp a then multiply!-by!-constant!-mod!-p(b,a)
   else if domainp b then multiply!-by!-constant!-mod!-p(a,b)
   else if mvar a=mvar b then plus!-mod!-p(
     plus!-mod!-p(times!-term!-mod!-p(lt a,b),
                  times!-term!-mod!-p(lt b,red a)),
     times!-mod!-p(red a,red b))
   else if ordop(mvar a,mvar b) then
     adjoin!-term(lpow a,times!-mod!-p(lc a,b),times!-mod!-p(red a,b))
   else adjoin!-term(lpow b,
        times!-mod!-p(a,lc b),times!-mod!-p(a,red b));

symbolic procedure times!-term!-mod!-p(term,b);
   % Multiply the given polynomial by the given term.
    if null b then nil
    else if domainp b then
        adjoin!-term(tpow term,
            multiply!-by!-constant!-mod!-p(tc term,b),nil)
    else if tvar term=mvar b then
	 adjoin!-term(mksp(tvar term,iplus2(tdeg term,ldeg b)),
                      times!-mod!-p(tc term,lc b),
                      times!-term!-mod!-p(term,red b))
    else if ordop(tvar term,mvar b) then
      adjoin!-term(tpow term,times!-mod!-p(tc term,b),nil)
    else adjoin!-term(lpow b,
      times!-term!-mod!-p(term,lc b),
      times!-term!-mod!-p(term,red b));

symbolic procedure difference!-mod!-p(a,b);
   plus!-mod!-p(a,minus!-mod!-p b);

symbolic procedure minus!-mod!-p a;
   if null a then nil
   else if domainp a then modular!-minus a
   else (lpow a .* minus!-mod!-p lc a) .+ minus!-mod!-p red a;


symbolic procedure reduce!-mod!-p a;
   % Converts a multivariate poly from normal into modular polynomial.
    if null a then nil
    else if domainp a then !*n2f modular!-number a
    else adjoin!-term(lpow a,reduce!-mod!-p lc a,reduce!-mod!-p red a);

symbolic procedure 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 MONIC-MOD-P"
   else multiply!-by!-constant!-mod!-p(a,
     modular!-reciprocal lc a);

symbolic procedure quotfail!-mod!-p(a,b);
   % Form quotient A/B, but complain if the division is not exact.
  begin scalar c;
    exact!-quotient!-flag:=t;
    c:=quotient!-mod!-p(a,b);
    if exact!-quotient!-flag then return c
    else errorf "QUOTIENT NOT EXACT (MOD P)"
  end;

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


symbolic procedure 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:=quotient!-mod!-p(lc a,lc b);
      if difference!-mod!-p(a,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)) .*
        quotient!-mod!-p(lc a,lc b);
   % That is the leading term of the quotient.  Now subtract term*b from
   % a.
      a:=plus!-mod!-p(red a,
                      times!-term!-mod!-p(negate!-term term,red b));
   % or a:=a-b*term given leading terms must cancel.
      return term .+ xquotient!-mod!-p(a,b,v)
    end;

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


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


symbolic procedure 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:=quotient!-mod!-p(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:=plus!-mod!-p(red a,
      multiply!-by!-constant!-mod!-p(red b,q))
    else
      a:=plus!-mod!-p(red a,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 xremainder!-mod!-p(a,b,v)
   end;

symbolic procedure 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 modular!-times(a,n)
   else adjoin!-term(lpow a,multiply!-by!-constant!-mod!-p(lc a,n),
     multiply!-by!-constant!-mod!-p(red a,n));

symbolic procedure 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
      ordered!-gcd!-mod!-p(a,b)
    else ordered!-gcd!-mod!-p(b,a) >>;

symbolic procedure ordered!-gcd!-mod!-p(a,b);
   % As above, but deg a > deg b.
  begin scalar steps;
    steps := 0;
top:
    a := reduce!-degree!-mod!-p(a,b);
    if null a then return 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 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:=modular!-quotient(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 plus!-mod!-p(red a,
      multiply!-by!-constant!-mod!-p(red b,q))
    else return plus!-mod!-p(red a,times!-term!-mod!-p(
            mksp(mvar b,w) .* q,red b))
   end;

symbolic procedure derivative!-mod!-p a;
   % Derivative of a wrt its main variable.
   if domainp a then nil
   else if ldeg a=1 then lc a
   else derivative!-mod!-p!-1(a,mvar a);

symbolic procedure derivative!-mod!-p!-1(a,v);
    if domainp a then nil
    else if not(mvar a=v) then nil
    else if ldeg a=1 then lc a
   else adjoin!-term(mksp(v,isub1 ldeg a),
                 multiply!-by!-constant!-mod!-p(lc a,
                                                modular!-number ldeg a),
                 derivative!-mod!-p!-1(red a,v));

symbolic procedure square!-free!-mod!-p a;
   % predicate that tests if a is square-free as a modular univariate
   % polynomial.
    domainp a or domainp gcd!-mod!-p(a,derivative!-mod!-p a);

symbolic procedure evaluate!-mod!-p(a,v,n);
   % Evaluate polynomial A at the point V=N.
    if domainp a then a
    else if n=0 then evaluate!-mod!-p(a,v,nil)
    else if v=nil then errorf "Variable=NIL in EVALUATE-MOD-P"
    else if mvar a=v then horner!-rule!-mod!-p(lc a,ldeg a,red a,n,v)
    else adjoin!-term(lpow a,
      evaluate!-mod!-p(lc a,v,n),
      evaluate!-mod!-p(red a,v,n));
 
symbolic procedure 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:=times!-mod!-p(v,expt!-mod!-p(n,degg));
                   plus!-mod!-p(a,v)>>
    else begin scalar newdeg;
      newdeg:=ldeg a;
      return horner!-rule!-mod!-p(if null n or zerop n then lc a
                                   else plus!-mod!-p(lc a,
         times!-mod!-p(v,expt!-mod!-p(n,idifference(degg,newdeg)))),
       newdeg,red a,n,var)
    end;

symbolic procedure 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:=expt!-mod!-p(a,car w);
     x:=times!-mod!-p(x,x);
     if not (cdr w = 0) then x:=times!-mod!-p(x,a);
     return x
    end;

symbolic procedure make!-bivariate!-mod!-p(u,imset,v);
   % Substitute into U for all variables in IMSET which should result in
   % a bivariate poly. One variable is M-IMAGE-VARIABLE and V is the
   % other U is modular multivariate with these two variables at top 2
   % levels - V at 2nd level.
  if domainp u then u
  else if mvar u = m!-image!-variable then
    adjoin!-term(lpow u,make!-univariate!-mod!-p(lc u,imset,v),
      make!-bivariate!-mod!-p(red u,imset,v))
  else make!-univariate!-mod!-p(u,imset,v);

symbolic procedure make!-univariate!-mod!-p(u,imset,v);
   % Substitute into U for all variables in IMSET giving a univariate
   % poly in V. U is modular multivariate with V at top level.
  if domainp u then u
  else if mvar u = v then
    adjoin!-term(lpow u,!*n2f evaluate!-in!-order!-mod!-p(lc u,imset),
      make!-univariate!-mod!-p(red u,imset,v))
  else !*n2f evaluate!-in!-order!-mod!-p(u,imset);

symbolic procedure evaluate!-in!-order!-mod!-p(u,imset);
   % Makes an image of u wrt imageset, imset, using Horner's rule.
   % Result should be purely numeric (and modular).
  if domainp u then !*d2n u
   else if mvar u=caar imset then
    horner!-rule!-in!-order!-mod!-p(
      evaluate!-in!-order!-mod!-p(lc u,cdr imset),ldeg u,red u,imset)
   else evaluate!-in!-order!-mod!-p(u,cdr imset);

symbolic procedure horner!-rule!-in!-order!-mod!-p(c,degg,a,vset);
   % C is running total and a is what is left.
  if domainp a then modular!-plus(!*d2n a,
    modular!-times(c,modular!-expt(cdar vset,degg)))
  else if not(mvar a=caar vset) then
    modular!-plus(
      evaluate!-in!-order!-mod!-p(a,cdr vset),
      modular!-times(c,modular!-expt(cdar vset,degg)))
  else begin scalar newdeg;
    newdeg:=ldeg a;
    return horner!-rule!-in!-order!-mod!-p(
      modular!-plus(
        evaluate!-in!-order!-mod!-p(lc a,cdr vset),
        modular!-times(c,
          modular!-expt(cdar vset,(idifference(degg,newdeg))))),
      newdeg,red a,vset)
  end;

symbolic procedure make!-modular!-symmetric a;
   % Input is a multivariate MODULAR poly A with nos in range 0->(p-1).
   % This folds it onto the symmetric range (-p/2)->(p/2).
    if null a then nil
    else if domainp a then
      if a>modulus!/2 then !*n2f(a - current!-modulus)
      else a
    else adjoin!-term(lpow a,make!-modular!-symmetric lc a,
      make!-modular!-symmetric red a);

endmodule;

end;