1module polydiv; % Enhanced polynomial division. 2 3% Redistribution and use in source and binary forms, with or without 4% modification, are permitted provided that the following conditions are met: 5% 6% * Redistributions of source code must retain the relevant copyright 7% notice, this list of conditions and the following disclaimer. 8% * Redistributions in binary form must reproduce the above copyright 9% notice, this list of conditions and the following disclaimer in the 10% documentation and/or other materials provided with the distribution. 11% 12% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 13% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 14% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 15% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR 16% CONTRIBUTORS 17% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 18% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 19% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 20% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 21% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 22% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 23% POSSIBILITY OF SUCH DAMAGE. 24% 25 26% Francis Wright, 1995, 2020. 27 28% Define (or redefine) the following polynomial division operators: 29% divide, poly_quotient and remainder (mod), 30% pseudo_divide, pseudo_quotient and pseudo_remainder. 31% All accept a list as first argument/operand and map over it. 32 33% =================================================================== 34 35% Enhanced algebraic-mode operators for performing polynomial division 36% over the current coefficient domain, based on the operator REMAINDER 37% previously defined in "packages/poly/polrep.red" by 38% put('remainder,'polyfn,'remf); 39 40% divide(p,q) or p divide q returns an algebraic list {quotient, 41% remainder} of p divided by q as polynomials over the current domain. 42 43% poly_quotient(p,q) or p poly_quotient q returns only the quotient. 44% remainder(p,q) or p mod q returns only the remainder. 45 46% An optional third argument (for the prefix forms) specifies the 47% variable to treat as the leading variable for the (effectively 48% univariate) polynomial division. 49 50 51% Interface code 52% ============== 53 54% Regular division: 55% ---------------- 56 57put('divide, 'psopfn, 'poly!-divide); 58symbolic procedure poly!-divide u; 59 poly!-divide!*(u, nil, nil); 60 61remprop('remainder, 'polyfn); 62put('remainder, 'psopfn, 'poly!-remainder); 63put('mod, 'psopfn, 'poly!-remainder); 64symbolic procedure poly!-remainder u; 65 poly!-divide!*(u, 'remainder, nil); 66 67put('poly_quotient, 'psopfn, 'poly!-quotient); 68symbolic procedure poly!-quotient u; 69 poly!-divide!*(u, 'quotient, nil); 70 71infix divide, mod, poly_quotient; 72% Set a relatively low precedence (see preclis!*): 73precedence divide, freeof; % higher than freeof, lower than + 74precedence poly_quotient, divide; 75precedence mod, poly_quotient; 76% With integer arguments, mod calls evalmod, defined in "modular.red". 77 78 79% Pseudo-division: 80% --------------- 81 82put('pseudo_divide, 'psopfn, 'pseudo!-divide); 83symbolic procedure pseudo!-divide u; 84 poly!-divide!*(u, nil, t); 85 86put('pseudo_remainder, 'psopfn, 'pseudo!-remainder); 87symbolic procedure pseudo!-remainder u; 88 poly!-divide!*(u, 'remainder, t); 89 90put('pseudo_div, 'psopfn, 'pseudo!-quotient); % deprecated; to be removed from manual 91put('pseudo_quotient, 'psopfn, 'pseudo!-quotient); 92symbolic procedure pseudo!-quotient u; 93 poly!-divide!*(u, 'quotient, t); 94 95 96fluid '(kord!*); 97 98symbolic procedure poly!-divide!*(u, fn, pseudo); 99 if rlistp cadr u then typerr(cadr u, 'polynomial) 100 % Accept an algebraic list as first argument/operand and map over it: 101 else if rlistp car u then 102 'list . for each v in cdar u collect poly!-divide!*!*(v . cdr u, fn, pseudo) 103 else poly!-divide!*!*(u, fn, pseudo); 104 105symbolic procedure poly!-divide!*!*(u, fn, pseudo); % u = (p, q, x) 106 % Return the quotient and remainder of p (pseudo-)divided by q. 107 % If specified, x is made the leading variable before dividing, 108 % otherwise the first variable found is used. 109 begin scalar p, q, x, new_kord; 110 if null cdr u then rederr "Divisor required"; 111 if length u > 3 then 112 rederr "Division operators take 2 or 3 arguments."; 113 if null (q := !*a2f cadr u) then rederr "Zero divisor"; 114 p := !*a2f car u; 115 % Use integer modulus? 116 if fn eq 'remainder and fixp p and fixp q then 117 if q>0 then return evalmod u else typerr(q, "modulus"); 118 if cddr u and (x := !*a2k caddr u) and 119 not(kord!* and x eq car kord!*) then << 120 new_kord := t; updkorder x; 121 p := reorder p; q := reorder q 122 >> where kord!* = kord!*; % preserve environment 123 u := if pseudo then pseudo!-qremf(p, q, x) else qremf(p, q); 124 p := car u; q := cdr u; 125 if new_kord then << 126 if not(fn eq 'remainder) then p := reorder p; 127 if not(fn eq 'quotient) then q := reorder q 128 >>; 129 return 130 if fn eq 'remainder then mk!*sq (q ./ 1) 131 else if fn eq 'quotient then mk!*sq (p ./ 1) 132 else {'list, mk!*sq (p ./ 1), mk!*sq (q ./ 1)} 133 end; 134 135 136% Pseudo-division code 137% ==================== 138 139symbolic procedure pseudo!-qremf(u, v, var); 140 % Returns quotient and remainder of u pseudo-divided by v wrt var. 141 % u and v are standard forms, var is a kernel or nil. 142 % If var = nil then var := first kernel found. 143 % Internally, polynomials are represented as coeff lists wrt var, 144 % i.e. as lists of forms. 145 % (Knuth 1981, Seminumerical Algorithms, Algorithm R, page 407.) 146 begin scalar no_var, m, n, k, q0, q, car_v, car_u, vv; 147 no_var := null var; 148 m := if domainp u or (var and not(mvar u eq var)) then 0 149 else << if not var then var := mvar u; ldeg u >>; 150 n := if domainp v or (var and not(mvar v eq var)) then 0 151 else << if not var then var := mvar v; ldeg v >>; 152 153 %% The following special-case code for n = 0 and m < n is not 154 %% necessary, but could be a cheap efficiency measure. 155 %% if zerop n then return multf(exptf(v,m), u) . nil; 156 %% if minusp(k := m - n) then return nil . u; 157 158 u := if zerop m then {u} else coeffs u; 159 v := if zerop n then {v} else coeffs v; 160 if no_var and not(domainp_list v and domainp_list u) then 161 msgpri("Main division variable selected is", var, 162 nil, nil, nil); 163 k := m - n; car_v := car v; 164 while k >= 0 do << 165 %% Compute the quotient q EFFICIENTLY. 166 %% q0 = (q_0 ... q_k) without powers of v_n 167 q0 := (car_u := car u) . q0; 168 vv := cdr v; 169 u := for each c in cdr u collect << 170 c := multf(c, car_v); 171 if vv then << 172 c := subtrf(c, multf(car_u, car vv)); 173 vv := cdr vv 174 >>; 175 c 176 >>; 177 k := k - 1 178 >>; 179 if q0 then << 180 %% Reverse q0 and multiply in powers of v_n: 181 q := car q0 . nil; vv := 1; % v_n^0 182 while (q0 := cdr q0) do 183 q := multf(car q0, (vv := multf(vv, car_v))) . q 184 >>; 185 return coeffs!-to!-form(q, var) . coeffs!-to!-form(u, var) 186 end; 187 188symbolic procedure coeffs!-to!-form(coeff_list, var); 189 % Convert a coefficient list in DESCENDING power order to a 190 % standard form wrt the specified variable var: 191 coeff_list and 192 coeffs!-to!-form1(coeff_list, var, length coeff_list - 1); 193 194symbolic procedure coeffs!-to!-form1(coeff_list, var, d); 195 if d > 0 then 196 ( if car coeff_list then 197 ((var .^ d) .* (car coeff_list)) .+ reductum 198 else reductum ) 199 where reductum = 200 coeffs!-to!-form1(cdr coeff_list, var, d - 1) 201 else car coeff_list; 202 203symbolic procedure domainp_list coeff_list; 204 % Returns true if argument is a list of domain elements: 205 null coeff_list or 206 (domainp car coeff_list and domainp_list cdr coeff_list); 207 208endmodule; 209 210end; 211