greatest common divisor

gcd_poly.py

import cmath


def gcd_poly(p, q):

    if len(p) < len(q):
        p, q = q, p
    if not q:
        return p

    # normalize to monic polynomials
    p = [x / p[0] for x in p]
    q = [x / q[0] for x in q]

    # subtract polynomials
    for i in range(len(q)):
        p[i] -= q[i]

    # trim leading zeros
    while p and abs(p[0]) < 1e-12:
        del p[0]

    return gcd_poly(p, q)


def quadratic_roots(a, b, c):

    d = cmath.sqrt(b ** 2 - 4 * a * c)
    return (-b + d) / (2 * a), (-b - d) / (2 * a)


def eval_poly(p, x):

    e = 0
    for i in p:
        e = e * x + i
    return e


# polynomials
p = [0.533507195114693, -0.419279629666331, 0.196266328799736, 0]
q = [1, -1.57178622333742, 1.35338686531122, -0.578227837482342, 0.135335283236613]

# greatest common divisor with roots
d = gcd_poly(p, q)
r1, r2 = quadratic_roots(*d)

print(d)
print(r1)
print(r2)

# test roots
print(eval_poly(p, r1))
print(eval_poly(q, r2))

gdc.py

def gcd(a, b):

    if a < b:
        a, b = b, a
    if b == 0:
        return a
    return gcd(a % b, b)


print(gcd(6 * 3 * 5, 8 * 5 * 7))