module rsolve;  % Rational and integer zeros of a univariate polynomial
                % using fast modular methods.

% Author: F.J.Wright@Maths.QMW.ac.uk
% Version 1.051, 16 Jan 1995

% ***** NOT QUITE COMPATIBLE WITH VERSION 1, 6 Jul 1994 *****

% More robust, compatible with solve, minor bug fixes.
% Multiplicity also computed.
% More flexible interface.

% Note the essential temporary patch at the end of this file,
% necessary with versions of REDUCE up to and including 3.5.

% DOCUMENTATION
% =============

% Defines operators r_solve and i_solve that take a single univariate
% polynomial (or polynomial equation) as argument, and optionally the
% variable as second argument, and return respectively the sets of
% rational and integer zeros.  Any denominator is completely ignored!
% See the test/demo file rsolve.tst for examples.

% Default output format is the same as used by solve (including
% respecting the multiplicities switch), but optional arguments allow
% other output formats (see below).  Solutions of degenerate equations
% are expressed by r_solve and i_solve using the operators ARBRAT
% (which is new) and ARBINT respectively.

% Computing only the integer zeros is slightly more efficient than
% extracting them from the rational zeros.  This algorithm appears to
% be faster than solve by a factor that depends on the example, but
% typically up to about 2, and gives more convenient output if only
% integer or rational zeros are required.

% The algorithm used is that described by R. Loos (1983):
% Computing rational zeros of integral polynomials by $p$-adic
% expansion.  {\it SIAM J. Computing}. {\bf 12}, 286--293.

% The switch TRSOLVE turns on tracing of the algorithm.


% Possible extensions:
%   improve root bound used;
%   map option;
%   coupled multivariate polynomials;
%   Gaussian integer and rational zeros.


% Note that modular methods do not appear to be helpful for computing
% integer square roots.  The function isqrt in arith.red, which
% effectively returns isqrt m = floor sqrt m computed using an integer
% Newton iteration, is significantly faster than any modular method
% that I have tried to develop!

create!-package('(rsolve),'(solve));

% Algebraic interface
% ===================

switch multiplicities;

global '(!!arbint multiplicities!*);
share multiplicities!*;

put('arbrat,'simpfn,'simpiden);         % Arbitrary rational operator

fluid '(!*i_solve);

put('i_solve, 'psopfn, 'i_solve!-eval);
%put('isolve, 'psopfn, 'i_solve!-eval);  % Later, maybe

symbolic procedure i_solve!-eval u;
   % Compute only integer zeros of a univariate polynomial.
   r_solve!-eval u where !*i_solve = t;

put('r_solve, 'psopfn, 'r_solve!-eval);

symbolic procedure r_solve!-eval u;
   % Algebraic interface.  Sets integer domain mode internally.
   %
   % Use: {R/I}_Solve(eqn, var, options)
   %
   % eqn: univariate equation with rational coefficients (required).
   % var: solution variable (optional).
   % options: optional keywords acting as local switches:-
   %   separate: (default) assign multiplicity list to the
   %     global variable root_multiplicities;
   %   expand or multiplicities: (default if multiplicities switch
   %     is on) expand the solution list to include multiple zeros
   %     multiple times;
   %   together: return each solution as a list whose second
   %     element is the multiplicity;
   %   nomul: do not compute multiplicities (thereby saving some time);
   %   noeqs: do not return univariate zeros as equations.

   if null u then rederr "r/i_solve called with no equations" else
   begin scalar var, mul, noeqs;
      % Process options, after setting defaults:
      mul := if !*multiplicities then 'expand else 'separate;
      begin scalar maybe_var;  maybe_var := t;
         for each x in cdr u do <<
            if x eq 'separate then mul := 'separate
            else if x memq '(expand multiplicities) then
               mul := 'expand
            else if x eq 'together then mul := t
            else if x eq 'nomul then mul := nil
            else if x eq 'noeqs then noeqs := t
            else if maybe_var then var := x
            else typerr(x, "optional r/i_solve argument");
            maybe_var := nil >>
      end;
      return
      (begin scalar result;
         % Domain modes have nasty global properties, so ...
         if old_dmode then setdmode(old_dmode, nil) where !*msg = nil;
         result := errorset!*({'r_solve!-eval1,
            mkquote car u, mkquote var, mkquote mul, mkquote noeqs}, t);
         if old_dmode then setdmode(old_dmode, t) where !*msg = nil;
         % The errorset outputs any error messages, so exit quietly:
         if errorp result then error1() else return car result
      end) where old_dmode = dmode!*;
    end;


% Tracing facility:
switch trsolve;

%% rlistat '(tr_write);
deflist('((tr_write rlis)), 'stat);  % flagged eval!

symbolic procedure tr_write u;
   if !*trsolve then
      << for each el in u do prin2 el; terpri() >>;

symbolic procedure r_solve!-eval1(f, var, mul, noeqs);
   % Algebraic interface to r_solve, ASSUMING integer domain mode.
   begin scalar x, u;
      f := numr simp!* !*eqn2a reval f; % Should check denr?
      tr_write "Simplified poly to solve: ", f;
      if not atom f then
         % Get mainvar and check f really is a polynomial:
         if eqcar(x := mvar f, 'list) then goto error;
      if var then                       % variable specified
         if not((var := !*a2k var) eq x) then
            if x then
               if not smember(var, f) then return makelist nil
               else goto error          % variable mis-match
            else if null f then <<      % degenerate equation
               multiplicities!* := makelist nil;
               return makelist { makeqn!-maybe(var,
                  {if !*i_solve then 'arbint else 'arbrat,
                     !!arbint:=!!arbint+1}, noeqs) } >>;
      if atom f then return makelist nil;
      u := f;
      while not atom u do
         u := if mvar u eq x and fixp lc u then red u else 'error;
      if not(u eq 'error) then
         return r_solve!-output(f, r_solve f, mul, noeqs);
   error:
      %% TypErr(prepf f, "univariate poly over Z in specified var")
      typerr(prepf f, "univariate polynomial over Z")
   end;


symbolic procedure r_solve!-output(f, zeros, mul, noeqs);
   % Assumes f is a univariate integer polynomial in standard form
   % and zeros is a list of integers or rational SQ forms.
   begin scalar x;  x := mvar f;
      return makelist
      if null mul then <<
         multiplicities!* := makelist nil;
         for each s in zeros collect makeqn!-maybe(x, !*n2a s, noeqs) >>
      else
      % Compute multiplicities and return composite list.  This simple
      % successive-differentiation algorithm is the most efficient for
      % zeros of multiplicity less than about 5, because it requires
      % essentially m polynomial evaluations per zero of multiplicty
      % m, whereas the next best algorithm based on evaluating
      % (x-a)f'/f at x=a requires one exact division, one gcd, 2
      % polynomial evaluations and one numerical division per zero.
      begin scalar solns, m;  m := 0;
         % Construct a list of the form ( (soln . mult) ... ):
         while zeros do <<
            f := diff(f,x);  m := m + 1;
            zeros := for each z in zeros join if
                  (if !*i_solve then zerop eval_uni_poly(f,z) % integer.
                  else null numr eval_uni_poly_q(f, z) ) % rational
               then {z} else << solns := (z . m) . solns; nil >>
            >>;
         solns := for each s in solns collect
            makeqn!-maybe(x, !*n2a car s, noeqs) . cdr s;
         multiplicities!* := nil;
         if mul eq 'separate then << % mul = separate (default)
            solns := for each s in solns collect
               <<multiplicities!* := cdr s . multiplicities!*; car s>>;
            multiplicities!* := reverse multiplicities!* >>
         else if mul eq 'expand then % cf. on multiplicities
            solns := for each s in solns join
               for i := 1 : cdr s collect car s
         else                           % mul = together
            solns := for each s in solns collect
               makelist {car s, cdr s};
         multiplicities!* := makelist multiplicities!*;
         return solns
      end
   end;

symbolic procedure makeqn!-maybe(x, s, noeqs);
   if noeqs then s else {'equal, x, s};

symbolic procedure !*n2a n;
   % Convert an integer or rational SQ to algebraic (prefix) form:
   if fixp n then n else !*q2a n;


% Modular arithmetic support, based on module alg/genmod
% ======================================================
% Avoids a few tests that are unnecessary for this application
% (I hope!).  Must be here because some definitions are smacros!

symbolic procedure mod!# n;
   % Reduce integer n modulo current!-modulus,
   % i.e. to range 0 <= n < current!-modulus:
   (if n_m < 0 then n_m + current!-modulus else n_m)
      where n_m = remainder(n, current!-modulus);

symbolic smacro procedure mod!+(a,b);
   % Sum of two positive modular numbers:
   general!-modular!-plus(a,b);

symbolic smacro procedure mod!*(a,b);
   % Product of two positive modular numbers:
   remainder(a*b, current!-modulus);

symbolic procedure mod!/(a,b);
   % Quotient of two positive modular numbers:
   mod!*(a, general!-reciprocal!-by!-gcd(current!-modulus,b,0,1));

symbolic smacro procedure mod!^(a,n);
   % Natural number power of a modular number:
   general!-modular!-expt(a,n);


% The main R_SOLVE procedure
% ==========================

fluid '(!*ezgcd !*gcd current!-modulus modulus!/2);

symbolic procedure r_solve f;
   % Returns the distinct RATIONAL zeros of a univariate integer
   % polynomial f by using quadratic Hensel lifting,
   % where f(x) = SUM_{i=0}^n a_i x^i = 0.
   % Assumes f is a standard form.
   % If !*i_solve then computes only INTEGER zeros.
   begin scalar x, f1, fm, f1m, l;
      scalar a_n, a_0, p, b, n, p!^n, !2!^r, !2!^r1, p2r, p2r1;
      scalar !*ezgcd, !*gcd, current!-modulus, modulus!/2;
      x := mvar f;  !*gcd := t;

      % Take squarefree factor (which removes numerical content)
      % and compute its derivative:
      f := (quotf!*(f, gcdf(f, diff(f,x))) where !*ezgcd = t);
      f1 := diff(f,x);
      tr_write "Squarefree poly: ", f;

      % Find abs leading coeff:
      a_n := abs lc f;
      % Find abs trailing coeff (not just constant term):
      a_0 := f;  while not atom a_0 and red a_0 do a_0 := red a_0;
      if not atom a_0 then a_0 := lc a_0;  a_0 := abs a_0;
      tr_write "Abs leading coeff: ", a_n,
            ",  abs trailing coeff: ", a_0;

      % Choose p as the smallest prime such that f
      % retains its degree and remains squarefree:
      % (Let the algebraic processor do the work here!)
      begin scalar old_mod, dmode!*;
         dmode!* := '!:mod!:;
         old_mod := setmod(p := 2);
         while zerop remainder(a_n, p) or
               gcdf(fm := resimpf f, resimpf f1) neq 1 do
            setmod(p := nextprime p);
         %% f1m := resimpf f1;
         setmod old_mod
      end;
      % Simplify back to integer domain mode:
      fm := resimpf fm;  %% f1m := resimpf f1m;
      tr_write "Prime modulus: ", p;
      tr_write "Poly mod ", p, ": ", fm;
      %% tr_write "Deriv mod ", p, ": ", f1m;

      % Find zero set of f mod p by simple exhaustive search:
      current!-modulus := p;
      for xx := 0 : p-1 do
         if zerop mod_eval_uni_poly(fm, xx) then l := xx . l;
      tr_write "Zero set mod ", p, ": ", l;
      if null l then return nil;

      % Determine modular prime power bound:
      b := 2 * if !*i_solve then a_0 else a_n*a_0;
      n := 1;  p!^n := p;
      while p!^n <= b do << n := n + 1;  p!^n := p * p!^n >>;
      tr_write "Modular prime power bound: ", n;

      % Lift to zeros of f mod p^N by quadratic Hensel lifting:
      !2!^r := 1;  p2r := p;            % r := 0
      while !2!^r < n do <<
         %% !2!^r1 := 2*!2!^r;  p2r1 := p2r^2;
         p2r1 := if (!2!^r1 := 2*!2!^r) < n then p2r^2
            else << !2!^r1 := n;  p!^n >>;
         tr_write "****************************************";
         tr_write "Lift modulo prime power ", p, "^", !2!^r1,
            " = ", p2r1;
         l := for each alpha_r in l collect
         begin integer fm, b_r, f1m, y;
            tr_write "Current zero mod ",p,"^",!2!^r,": ",alpha_r;
            current!-modulus := p2r1;
            fm := mod_eval_uni_poly(f, alpha_r);
            tr_write "fm: ", fm;
            if zerop fm then return alpha_r;
            b_r := fm/p2r;  % exact (non-modular) division
            current!-modulus := p2r;
            f1m := mod_eval_uni_poly(f1, alpha_r);
            tr_write "b_r: ", b_r, ",  f1m: ", f1m;
            % y := - b_r / f1m;
            % Must do this division mod p2r:
            y := mod!/(p2r - b_r, f1m);
            tr_write "y = -b_r/f1m mod ", p2r, ": ", y;
            return remainder(alpha_r + p2r*y, p2r1)
         end;
         tr_write "Zero set mod ", p, "^", !2!^r1, ": ", l;
         !2!^r := !2!^r1;  p2r := p2r1
      >>;

      return if !*i_solve then
         % Check for genuine (non-modular) integer roots:
         for each alpha in l join
            (if zerop eval_uni_poly(f, a) then {a})
               where a = balance_mod(alpha, p2r)
      else
         %% Reconstruct, using alpha = s/t mod p^N and t | a_n, so
         %% s/t = (alpha * a_n bal_mod p^N)/a_n canonicalized
         %% (because a_n = a_n bal_mod p^N by choice of N)
         %% and check genuine rational roots, because e.g.
         %% x^2 + 2 has zeros {1,2} mod 3, but no rational zeros.
         %% balance_mod is an odd function, i.e. commutes with sign
         %% change, so we can use a_n := abs lc f as computed above:
         for each alpha in l join <<
            if zerop alpha then alpha := nil ./ 1
            else <<
               alpha := balance_mod(remainder(alpha*a_n, p2r), p2r);
               alpha := ((alpha/g)./(a_n/g) where g = gcdn(alpha, a_n))
            >>;
            if null numr eval_uni_poly_q(f, alpha) then {alpha} >>
   end;


% Support procedures
% ==================

symbolic procedure resimpf f;
   % Resimplify a standard form.
   numr subf1(f, nil) where varstack!*=nil;

symbolic procedure balance_mod(n, m);
   % Shift n to balanced range modulo m
   % [-floor((m-1)/2), ceil((m-1)/2)]
   % (cf. modprep!: in POLY.RED)
   if n + n > m then n - m else n;

symbolic procedure eval_uni_poly(f, a);
   % Evaluate univariate integer polynomial
   % f(x) at x = a using Horner's rule.
   % f : standard form; a : integer
   % Returns an integer.
   if atom f then f else
   begin scalar r, d;
      r := lc f;  d := ldeg f;  f := red f;
      while not atom f do <<
         r := r*a^(d - ldeg f) + lc f;
         d := ldeg f;  f := red f
      >>;
      r := r*a^d;  if f then r := r + f;
      return r
   end;

symbolic procedure eval_uni_poly_q(f, a);
   % Evaluate univariate integer polynomial
   % f(x) at x = a using Horner's rule.
   % f : standard form; a : rational number as standard quotient
   % Returns a rational number as a standard quotient.
   if atom f then f ./ 1 else
   begin scalar r, d;
      r := lc f ./ 1;  d := ldeg f;  f := red f;
      while not atom f do <<
         r := addsq(multsq(r, exptsq(a,(d - ldeg f))), lc f ./ 1);
         d := ldeg f;  f := red f
      >>;
      return addsq(multsq(r, exptsq(a,d)), f ./ 1);
   end;

symbolic procedure mod_eval_uni_poly(f, a);
   % Evaluate modulo current!-modulus the univariate
   % integer polynomial f(x) at x = a using Horner's rule.
   % f : standard form; a : cardinal, 0 <= a < current!-modulus.
   if atom f then mod!# f else
   begin scalar r, d;
      r := mod!# lc f;  d := ldeg f;  f := red f;
      while not atom f do <<
         r := mod!+(mod!*(r, mod!^(a, d - ldeg f)), mod!# lc f);
         d := ldeg f;  f := red f
      >>;
      r := mod!*(r, mod!^(a,d));  if f then r := mod!+(r, mod!# f);
      return r
   end;


% Temporary patch to procedure in alg/genmod
% ==========================================
% Remove when fixed elsewhere.

% 0 is a valid modular number, but
% general!-modular!-expt(0,n) returns 1 instead of 0.

symbolic procedure general!-modular!-expt(a,n);
   % Computes a**n modulo current-modulus.  Uses Fermat's Little
   % Theorem where appropriate for a prime modulus.
    if a=0 then 0 else  % FJW
    if n=0 then 1
    else if n=1 then a
    else if n>=current!-modulus-1 and primep current!-modulus
     then general!-modular!-expt(a,remainder(n,current!-modulus-1))
    else begin scalar x;
     x:=general!-modular!-expt(a,n/2);
     x:=general!-modular!-times(x,x);
     if not evenp n then x:=general!-modular!-times(x,a);
     return x
    end;

endmodule;


end;
