<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def polynomial(A, x, b, coefficients, iterations=1):
 """Apply a polynomial smoother to the system Ax=b. Parameters A : sparse matrix Sparse NxN matrix x : ndarray Approximate solution (length N) b : ndarray Right-hand side (length N) coefficients : array_like Coefficients of the polynomial. See Notes section for details. iterations : int Number of iterations to perform Returns ------- Nothing, x will be modified in place. Notes ----- The smoother has the form x[:] = x + p(A) (b - A*x) where p(A) is a polynomial in A whose scalar coefficients are specified (in descending order) by argument 'coefficients'. - Richardson iteration p(A) = c_0: polynomial_smoother(A, x, b, [c_0]) - Linear smoother p(A) = c_1*A + c_0: polynomial_smoother(A, x, b, [c_1, c_0]) - Quadratic smoother p(A) = c_2*A^2 + c_1*A + c_0: polynomial_smoother(A, x, b, [c_2, c_1, c_0]) Here, Horner's Rule is applied to avoid computing A^k directly. For efficience, the method detects the case x = 0 one matrix-vector product is avoided (since (b - A*x) is b). Examples -------- """

    A, x, b = make_system(A, x, b, formats=None)

    for i in range(iterations):
        from pyamg.util.linalg import norm

        if norm(x) == 0:
            residual = b
        else:
            residual = (b - A*x)

        h = coefficients[0]*residual

        for c in coefficients[1:]:
            h = c*residual + A*h

        x += h