import numpy as np

def polygcd(a,b):
    '''return monic GCD of polynomials a and b'''
    pa = a
    pb = b
    M = lambda x: x/x[0]
    # find gcd of a and b
    while len(pb) > 1 or pb[0] != 0:
        # Danger Will Robinson! requires numerical equality
        q,r = np.polydiv(pa,M(pb))
        pa = pb
        pb = r
    return M(pa)

def polylcm(a,b):
    '''return (Ka,Kb,c) such that c = LCM(a,b) = Ka*a = Kb*b'''        
    gcd = polygcd(a,b)
    Ka,_ = np.polydiv(b,gcd)
    Kb,_ = np.polydiv(a,gcd)
    return (Ka,Kb,np.polymul(Ka,a))