Polynomial long division
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 (e.g. see here), then the i-th element keeps the coefficient of xi, 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>