Algorithms to multiply, divide, add, and subtract polynomials. Also some various functions to go with it. Note that we have that a polynomial is a list [t1, t2, ...] where we have t1 + xt2 + x^2 t3 + ...

polynomialdivision.py

import warnings

#this is just for notational simplicity
def degree(f):
    try:
        return len(f)-1
    except:
        return 0
# -----------------------------------------------------------------------------------------------------------------------
#also for notational simplicity
def lead(f):
    try:
        return f[len(f) -1]
    except:
        return 0
# -----------------------------------------------------------------------------------------------------------------------
#generic function to copy lists easily (i.e. copy polynomials)
def copy_list(f):
    new_poly = []
    for i in range(len(f)):
        new_poly.append(f[i])
    return new_poly
# -----------------------------------------------------------------------------------------------------------------------
#If leading coefficient is 0, eliminates it
def check(f):
    new_poly = copy_list(f)
    if len(new_poly) == 0:
        return 0
    elif new_poly[len(new_poly)-1] == 0:
        new_poly.pop()
        return check(new_poly)
    else:
        return f
# -----------------------------------------------------------------------------------------------------------------------
#Prints a polynomial in a readable manner
def print_polynomial(f):
    stri = ""
    if len(f) > 1:
        for i in range(len(f)-1):
            if f[i] != 0:
                if i == 0:
                    stri += str(f[i])
                    stri += " + "
                elif f[i] == 1:
                    stri += "x^"
                    stri += str(i)
                    stri += " + "
                else:
                    stri += str(f[i])
                    stri += "x^"
                    stri += str(i)
                    stri += " + "
        if f[len(f)-1] != 1:
            stri += str(f[len(f)-1])
            stri += "x^"
            stri += str(len(f)-1)
        else:
            stri += "x^"
            stri += str(len(f)-1)
    else:
        stri += str(f[len(f)-1])
    print(stri)

# -----------------------------------------------------------------------------------------------------------------------

#Multiplies a scalar (i.e. a constant) to all the values of a polynomial
def scalar_mult_polynomial(scalar, f):
    #Multiply a scalar to all values of polynomials
    new_poly = []
    for i in range(len(f)):
        new_poly.append(scalar* f[i])
    return new_poly

# -----------------------------------------------------------------------------------------------------------------------

#subtracts g from f
def polynomial_subtraction(f,g):
    #subtracting g from f
    new_poly = []
    if type(f) == int:
        f = convert_int(f)
    if type(g) == int:
        g = convert_int(g)
    if len(f) > len(g):
        for i in range(len(g)):
            new_poly.append(f[i] - g[i])
        for j in range(len(g), len(f)):
            new_poly.append(f[j])
    elif len(f) == len(g):
        for i in range(len(f)):
            new_poly.append(f[i] - g[i])
    else:
        for i in range(len(f)):
            new_poly.append(f[i] - g[i])
        for j in range(len(f), len(g)):
            new_poly.append(-g[j])
    return new_poly

# -----------------------------------------------------------------------------------------------------------------------

#adds two polynomials together
def polynomial_addition(f,g):
    #adds g to f
    if type(f) != list:
        f = [f]
    if type(g) != list:
        g = [g]
    new_poly = []
    if len(f) > len(g):
        for i in range(len(g)):
            new_poly.append(f[i] + g[i])
        for j in range(len(g), len(f)):
            new_poly.append(f[j])
    elif len(f) == len(g):
        for i in range(len(f)):
            new_poly.append(f[i] + g[i])
    else:
        for i in range(len(f)):
            new_poly.append(f[i] + g[i])
        for j in range(len(f), len(g)):
            new_poly.append(-g[j])
    return new_poly

# -----------------------------------------------------------------------------------------------------------------------

# Uses a general form of foil to multiply polynomials (see - convolution)
def multiply_polynomials(f,g):
    new_poly = []
    for i in range(len(f) * len(g)):
        new_poly.append(0)
    for i in range(len(f)):
        for j in range(len(g)):
            index = i+j
            new_poly[index] += f[i] * g[j]
    new_poly = check(new_poly)
    return new_poly

# -----------------------------------------------------------------------------------------------------------------------

#Takes a polynomial to a power.
def polynomial_power(f, power):
    x = copy_list(f)
    for i in range(power-1):
        x = multiply_polynomials(f, x)
    return x

# -----------------------------------------------------------------------------------------------------------------------

#Divides polynomials utilizing the Euclidean algorithm. Returns q and r, where r is the remainder and q is the polynomial divisor.
def divide_polynomials(f,g):
    if degree(g) == 0 and g[0] == 0:
        warnings.warn("g cannot be the 0 polynomial")
    q = [0]
    r = f
    while r != 0 and degree(r) >= degree(g):
        ti = []
        t = float(lead(r)) / float(lead(g))
        x = degree(r) - degree(g)
        while degree(q) < x:
            q.append(0)
        while degree(ti) < x:
            ti.append(0)
        if x == 0:
            q.append(0)
            ti.append(0)
        q[x] = t
        ti[x] = t
        r = polynomial_subtraction(r, multiply_polynomials(ti, g))

        r = check(r)
    return (q,r)

# -----------------------------------------------------------------------------------------------------------------------

#converts a polynomial to a monic polynomial (i.e. the leading coefficient is 1)
def convert_to_monic(polynomial):
    if lead(polynomial) != 1:
        new_poly = []
        for i in range(len(polynomial)):
            new_poly.append(polynomial[i]/lead(polynomial))
        return new_poly
    else:
        return polynomial

# -----------------------------------------------------------------------------------------------------------------------

#Converts an integer to a polynomial
def convert_int(x):
    return [x]

# -----------------------------------------------------------------------------------------------------------------------

#Find gcd of polynomials
def gcd(f,g):
    if degree(g) == 0 and g[0] == 0:
        warnings.warn("g cannot be the 0 polynomial")
    if degree(f) >= degree(g):
        ma = f
        mi = g
    else:
        ma = g
        mi = f
    q, r = divide_polynomials(ma,mi)
    if r== 0:
        return mi
    else:
        r1 = convert_to_monic(r)
        return gcd(mi, r1)

# -----------------------------------------------------------------------------------------------------------------------
#Runs a test case for example.
if __name__ == "__main__":
    f=[-24,2,1]
    g= [-18,3,1]
    print_polynomial(gcd(f,g))

For more on the algorithm, see https://en.wikipedia.org/wiki/Polynomial_long_division