Polynomial long division: Difference between revisions
m (tighten up some formatting) |
m (a bit more formatting of the header matter) |
||
Line 1: | Line 1: | ||
{{task|Classic CS problems and programs}}{{Wikipedia}} |
{{task|Classic CS problems and programs}}{{Wikipedia}} |
||
<cite>In algebra, [[wp:Polynomial long division|polynomial long division]] is an algorithm for dividing a polynomial by another polynomial of the same or lower degree.</cite> |
:<cite>In algebra, [[wp:Polynomial long division|polynomial long division]] is an algorithm for dividing a polynomial by another polynomial of the same or lower degree.</cite> |
||
Let us suppose a polynomial is represented by a vector (e. |
Let us suppose a polynomial is represented by a vector, <math>x</math> (i.e., an ordered collection of [[wp:Coefficient|coefficients]]) so that the <math>i</math><sup>th</sup> element keeps the coefficient of <math>x^i</math>, 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. |
||
from left) followed by a multiplication of each element by the coefficient of the monomial. |
|||
Then a pseudocode for the polynomial long division could be: |
Then a pseudocode for the polynomial long division could be: |
||
Line 10: | Line 9: | ||
'''return''' the index of the last non-zero element of '''P'''; |
'''return''' the index of the last non-zero element of '''P'''; |
||
if all elements are 0, return -∞ |
if all elements are 0, return -∞ |
||
polynomial_long_division('''N''', '''D''') ''returns'' ('''q''', '''r'''): |
polynomial_long_division('''N''', '''D''') ''returns'' ('''q''', '''r'''): |
||
<span class="co1">// '''N''', '''D''', '''q''', '''r''' are vectors</span> |
<span class="co1">// '''N''', '''D''', '''q''', '''r''' are vectors</span> |
||
Line 33: | Line 32: | ||
* Error handling (for allocations or for wrong inputs) is not mandatory. |
* Error handling (for allocations or for wrong inputs) is not mandatory. |
||
'''Example for clarification''' |
'''Example for clarification''' |
||
<br> |
|||
This example is from Wikipedia, but changed to show how the given pseudocode works. |
This example is from Wikipedia, but changed to show how the given pseudocode works. |
||
Revision as of 21:47, 17 June 2009
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.
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.
Then a pseudocode for the polynomial long division 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.
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>
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>
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]
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>