def divide(dividend, divisor, quotient=None):
    d1, d2 = len(dividend), len(divisor)
    if quotient is None:
        terms = d1 - d2 + 1
        quotient = [0]*terms
    pos = 0
    for i in range(d1):
        if dividend[i] != 0:
            pos = i
            break
    power = d1 - 1 - pos
    factor, mini = dividend[pos] / divisor[0], d2 - 1
    if power - mini >= 0:
        quotient[len(quotient) - (power - mini) - 1] += factor
        for i in range(d2):
            cur = divisor[i] * factor
            curpower = d2 - 1 - i + power - mini
            dividend[d1 - curpower - 1] -= cur
        for j in range(d1 - d2 + 1):
            if dividend[j] != 0:
                return divide(dividend, divisor, quotient)
        return "Quotient is: %s (remainder:%s) " % (pol(quotient), pol(dividend))
    else:
        return "Quotient is: %s (remainder:%s) " % (pol(quotient), pol(dividend))


def pol(x):
    val = "  "
    for i in range(len(x)):
        if x[i] != 0:
            if i == len(x) - 1:
                val += "+ " + str(int(x[i]))
            else:
                val += "+ " + str(int(x[i])) + " x^" + str(int(len(x)) - i - 1) + " "
    if val == "  ":
        return " 0"
    return val.replace("  +", "").replace("+ -", "- ").replace(" 1 ", "").replace("x^1", "x").replace(" x", "x")


print(divide([1, 3, 3, 1], [1, 1]))  # Divides (x^3 + 3x^2 + 3x + 1) by (x + 1)
print(divide([2, -5, -1], [1, -3]))  # Divides (2x^2 - 5x - 1) by (x - 3)
print(divide([4, 3, 0, 2, 1], [1, 1, 2]))  # Divides (4x^4 + 3x^3 + 2x + 1) by (x^2 + x + 2)
print(divide([2, 13, 15, 0, 0], [2, 3]))  # Divides (2x^4 + 13x^3 + 15x^2) by (2x + 3)