Artifact 6d60efd890ffc4453b06378f82c68bb5e88b1aa6d9a7178710638d138f35f624:
- Executable file
r36/src/rsolve.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 16897) [annotate] [blame] [check-ins using] [more...]
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;