module polydiv; % Enhanced polynomial division.
% F.J.Wright@Maths.QMW.ac.uk, 6 Nov 1995.
% Defines (or redefines) the following polynomial division operators:
% divide, div and remainder (mod),
% pseudo_divide, pseudo_quotient (pseudo_div) and pseudo_remainder.
% However, for now, div has been commented out, since it conflicts with
% other packages (avector and fide).
% ===================================================================
% Enhanced algebraic-mode operators for performing polynomial division
% over the current coefficient domain, based on the operator REMAINDER
% currently defined in poly.red by put('remainder,'polyfn,'remf);
% divide(p,q) or p divide q returns an algebraic list {quotient,
% remainder} of p divided by q as polynomials over the current domain.
% div(p,q) or p div q returns only the quotient.
% remainder(p,q) or p mod q returns only the remainder.
% div and mod are the operator names used in Pascal for Euclidean
% (integer) division.
% An optional third argument (for the prefix forms) specifies the
% variable to treat as the leading variable for the (effectively
% univariate) polynomial division.
% Interface code
% ==============
% Regular division:
% ----------------
put('divide, 'psopfn, 'poly!-divide);
symbolic procedure poly!-divide u;
poly!-divide!*(u, nil, nil);
remprop('remainder, 'polyfn);
put('remainder, 'psopfn, 'poly!-remainder);
put('mod, 'psopfn, 'poly!-remainder); % name from Pascal
symbolic procedure poly!-remainder u;
poly!-divide!*(u, 'remainder, nil);
% put('div, 'psopfn, 'poly!-quotient); % name from Pascal
symbolic procedure poly!-quotient u;
poly!-divide!*(u, 'quotient, nil);
infix divide, mod;
% infix div;
% Set a relatively low precedence:
precedence divide, freeof; % higher than freeof, lower than +
% precedence div, divide;
% precedence mod, div;
% Pseudo-division:
% ---------------
put('pseudo_divide, 'psopfn, 'pseudo!-divide);
symbolic procedure pseudo!-divide u;
poly!-divide!*(u, nil, t);
put('pseudo_remainder, 'psopfn, 'pseudo!-remainder);
symbolic procedure pseudo!-remainder u;
poly!-divide!*(u, 'remainder, t);
put('pseudo_div, 'psopfn, 'pseudo!-quotient);
put('pseudo_quotient, 'psopfn, 'pseudo!-quotient);
symbolic procedure pseudo!-quotient u;
poly!-divide!*(u, 'quotient, t);
fluid '(kord!*);
symbolic procedure poly!-divide!*(u, fn, pseudo); % u = (p, q, x)
% Returns the quotient and remainder of p (pseudo-)divided by q.
% If specified, x is made the leading variable before dividing,
% otherwise the first variable found is used.
begin scalar p, q, x, new_kord;
if null cdr u then RedErr "Divisor required";
if length u > 3 then
RedErr "Division operators take 2 or 3 arguments.";
if null (q := !*a2f cadr u) then RedErr "Zero divisor";
p := !*a2f car u;
if cddr u and (x := !*a2k caddr u) and
not(kord!* and x eq car kord!*) then <<
new_kord := t; updkorder x;
p := reorder p; q := reorder q
>> where kord!* = kord!*; % preserve environment
u := if pseudo then pseudo!-qremf(p, q, x) else qremf(p, q);
p := car u; q := cdr u;
if new_kord then <<
if not(fn eq 'remainder) then p := reorder p;
if not(fn eq 'quotient) then q := reorder q
>>;
return
if fn eq 'remainder then mk!*sq (q ./ 1)
else if fn eq 'quotient then mk!*sq (p ./ 1)
else {'list, mk!*sq (p ./ 1), mk!*sq (q ./ 1)}
end;
% Pseudo-division code
% ====================
symbolic procedure pseudo!-qremf(u, v, var);
% Returns quotient and remainder of u pseudo-divided by v wrt var.
% u and v are standard forms, var is a kernel or nil.
% If var = nil then var := first kernel found.
% Internally, polynomials are represented as coeff lists wrt var,
% i.e. as lists of forms.
% (Knuth 1981, Seminumerical Algorithms, Algorithm R, page 407.)
begin scalar no_var, m, n, k, q0, q, car_v, car_u, vv;
no_var := null var;
m := if domainp u or (var and not(mvar u eq var)) then 0
else << if not var then var := mvar u; ldeg u >>;
n := if domainp v or (var and not(mvar v eq var)) then 0
else << if not var then var := mvar v; ldeg v >>;
%% The following special-case code for n = 0 and m < n is not
%% necessary, but could be a cheap efficiency measure.
%% if zerop n then return multf(exptf(v,m), u) . nil;
%% if minusp(k := m - n) then return nil . u;
u := if zerop m then {u} else coeffs u;
v := if zerop n then {v} else coeffs v;
if no_var and not(domainp_list v and domainp_list u) then
msgpri("Main division variable selected is", var,
nil, nil, nil);
k := m - n; car_v := car v;
while k >= 0 do <<
%% Compute the quotient q EFFICIENTLY.
%% q0 = (q_0 ... q_k) without powers of v_n
q0 := (car_u := car u) . q0;
vv := cdr v;
u := for each c in cdr u collect <<
c := multf(c, car_v);
if vv then <<
c := subtrf(c, multf(car_u, car vv));
vv := cdr vv
>>;
c
>>;
k := k - 1
>>;
if q0 then <<
%% Reverse q0 and multiply in powers of v_n:
q := car q0 . nil; vv := 1; % v_n^0
while (q0 := cdr q0) do
q := multf(car q0, (vv := multf(vv, car_v))) . q
>>;
return coeffs!-to!-form(q, var) . coeffs!-to!-form(u, var)
end;
symbolic procedure coeffs!-to!-form(coeff_list, var);
% Convert a coefficient list in DESCENDING power order to a
% standard form wrt the specified variable var:
coeff_list and
coeffs!-to!-form1(coeff_list, var, length coeff_list - 1);
symbolic procedure coeffs!-to!-form1(coeff_list, var, d);
if d > 0 then
( if car coeff_list then
((var .^ d) .* (car coeff_list)) .+ reductum
else reductum )
where reductum =
coeffs!-to!-form1(cdr coeff_list, var, d - 1)
else car coeff_list;
symbolic procedure domainp_list coeff_list;
% Returns true if argument is a list of domain elements:
null coeff_list or
(domainp car coeff_list and domainp_list cdr coeff_list);
endmodule;
end;