'Construct a minimal domain that contains elements of ``K0`` and ``K1``.
Known domains (from smallest to largest):
- ``GF(p)``
- ``ZZ``
- ``QQ``
- ``RR(prec, tol)``
- ``CC(prec, tol)``
- ``ALG(a, b, c)``
- ``K[x, y, z]``
- ``K(x, y, z)``
- ``EX``'

def unify(K0, K1, symbols=None):

    if (symbols is not None):
        return K0.unify_with_symbols(K1, symbols)
    if (K0 == K1):
        return K0
    if K0.is_EX:
        return K0
    if K1.is_EX:
        return K1
    if (K0.is_Composite or K1.is_Composite):
        K0_ground = (K0.dom if K0.is_Composite else K0)
        K1_ground = (K1.dom if K1.is_Composite else K1)
        K0_symbols = (K0.symbols if K0.is_Composite else ())
        K1_symbols = (K1.symbols if K1.is_Composite else ())
        domain = K0_ground.unify(K1_ground)
        symbols = _unify_gens(K0_symbols, K1_symbols)
        order = (K0.order if K0.is_Composite else K1.order)
        if (((K0.is_FractionField and K1.is_PolynomialRing) or (K1.is_FractionField and K0.is_PolynomialRing)) and ((not K0_ground.is_Field) or (not K1_ground.is_Field)) and domain.is_Field):
            domain = domain.get_ring()
        if (K0.is_Composite and ((not K1.is_Composite) or K0.is_FractionField or K1.is_PolynomialRing)):
            cls = K0.__class__
        else:
            cls = K1.__class__
        from sympy.polys.domains.old_polynomialring import GlobalPolynomialRing
        if (cls == GlobalPolynomialRing):
            return cls(domain, symbols)
        return cls(domain, symbols, order)
    def mkinexact(cls, K0, K1):
        prec = max(K0.precision, K1.precision)
        tol = max(K0.tolerance, K1.tolerance)
        return cls(prec=prec, tol=tol)
    if (K0.is_ComplexField and K1.is_ComplexField):
        return mkinexact(K0.__class__, K0, K1)
    if (K0.is_ComplexField and K1.is_RealField):
        return mkinexact(K0.__class__, K0, K1)
    if (K0.is_RealField and K1.is_ComplexField):
        return mkinexact(K1.__class__, K1, K0)
    if (K0.is_RealField and K1.is_RealField):
        return mkinexact(K0.__class__, K0, K1)
    if (K0.is_ComplexField or K0.is_RealField):
        return K0
    if (K1.is_ComplexField or K1.is_RealField):
        return K1
    if (K0.is_AlgebraicField and K1.is_AlgebraicField):
        return K0.__class__(K0.dom.unify(K1.dom), *_unify_gens(K0.orig_ext, K1.orig_ext))
    elif K0.is_AlgebraicField:
        return K0
    elif K1.is_AlgebraicField:
        return K1
    if K0.is_RationalField:
        return K0
    if K1.is_RationalField:
        return K1
    if K0.is_IntegerRing:
        return K0
    if K1.is_IntegerRing:
        return K1
    if (K0.is_FiniteField and K1.is_FiniteField):
        return K0.__class__(max(K0.mod, K1.mod, key=default_sort_key))
    from sympy.polys.domains import EX
    return EX
