# long division
def divide(A, B):
    dA = A.degree()
    dB = B.degree()

    if (dB > dA):
        return 0, A
    else:
        q = x^(dA-dB)*A.leading_coefficient()/B.leading_coefficient()
        newA = A - q*B
        qr = divide(newA, B)
        return qr[0]+q, qr[1]


# truncation up to x^(deg-1) inclusive
def trunc(B, deg):
    return B.parent()(B.coefficients()[0:deg])


# Newton iteration for inverse
# computes up to x^(deg-1) inclusive
def inverseNewton(B, deg):
    C = 1/B[0]
    i = 1
    while (i < deg):
        i = 2*i
        C = trunc(C * (2 - trunc(B, i)*C), i)
    return trunc(C, deg)


# division using Newton iteration
def divideNewton(A, B):
    dA = A.degree()
    dB = B.degree()
    revB = B.reverse()
    inv_revB = inverseNewton(revB, dA-dB+1)
    revQ = trunc(inv_revB*A.reverse(), dA-dB+1)
    Q = revQ.reverse()
    R = A-B*Q
    return Q, R

#################################
## test
#################################

# base ring
U.<x> = PolynomialRing(QQ)

# input polynomials
A = 28*x^5 + x^4 - 17*x^3 + 22*x^2 - 44*x + 10
B = x^2 + 10*x - 1 


# long division
Q, R = divide(A, B)
print(A == (Q*B+R))


# division using Newton
Q2, R2 = divideNewton(A, B)
print(Q==Q2)
print(R==R2)