def primitive_element(extension, **args):
    'Construct a common number field for all extensions.\n\n    References\n    ==========\n\n    * :cite:`Yokoyama1989primitive`\n    * :cite:`Arno1996alg`\n\n    '
    if (not extension):
        raise ValueError("can't compute primitive element for empty extension")
    extension = list(uniq(extension))
    x = Dummy('x')
    domain = args.get('domain', QQ)
    F = [minimal_polynomial(e, domain=domain) for e in extension]
    Y = [p.gen for p in F]
    for u in range(1, (((len(F) - 1) * math.prod((f.degree() for f in F))) + 1)):
        coeffs = [(u ** n) for n in range(len(Y))]
        f = (x - sum(((c * y) for (c, y) in zip(coeffs, Y))))
        (*H, g) = groebner((F + [f]), *(Y + [x]), domain=domain)
        for (i, (h, y)) in enumerate(zip(H, Y)):
            H[i] = (y - h).eject(*Y).retract(field=True)
            if (not (H[i].domain.is_RationalField or H[i].domain.is_AlgebraicField)):
                break
            else:
                H[i] = H[i].set_domain(domain)
        else:
            g = g.eject(*Y).set_domain(domain)
            break
    else:
        if (len(F) == 1):
            (g, coeffs, H) = (F[0].replace(x), [Integer(1)], [x.as_poly(domain=domain)])
        else:
            raise RuntimeError('run out of coefficient configurations')
    (_, factors) = factor_list(g, domain=domain)
    t = sum(((c * e) for (c, e) in zip(coeffs, extension)))
    g = _choose_factor(factors, x, t, dom=domain)
    H = [h.rem(g).rep.all_coeffs() for h in H]
    (_, g) = PurePoly(g).clear_denoms(convert=True)
    if (g.LC() != 1):
        for d in divisors(g.LC())[1:]:
            new_g = ((g.compose((g.gen / d).as_poly()) * (d ** g.degree())) // d)
            (_, new_g) = new_g.monic().clear_denoms(convert=True)
            if (new_g.LC() == 1):
                g = new_g
                H = [[(c / (d ** n)) for (n, c) in enumerate(h)] for h in H]
                coeffs = [(c * d) for c in coeffs]
                break
    return (g, list(coeffs), H)
Construct a common number field for all extensions.

References
==========

* :cite:`Yokoyama1989primitive`
* :cite:`Arno1996alg`
def primitive_element(extension, **args):
    'Construct a common number field for all extensions.\n\n    References\n    ==========\n\n    * :cite:`Yokoyama1989primitive`\n    * :cite:`Arno1996alg`\n\n    '
    if (not extension):
        raise ValueError("can't compute primitive element for empty extension")
    extension = list(uniq(extension))
    x = Dummy('x')
    domain = args.get('domain', QQ)
    F = [minimal_polynomial(e, domain=domain) for e in extension]
    Y = [p.gen for p in F]
    for u in range(1, (((len(F) - 1) * math.prod((f.degree() for f in F))) + 1)):
        coeffs = [(u ** n) for n in range(len(Y))]
        f = (x - sum(((c * y) for (c, y) in zip(coeffs, Y))))
        (*H, g) = groebner((F + [f]), *(Y + [x]), domain=domain)
        for (i, (h, y)) in enumerate(zip(H, Y)):
            H[i] = (y - h).eject(*Y).retract(field=True)
            if (not (H[i].domain.is_RationalField or H[i].domain.is_AlgebraicField)):
                break
            else:
                H[i] = H[i].set_domain(domain)
        else:
            g = g.eject(*Y).set_domain(domain)
            break
    else:
        if (len(F) == 1):
            (g, coeffs, H) = (F[0].replace(x), [Integer(1)], [x.as_poly(domain=domain)])
        else:
            raise RuntimeError('run out of coefficient configurations')
    (_, factors) = factor_list(g, domain=domain)
    t = sum(((c * e) for (c, e) in zip(coeffs, extension)))
    g = _choose_factor(factors, x, t, dom=domain)
    H = [h.rem(g).rep.all_coeffs() for h in H]
    (_, g) = PurePoly(g).clear_denoms(convert=True)
    if (g.LC() != 1):
        for d in divisors(g.LC())[1:]:
            new_g = ((g.compose((g.gen / d).as_poly()) * (d ** g.degree())) // d)
            (_, new_g) = new_g.monic().clear_denoms(convert=True)
            if (new_g.LC() == 1):
                g = new_g
                H = [[(c / (d ** n)) for (n, c) in enumerate(h)] for h in H]
                coeffs = [(c * d) for c in coeffs]
                break
    return (g, list(coeffs), H)