Polynomial long division
Let us suppose a polynomial is represented by a vector, (i.e., an ordered collection of coefficients) so that the th element keeps the coefficient of , and the multiplication by a monomial is a shift of the vector's elements "towards right" (injecting zeros from left) followed by a multiplication of each element by the coefficient of the monomial.
You are encouraged to solve this task according to the task description, using any language you may know.
This page uses content from Wikipedia. The original article was at Polynomial long division. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance) |
- In algebra, polynomial long division is an algorithm for dividing a polynomial by another polynomial of the same or lower degree.
Then a pseudocode for the polynomial long division using the conventions described above could be:
degree(P): return the index of the last non-zero element of P; if all elements are 0, return -∞ polynomial_long_division(N, D) returns (q, r): // N, D, q, r are vectors if degree(D) < 0 then error if degree(N) ≥ degree(D) then q ← 0 while degree(N) ≥ degree(D) d ← D shifted right by (degree(N) - degree(D)) q(degree(N) - degree(D)) ← N(degree(N)) / d(degree(d)) // by construction, degree(d) = degree(N) of course d ← d * q(degree(N) - degree(D)) N ← N - d endwhile r ← N else q ← 0 r ← N endif return (q, r)
Note: vector * scalar
multiplies each element of the vector by the scalar; vectorA - vectorB
subtracts each element of the vectorB from the element of the vectorA with "the same index". The vectors in the pseudocode are zero-based.
- Error handling (for allocations or for wrong inputs) is not mandatory.
- Conventions can be different; in particular, note that if the first coefficient in the vector is the highest power of x for the polynomial represented by the vector, then the algorithm becomes simpler.
Example for clarification
This example is from Wikipedia, but changed to show how the given pseudocode works.
0 1 2 3 ---------------------- N: -42 0 -12 1 degree = 3 D: -3 1 0 0 degree = 1 d(N) - d(D) = 2, so let's shift D towards right by 2: N: -42 0 -12 1 d: 0 0 -3 1 N(3)/d(3) = 1, so d is unchanged. Now remember that "shifting by 2" is like multiplying by x2, and the final multiplication (here by 1) is the coefficient of this monomial. Let's store this into q: 0 1 2 --------------- q: 0 0 1 now compute N - d, and let it be the "new" N, and let's loop N: -42 0 -9 0 degree = 2 D: -3 1 0 0 degree = 1 d(N) - d(D) = 1, right shift D by 1 and let it be d N: -42 0 -9 0 d: 0 -3 1 0 * -9/1 = -9 q: 0 -9 1 d: 0 27 -9 0 N ← N - d N: -42 -27 0 0 degree = 1 D: -3 1 0 0 degree = 1 looping again... d(N)-d(D)=0, so no shift is needed; we multiply D by -27 (= -27/1) storing the result in d, then q: -27 -9 1 and N: -42 -27 0 0 - d: 81 -27 0 0 = N: -123 0 0 0 (last N) d(N) < d(D), so now r ← N, and the result is: 0 1 2 ------------- q: -27 -9 1 → x2 - 9x - 27 r: -123 0 0 → -123
C
<lang c>#include <stdio.h>
- include <stdlib.h>
- include <stdarg.h>
- include <assert.h>
- include <gsl/gsl_vector.h>
- define MAX(A,B) (((A)>(B))?(A):(B))
void reoshift(gsl_vector *v, int h) {
if ( h > 0 ) { gsl_vector *temp = gsl_vector_alloc(v->size); gsl_vector_view p = gsl_vector_subvector(v, 0, v->size - h); gsl_vector_view p1 = gsl_vector_subvector(temp, h, v->size - h); gsl_vector_memcpy(&p1.vector, &p.vector); p = gsl_vector_subvector(temp, 0, h); gsl_vector_set_zero(&p.vector); gsl_vector_memcpy(v, temp); gsl_vector_free(temp); }
}
gsl_vector *poly_long_div(gsl_vector *n, gsl_vector *d, gsl_vector **r) {
gsl_vector *nt = NULL, *dt = NULL, *rt = NULL, *d2 = NULL, *q = NULL; int gn, gt, gd;
if ( (n->size >= d->size) && (d->size > 0) && (n->size > 0) ) { nt = gsl_vector_alloc(n->size); assert(nt != NULL); dt = gsl_vector_alloc(n->size); assert(dt != NULL); rt = gsl_vector_alloc(n->size); assert(rt != NULL); d2 = gsl_vector_alloc(n->size); assert(d2 != NULL); gsl_vector_memcpy(nt, n); gsl_vector_set_zero(dt); gsl_vector_set_zero(rt); gsl_vector_view p = gsl_vector_subvector(dt, 0, d->size); gsl_vector_memcpy(&p.vector, d); gsl_vector_memcpy(d2, dt); gn = n->size - 1; gd = d->size - 1; gt = 0;
while( gsl_vector_get(d, gd) == 0 ) gd--; while ( gn >= gd ) { reoshift(dt, gn-gd); double v = gsl_vector_get(nt, gn)/gsl_vector_get(dt, gn); gsl_vector_set(rt, gn-gd, v); gsl_vector_scale(dt, v); gsl_vector_sub(nt, dt); gt = MAX(gt, gn-gd); while( (gn>=0) && (gsl_vector_get(nt, gn) == 0.0) ) gn--; gsl_vector_memcpy(dt, d2); }
q = gsl_vector_alloc(gt+1); assert(q != NULL); p = gsl_vector_subvector(rt, 0, gt+1); gsl_vector_memcpy(q, &p.vector); if ( r != NULL ) { if ( (gn+1) > 0 ) {
*r = gsl_vector_alloc(gn+1); assert( *r != NULL ); p = gsl_vector_subvector(nt, 0, gn+1); gsl_vector_memcpy(*r, &p.vector);
} else {
*r = gsl_vector_alloc(1); assert( *r != NULL ); gsl_vector_set_zero(*r);
} } gsl_vector_free(nt); gsl_vector_free(dt); gsl_vector_free(rt); gsl_vector_free(d2); return q; } else { q = gsl_vector_alloc(1); assert( q != NULL ); gsl_vector_set_zero(q); if ( r != NULL ) { *r = gsl_vector_alloc(n->size); assert( *r != NULL ); gsl_vector_memcpy(*r, n); } return q; }
}
void poly_print(gsl_vector *p) {
int i; for(i=p->size-1; i >= 0; i--) { if ( i > 0 ) printf("%lfx^%d + ",
gsl_vector_get(p, i), i);
else printf("%lf\n", gsl_vector_get(p, i)); }
}
gsl_vector *create_poly(int d, ...) {
va_list al; int i; gsl_vector *r = NULL;
va_start(al, d); r = gsl_vector_alloc(d); assert( r != NULL ); for(i=0; i < d; i++) gsl_vector_set(r, i, va_arg(al, double));
return r;
}</lang>
<lang c>int main() {
int i; gsl_vector *q, *r; gsl_vector *nv, *dv; //nv = create_poly(4, -42., 0., -12., 1.); //dv = create_poly(2, -3., 1.); //nv = create_poly(3, 2., 3., 1.); //dv = create_poly(2, 1., 1.); nv = create_poly(4, -42., 0., -12., 1.); dv = create_poly(3, -3., 1., 1.);
q = poly_long_div(nv, dv, &r);
poly_print(q); poly_print(r);
gsl_vector_free(q); gsl_vector_free(r);
return 0;
}</lang>
E
This program has some unnecessary features contributing to its length:
- It creates polynomial objects rather than performing its operations directly on arrays.
- It includes code for printing polynomials nicely.
- It prints the intermediate steps of the division.
pragma.syntax("0.9") pragma.enable("accumulator") def superscript(x, out) { if (x >= 10) { superscript(x // 10) } out.print("⁰¹²³⁴⁵⁶⁷⁸⁹"[x %% 10]) } def makePolynomial(initCoeffs :List) { def degree := { var i := initCoeffs.size() - 1 while (i >= 0 && initCoeffs[i] <=> 0) { i -= 1 } if (i < 0) { -Infinity } else { i } } def coeffs := initCoeffs(0, degree + 1) def polynomial { /** Print the polynomial (not necessary for the task) */ to __printOn(out) { out.print("(λx. ") var first := true for i in (0..!(coeffs.size())).descending() { def coeff := coeffs[i] if (coeff <=> 0) { continue } if (!first) { out.print(" ") } if (coeff <=> 1 && !(i <=> 0)) { # no coefficient written if it's 1 and not the constant term } else if (first) { out.print(coeff) } else if (coeff > 0) { out.print("+ ", coeff) } else { out.print("- ", -coeff) } if (i <=> 0) { # no x if it's the constant term } else if (i <=> 1) { out.print("x") } else { out.print("x"); superscript(i, out) } first := false } out.print(")") } /** Evaluate the polynomial (not necessary for the task) */ to run(x) { return accum 0 for i => c in coeffs { _ + c * x**i } } to degree() { return degree } to coeffs() { return coeffs } to highestCoeff() { return coeffs[degree] } /** Could support another polynomial, but not part of this task. Computes this * x**power. */ to timesXToThe(power) { return makePolynomial([0] * power + coeffs) } /** Multiply (by a scalar only). */ to multiply(scalar) { return makePolynomial(accum [] for x in coeffs { _.with(x * scalar) }) } /** Subtract (by another polynomial only). */ to subtract(other) { def oc := other.coeffs() :List return makePolynomial(accum [] for i in 0..(coeffs.size().max(oc.size())) { _.with(coeffs.fetch(i, fn{0}) - oc.fetch(i, fn{0})) }) } /** Polynomial long division. */ to quotRem(denominator, trace) { var numerator := polynomial require(denominator.degree() >= 0) if (numerator.degree() < denominator.degree()) { return [makePolynomial([]), denominator] } else { var quotientCoeffs := [0] * (numerator.degree() - denominator.degree()) while (numerator.degree() >= denominator.degree()) { trace.print(" ", numerator, "\n") def qCoeff := numerator.highestCoeff() / denominator.highestCoeff() def qPower := numerator.degree() - denominator.degree() quotientCoeffs with= (qPower, qCoeff) def d := denominator.timesXToThe(qPower) * qCoeff trace.print("- ", d, " (= ", denominator, " * ", qCoeff, "x"); superscript(qPower, trace); trace.print(")\n") numerator -= d trace.print(" -------------------------- (Quotient so far: ", makePolynomial(quotientCoeffs), ")\n") } return [makePolynomial(quotientCoeffs), numerator] } } } return polynomial }
<lang e>def n := makePolynomial([-42, 0, -12, 1]) def d := makePolynomial([-3, 1]) println("Numerator: ", n) println("Denominator: ", d) def [q, r] := n.quotRem(d, stdout) println("Quotient: ", q) println("Remainder: ", r) </lang>
Output:
Numerator: (λx. x³ - 12x² - 42) Denominator: (λx. x - 3) (λx. x³ - 12x² - 42) - (λx. x³ - 3.0x²) (= (λx. x - 3) * 1.0x²) -------------------------- (Quotient so far: (λx. x²)) (λx. -9.0x² - 42.0) - (λx. -9.0x² + 27.0x) (= (λx. x - 3) * -9.0x¹) -------------------------- (Quotient so far: (λx. x² - 9.0x)) (λx. -27.0x - 42.0) - (λx. -27.0x + 81.0) (= (λx. x - 3) * -27.0x⁰) -------------------------- (Quotient so far: (λx. x² - 9.0x - 27.0)) Quotient: (λx. x² - 9.0x - 27.0) Remainder: (λx. -123.0)
Fortran
<lang fortran>module Polynom
implicit none
contains
subroutine poly_long_div(n, d, q, r) real, dimension(:), intent(in) :: n, d real, dimension(:), intent(out), allocatable :: q real, dimension(:), intent(out), allocatable, optional :: r
real, dimension(:), allocatable :: nt, dt, rt integer :: gn, gt, gd
if ( (size(n) >= size(d)) .and. (size(d) > 0) .and. (size(n) > 0) ) then allocate(nt(size(n)), dt(size(n)), rt(size(n)))
nt = n dt = 0 dt(1:size(d)) = d rt = 0 gn = size(n)-1 gd = size(d)-1 gt = 0
do while ( d(gd+1) == 0 ) gd = gd - 1 end do
do while( gn >= gd ) dt = eoshift(dt, -(gn-gd)) rt(gn-gd+1) = nt(gn+1) / dt(gn+1) nt = nt - dt * rt(gn-gd+1) gt = max(gt, gn-gd) do gn = gn - 1 if ( nt(gn+1) /= 0 ) exit end do dt = 0 dt(1:size(d)) = d end do
allocate(q(gt+1)) q = rt(1:gt+1) if ( present(r) ) then if ( (gn+1) > 0 ) then allocate(r(gn+1)) r = nt(1:gn+1) else allocate(r(1)) r = 0.0 end if end if deallocate(nt, dt, rt) else allocate(q(1)) q = 0 if ( present(r) ) then allocate(r(size(n))) r = n end if end if
end subroutine poly_long_div
subroutine poly_print(p) real, dimension(:), intent(in) :: p
integer :: i
do i = size(p), 1, -1 if ( i > 1 ) then write(*, '(F0.2,"x^",I0," + ")', advance="no") p(i), i-1 else write(*, '(F0.2)') p(i) end if end do
end subroutine poly_print
end module Polynom</lang>
<lang fortran>program PolyDivTest
use Polynom implicit none
real, dimension(:), allocatable :: q real, dimension(:), allocatable :: r
!! three tests from Wikipedia, plus an extra !call poly_long_div( (/ -3., 1. /), (/ -42., 0.0, -12., 1. /), q, r) call poly_long_div( (/ -42., 0.0, -12., 1. /), (/ -3., 1. /), q, r) !call poly_long_div( (/ -42., 0.0, -12., 1. /), (/ -3., 1., 1. /), q, r) !call poly_long_div( (/ 2., 3., 1. /), (/ 1., 1. /), q, r)
call poly_print(q) call poly_print(r) deallocate(q, r)
end program PolyDivTest</lang>
OCaml
First define some utility operations on polynomials as lists (with highest power coefficient first). <lang ocaml> let rec shift n l = if n <= 0 then l else shift (pred n) (l @ [0.0]) let rec pad n l = if n <= 0 then l else pad (pred n) (0.0 :: l) let rec scale k l = List.map (( *.) k) l let rec norm = function | 0.0 :: tl -> norm tl | x -> x let deg l = List.length (norm l) - 1
let zip op p q =
let gap = (List.length p) - (List.length q) in if gap = 0 then List.map2 op p q else if gap > 0 then List.map2 op p (pad gap q) else List.map2 op (pad (-gap) p) q
let sub = zip (-.) let add = zip (+.) </lang> Then the main polynomial division function <lang ocaml> let polydiv f g =
let rec aux f s q = let gap = (deg f) - (deg s) in if gap < 0 then (q, f) else let k = (List.hd f) /. (List.hd s) in let q' = add q (scale k (shift gap [1.0])) in let f' = norm (List.tl (sub f (scale k (shift gap s)))) in aux f' s q' in aux (norm f) (norm g) []
</lang> For output we need a pretty-printing function <lang ocaml> let str_poly l =
let term v p = match (v, p) with | ( _, 0) -> string_of_float v | (1.0, 1) -> "x" | ( _, 1) -> (string_of_float v) ^ "*x" | (1.0, _) -> "x^" ^ (string_of_int p) | _ -> (string_of_float v) ^ "*x^" ^ (string_of_int p) in let rec terms = function | [] -> [] | h :: t -> if h = 0.0 then (terms t) else (term h (List.length t)) :: (terms t) in String.concat " + " (terms l)
</lang> and the test driver <lang ocaml> let _ =
let f = [1.0; -4.0; 6.0; 5.0; 3.0] and g = [1.0; 2.0; 1.0] in let q, r = polydiv f g in begin Printf.printf "f = %s\ng = %s\n\n" (str_poly f) (str_poly g); Printf.printf "q = %s\nr = %s\n" (str_poly q) (str_poly r) end
</lang> which gives the output:
f = x^4 + -4.*x^3 + 6.*x^2 + 5.*x + 3. g = x^2 + 2.*x + 1. q = x^2 + -6.*x + 17. r = -23.*x + -14.
Octave
Octave has already facilities to divide two polynomials (deconv(n,d)
); and the reason to adopt the convention of keeping the highest power coefficient first, is to make the code compatible with builtin functions: we can use polyout to output the result.
<lang octave>function [q, r] = poly_long_div(n, d)
gd = length(d); pv = zeros(1, length(n)); pv(1:gd) = d; if ( length(n) >= gd ) q = []; while ( length(n) >= gd ) q = [q, n(1)/pv(1)]; n = n - pv .* (n(1)/pv(1)); n = shift(n, -1); % tn = n(1:length(n)-1); % eat the higher power term n = tn; % tp = pv(1:length(pv)-1); pv = tp; % make pv the same length of n endwhile r = n; else q = [0]; r = n; endif
endfunction
[q, r] = poly_long_div([1,-12,0,-42], [1,-3]); polyout(q, 'x'); polyout(r, 'x'); disp(""); [q, r] = poly_long_div([1,-12,0,-42], [1,1,-3]); polyout(q, 'x'); polyout(r, 'x'); disp(""); [q, r] = poly_long_div([1,3,2], [1,1]); polyout(q, 'x'); polyout(r, 'x'); disp(""); [q, r] = poly_long_div([1,3], [1,-12,0,-42]); polyout(q, 'x'); polyout(r, 'x');</lang>
Perl
This solution keeps the highest power coefficient first, like OCaml solution and Octave solution.
<lang perl>use strict; use List::Util qw(min);
sub poly_long_div {
my ($rn, $rd) = @_; my @n = @$rn; my $gd = scalar(@$rd); if ( scalar(@n) >= $gd ) {
my @q = (); while ( scalar(@n) >= $gd ) { my $piv = $n[0]/$rd->[0]; push @q, $piv; $n[$_] -= $rd->[$_] * $piv foreach ( 0 .. min(scalar(@n), $gd)-1 ); shift @n; } return ( \@q, \@n );
} else {
return ( [0], $rn );
}
}</lang>
<lang perl>sub poly_print {
my @c = @_; my $l = scalar(@c); for(my $i=0; $i < $l; $i++) {
print $c[$i]; print "x^" . ($l-$i-1) . " + " if ($i < ($l-1));
} print "\n";
}</lang>
<lang perl>my ($q, $r);
($q, $r) = poly_long_div([1, -12, 0, -42], [1, -3]); poly_print(@$q); poly_print(@$r); print "\n"; ($q, $r) = poly_long_div([1,-12,0,-42], [1,1,-3]); poly_print(@$q); poly_print(@$r); print "\n"; ($q, $r) = poly_long_div([1,3,2], [1,1]); poly_print(@$q); poly_print(@$r); print "\n";
- the example from the OCaml solution
($q, $r) = poly_long_div([1,-4,6,5,3], [1,2,1]); poly_print(@$q); poly_print(@$r);</lang>
Python
<lang python># -*- coding: utf-8 -*-
from itertools import izip
def degree(poly):
deg = len(poly) - 1 while deg >= 0 and poly[deg] == 0: del poly[deg] # normalize deg -= 1 return deg
def poly_div(N, D):
dD = degree(D) dN = degree(N) if dD < 0: raise ZeroDivisionError if dN >= dD: q = [0] * dN while dN >= dD: d = [0]*(dN - dD) + D mult = q[dN - dD] = N[-1] / float(d[-1]) d = [coeff*mult for coeff in d] N = [coeffN - coeffd for coeffN, coeffd in izip(N, d)] dN = degree(N) r = N else: q = [0] r = N return q, r
if __name__ == '__main__':
print "POLYNOMIAL LONG DIVISION" N = [-42, 0, -12, 1] D = [-3, 1, 0, 0] print " %s / %s =" % (N,D), print " %s remainder %s" % poly_div(N, D)</lang>
Sample output:
POLYNOMIAL LONG DIVISION [-42, 0, -12, 1] / [-3, 1, 0, 0] = [-27.0, -9.0, 1.0] remainder [-123.0]
R
<lang R>polylongdiv <- function(n,d) {
gd <- length(d) pv <- vector("numeric", length(n)) pv[1:gd] <- d if ( length(n) >= gd ) { q <- c() while ( length(n) >= gd ) { q <- c(q, n[1]/pv[1]) n <- n - pv * (n[1]/pv[1]) n <- n[2:length(n)] pv <- pv[1:(length(pv)-1)] } list(q=q, r=n) } else { list(q=c(0), r=n) }
}
- an utility function to print polynomial
print.polynomial <- function(p) {
i <- length(p)-1 for(a in p) { if ( i == 0 ) { cat(a, "\n") } else { cat(a, "x^", i, " + ", sep="") } i <- i - 1 }
}
r <- polylongdiv(c(1,-12,0,-42), c(1,-3)) print.polynomial(r$q) print.polynomial(r$r)</lang>
Ruby
Implementing the algorithm given in the task description: <lang ruby>def polynomial_long_division(numerator, denominator)
dd = degree(denominator) raise ArgumentError, "denominator is zero" if dd < 0 if dd == 0 return [multiply(numerator, 1.0/denominator[0]), [0]*numerator.length] end q = Array.new(numerator.length) {0} while (dn = degree(numerator)) >= dd d = shift_right(denominator, dn - dd) q[dn-dd] = numerator[dn] / d[degree(d)] d = multiply(d, q[dn-dd]) numerator = subtract(numerator, d) end [q, numerator]
end
def degree(ary)
idx = ary.reverse.find_index {|x| x.nonzero?} idx.nil? ? -1 : ary.length - 1 - idx
end
def shift_right(ary, n)
[0]*n + ary[0, ary.length - n]
end
def subtract(a1, a2)
a1.zip(a2).collect {|v1,v2| v1 - v2}
end
def multiply(ary, num)
ary.collect {|x| x * num}
end
f = [-42, 0, -12, 1] g = [-3, 1, 0, 0] q, r = polynomial_long_division(f, g) p [f, g, q, r]
- => [[-42, 0, -12, 1], [-3, 1, 0, 0], [-27, -9, 1, 0, 0], [-123, 0, 0, 0, 0]]
g = [-3, 1, 1, 0] q, r = polynomial_long_division(f, g) p [f, g, q, r]
- => [[-42, 0, -12, 1], [-3, 1, 0, 0], [-13, 1, 0, 0, 0], [-81, 16, 0, 0, 0]]</lang>
Implementing the algorithms on the wikipedia page -- uglier code but nicer user interface <lang ruby>def polynomial_division(f, g)
if g.length == 0 or (g.length == 1 and g[0] == 0) raise ArgumentError, "denominator is zero" elsif g.length == 1 [f.collect {|x| Float(x)/g[0]}, [0]] elsif g.length == 2 synthetic_division(f, g) else higher_degree_synthetic_division(f, g) end
end
def synthetic_division(f, g)
board = [f] << Array.new(f.length) << Array.new(f.length) board[2][0] = board[0][0] 1.upto(f.length - 1).each do |i| board[1][i] = board[2][i-1] * -g[1] board[2][i] = board[0][i] + board[1][i] end [board[2][0..-2], [board[2][-1]]]
end
- an ugly mess of array index arithmetic
- http://en.wikipedia.org/w/index.php?title=Polynomial_long_division#Higher_degree_synthetic_division
def higher_degree_synthetic_division(f, g)
# [use] the negative coefficients of the denominator following the leading term lhs = g[1..-1].collect {|x| -x} board = [f] q = [] 1.upto(f.length - lhs.length).each do |i| n = 2*i - 1 # underline the leading coefficient of the right-hand side, multiply it by # the left-hand coefficients and write the products beneath the next columns # on the right. q << board[n-1][i-1] board << Array.new(f.length).fill(0, i) # row n (lhs.length).times do |j| board[n][i+j] = q[-1]*lhs[j] end # perform an addition board << Array.new(f.length).fill(0, i) # row n+1 (lhs.length + 1).times do |j| board[n+1][i+j] = board[n-1][i+j] + board[n][i+j] if i+j < f.length end end # the remaining numbers in the bottom row correspond to the coefficients of the remainder r = board[-1].find_all {|x| not x.nil?} q = [0] if q.empty? [q, r]
end
f = [1, -12, 0, -42] g = [1, -3] q, r = polynomial_division(f, g) p [f, g, q, r]
- => [[1, -12, 0, -42], [1, -3], [1, -9, -27], [-123]]
g = [1, 1, -3] q, r = polynomial_division(f, g) p [f, g, q, r]
- => [[1, -12, 0, -42], [1, 1, -3], [1, -13], [16, -81]]</lang>
Best of both worlds:
<lang ruby>def tcl_polynomial_division(f, g)
if g.length == 0 or (g.length == 1 and g[0] == 0) raise ArgumentError, "denominator is zero" end return [[0], f] if f.length < g.length q = [] n, d = f.dup, g while n.length >= d.length q << Float(n[0]) / d[0] n[0, d.length].zip(d).each_with_index do |pair, i| n[i] = Float(pair[0]) - q[-1] * pair[1] end n.shift end q = [0] if q.empty? n = [0] if n.empty? [q, n]
end
f = [1, -12, 0, -42] g = [1, -3] q, r = polynomial_division(f, g) p [f, g, q, r]
- => [[1, -12, 0, -42], [1, -3], [1, -9, -27], [-123]]
g = [1, 1, -3] q, r = polynomial_division(f, g) p [f, g, q, r]
- => [[1, -12, 0, -42], [1, 1, -3], [1, -13], [16, -81]]</lang>
Smalltalk
<lang smalltalk>Object subclass: Polynomial [
|coeffs| Polynomial class >> new [ ^ super basicNew init ] init [ coeffs := OrderedCollection new. ^ self ] Polynomial class >> newWithCoefficients: coefficients [ |r| r := super basicNew. ^ r initWithCoefficients: coefficients ] initWithCoefficients: coefficients [ coeffs := coefficients asOrderedCollection. ^ self ] / denominator [ |n q| n := self deepCopy. self >= denominator ifTrue: [ q := Polynomial new. [ n >= denominator ] whileTrue: [ |piv| piv := (n coeff: 0) / (denominator coeff: 0).
q addCoefficient: piv. n := n - (denominator * piv). n clean
]. ^ { q . (n degree) > 0 ifTrue: [ n ] ifFalse: [ n addCoefficient: 0. n ] } ] ifFalse: [ ^ { Polynomial newWithCoefficients: #( 0 ) . self deepCopy } ] ] * constant [ |r| r := self deepCopy. 1 to: (coeffs size) do: [ :i | r at: i put: ((r at: i) * constant) ]. ^ r ] at: index [ ^ coeffs at: index ] at: index put: obj [ ^ coeffs at: index put: obj ] >= anotherPoly [ ^ (self degree) >= (anotherPoly degree) ] degree [ ^ coeffs size ] - anotherPoly [ "This is not a real subtraction between Polynomial: it is an internal method ..." |a| a := self deepCopy. 1 to: ( (coeffs size) min: (anotherPoly degree) ) do: [ :i | a at: i put: ( (a at: i) - (anotherPoly at: i) ) ]. ^ a ] coeff: index [ ^ coeffs at: (index + 1) ] addCoefficient: coeff [ coeffs add: coeff ] clean [ [ (coeffs size) > 0 ifTrue: [ (coeffs at: 1) = 0 ] ifFalse: [ false ] ] whileTrue: [ coeffs removeFirst ]. ] display [ 1 to: (coeffs size) do: [ :i | (coeffs at: i) display. i < (coeffs size) ifTrue: [ ('x^%1 + ' % {(coeffs size) - i} ) display ] ] ] displayNl [ self display. Character nl display ]
].</lang>
<lang smalltalk>|res| res := OrderedCollection new.
res add: ((Polynomial newWithCoefficients: #( 1 -12 0 -42) ) /
(Polynomial newWithCoefficients: #( 1 -3 ) )) ; add: ((Polynomial newWithCoefficients: #( 1 -12 0 -42) ) / (Polynomial newWithCoefficients: #( 1 1 -3 ) )).
res do: [ :o |
(o at: 1) display. ' with rest: ' display. (o at: 2) displayNl
] </lang>
Tcl
<lang tcl># poldiv - Divide two polynomials n and d.
- Result is a list of two polynomials, q and r, where n = qd + r
- and the degree of r is less than the degree of b.
- Polynomials are represented as lists, where element 0 is the
- x**0 coefficient, element 1 is the x**1 coefficient, and so on.
proc poldiv {a b} {
# Toss out leading zero coefficients efficiently while {[lindex $a end] == 0} {set a [lrange $a[set a {}] 0 end-1]} while {[lindex $b end] == 0} {set b [lrange $b[set b {}] 0 end-1]} if {[llength $a] < [llength $b]} { return [list 0 $a] }
# Rearrange the terms to put highest powers first set n [lreverse $a] set d [lreverse $b]
# Carry out classical long division, accumulating quotient coefficients # in q, and replacing n with the remainder. set q {} while {[llength $n] >= [llength $d]} { set qd [expr {[lindex $n 0] / [lindex $d 0]}] set i 0 foreach nd [lrange $n 0 [expr {[llength $d] - 1}]] dd $d { lset n $i [expr {$nd - $qd * $dd}] incr i } lappend q $qd set n [lrange $n 1 end] }
# Return quotient and remainder, constant term first return [list [lreverse $q] [lreverse $n]]
}
- Demonstration
lassign [poldiv {-42. 0. -12. 1.} {-3. 1. 0. 0.}] Q R puts [list Q = $Q] puts [list R = $R]</lang>
Ursala
The input is a pair of lists of coefficients in order of increasing degree. Trailing zeros can be omitted. The output is a pair of lists (q,r), the quotient and remainder polynomial coefficients. This is a straightforward implementation of the algorithm in terms of list operations (fold, zip, map, distribute, etc.) instead of array indexing, hence not unnecessarily verbose. <lang>#import std
- import flo
polydiv =
zeroid~-l~~; leql?rlX\~&NlX ^H\(@rNrNSPXlHDlS |\ :/0.) @NlX //=> ?(
@lrrPX ==!| zipp0.; @x not zeroid+ ==@h->hr ~&t, (^lryPX/~&lrrl2C minus^*p/~&rrr times*lrlPD)^/div@bzPrrPlXO ~&, @r ^|\~& ~&i&& :/0.)</lang>
test program: <lang Ursala>#cast %eLW
example = polydiv(<-42.,0.,-12.,1.>,<-3.,1.,0.,0.>)</lang> output:
( <-2.700000e+01,-9.000000e+00,1.000000e+00>, <-1.230000e+02>)