Division of Polynomials in python

Question

I am stuck with division of polynomials in python. Here is code that I modified. The while loop couldnt work. This code only output the original L as r. If I remove the while loop, only the remainder from first time division was outputted. I tried a bunch of ways to make it work, but all failed. Any suggestions will be greatly appreciated. Thanks!

def GetDegree(poly):
    while poly and poly[-1] == 0:
        poly.pop()   # normalize
    return len(poly)-1

def division(p1,p2):
    d1 = GetDegree(p1)
    d2 = GetDegree(p2)
    if d2 < 0 or d1<0:
        raise ZeroDivisionError
    if d1 > d2: 
        S,L = p2,p1#obtain S with lower degree, L with higher degree
    else: 
        S,L = p1,p2
    d1 = GetDegree(L)
    d2 = GetDegree(S)
    while d1>0:
            q = [0]*d1
        d = [0]*(d1 - d2) + S#shift short towards right by d1-d2
        mult = q[d1 - d2] = L[-1] / float(d[-1])#get the result by dividing the first term of the dividend by the highest term of the divisor
        d = [coef*mult for coef in d]#multiply the above result by short
        L = [fabs( coefL - coefd ) for coefL, coefd in zip(L, d)]#return a new long by subtracting long with d
        d1 = GetDegree(L)#return new d1
    r = L#return new long and keeping looping for there is no variable left and return as remainder
    return r

I want to input any random polynomials for the computation. However, when I modified it, the results still not right. Here is the test that I ran: num:[2,1,1,1] den:[1,1,2]. Print result was: quote:[0.25,0.5], rem:[1.75,0.25]. Here is the code that I modified for the case of input, based on the answer from PM 2Ring:

    def normalize(poly):
        while poly and poly[-1] == 0:
            poly.pop()
        if poly == []:
            poly.append(0)


    def poly_divmod(num, den):
        #Create normalized copies of the args
        num = num[:]
        normalize(num)
        den = den[:]
        normalize(den)

        if len(num) >= len(den):
            #Shift den towards right so it's the same degree as num
            shiftlen = len(num) - len(den)
            den = [0] * shiftlen + den
        else:
            return [0], num

        quot = []
        divisor = float(den[-1])
        for i in range(shiftlen + 1):
            #Get the next coefficient of the quotient.
            mult = num[-1] / divisor
            quot = [mult] + quot

            #Subtract mult * den from num, but don't bother if mult == 0
            #Note that when i==0, mult!=0; so quot is automatically normalized.
            if mult != 0:
                d = [mult * u for u in den]
                num = [u - v for u, v in zip(num, d)]

            num.pop()
            den.pop(0)

        normalize(num)
        return quot, num


    def test(num, den):
        print ("%s / %s ->" % (num, den))
        q, r = poly_divmod(num, den)
        print ("quot: %s, rem: %s\n" % (q, r))
        return q, r


    def main():
        degree = int(input('Enter the degree of your polynomial 1:'))
        num = []
        for i in range (0,degree+1):
            coefficient = int(input('Enter the coefficient for x^ %i ? ' %i))
            num.append(coefficient)
        degree = int(input('Enter the degree of your polynomial 2:'))
        den = []
        for i in range (0,degree+1):
            coefficient = int(input('Enter the coefficient for x^ %i ? ' %i))
            den.append(coefficient)
        test(num, den)



    if __name__ == '__main__':
        main()

Solution 1

I've modified your code slightly, so now it returns the quotient and remainder.

FWIW, it would be fairly easy to create a polynomial class, and then you could do polynomial arithmetic using standard operators and functions...

#! /usr/bin/env python

''' Polynomial long division

From http://stackoverflow.com/questions/26173058/division-of-polynomials-in-python

A polynomial is represented by a list of its coefficients, eg
5*x**3 + 4*x**2 + 1 -> [1, 0, 4, 5]

Modified by PM 2Ring 2014.10.03
'''

def normalize(poly):
    while poly and poly[-1] == 0:
        poly.pop()
    if poly == []:
        poly.append(0)


def poly_divmod(num, den):
    #Create normalized copies of the args
    num = num[:]
    normalize(num)
    den = den[:]
    normalize(den)

    if len(num) >= len(den):
        #Shift den towards right so it's the same degree as num
        shiftlen = len(num) - len(den)
        den = [0] * shiftlen + den
    else:
        return [0], num

    quot = []
    divisor = float(den[-1])
    for i in xrange(shiftlen + 1):
        #Get the next coefficient of the quotient.
        mult = num[-1] / divisor
        quot = [mult] + quot

        #Subtract mult * den from num, but don't bother if mult == 0
        #Note that when i==0, mult!=0; so quot is automatically normalized.
        if mult != 0:
            d = [mult * u for u in den]
            num = [u - v for u, v in zip(num, d)]

        num.pop()
        den.pop(0)

    normalize(num)
    return quot, num


def test(num, den):
    print "%s / %s ->" % (num, den)
    q, r = poly_divmod(num, den)
    print "quot: %s, rem: %s\n" % (q, r)
    return q, r


def main():
    num = [1, 5, 10, 10, 5, 1]
    den = [1, 2, 1]
    test(num, den)

    num = [5, 16, 10, 22, 7, 11, 1, 3]
    den = [1, 2, 1, 3]

    quot = [5, 1, 3, 0, 1]
    rem = [0, 5]

    q, r = test(num, den)
    assert quot == q
    assert rem == r


if __name__ == '__main__':
    main()

Solution 2

In case you're open to using an external library (I saw sympy mentioned above), numpy can easily solve this for you. numpy.polydiv is what you'd need. Example: https://numpy.org/doc/stable/reference/generated/numpy.polydiv.html

Comments

[Orangeblue]

Hi, Thanks very much for the help! I want to input any random polynomials for the computation. However, when I modified it, the results still not right. Here is the test that I ran: num:[2,1,1,1] den:[1,1,2]. Print result was: quote:[0.25,0.5], rem:[1.75,0.25].

[Orangeblue]

The code was only a part of all codes that suppose to run all add/sub/mult/division operation of polynomials. I have the polynomials class.

[PM 2Ring]

You had me worried for a moment, Orangeblue! But I just checked it, and it is correct. I verified it by doing the division in a polynomial class I wrote about 6 years ago, and verifying that num == quot * den + rem. And to double-check that my Poly class isn't buggy I also did the division with pencil & paper. (x³ + x² + x + 2) / (2x² + x + 1) = 0.5x + 0.25, remainder 0.25x + 1.75.

[DatRid]

Hi ZeldasLizard. Can you please extend your answer and describe what you changed and why, and how it helps the questioner ? Thank you! A side note: For formatting code, select the code and press CRTL + K