def _separate_sq(p):
    '\n    Helper function for ``_minimal_polynomial_sq``.\n\n    It selects a rational ``g`` such that the polynomial ``p``\n    consists of a sum of terms whose surds squared have gcd equal to ``g``\n    and a sum of terms with surds squared prime with ``g``;\n    then it takes the field norm to eliminate ``sqrt(g)``\n\n    See simplify.simplify.split_surds and polytools.sqf_norm.\n\n    Examples\n    ========\n\n    >>> p = -x + sqrt(2) + sqrt(3) + sqrt(7)\n    >>> p = _separate_sq(p)\n    >>> p\n    -x**2 + 2*sqrt(3)*x + 2*sqrt(7)*x - 2*sqrt(21) - 8\n    >>> p = _separate_sq(p)\n    >>> p\n    -x**4 + 4*sqrt(7)*x**3 - 32*x**2 + 8*sqrt(7)*x + 20\n    >>> p = _separate_sq(p)\n    >>> p\n    -x**8 + 48*x**6 - 536*x**4 + 1728*x**2 - 400\n\n    '

    def is_sqrt(expr):
        return (expr.is_Pow and (expr.exp == Rational(1, 2)))
    p = p.doit()
    a = []
    for y in p.args:
        if (not y.is_Mul):
            if is_sqrt(y):
                a.append((Integer(1), (y ** 2)))
            elif y.is_Atom:
                a.append((y, Integer(1)))
            else:
                raise NotImplementedError
        else:
            sifted = sift(y.args, is_sqrt)
            a.append((Mul(*sifted[False]), (Mul(*sifted[True]) ** 2)))
    a.sort(key=(lambda z: z[1]))
    if (a[(- 1)][1] == 1):
        return p
    surds = [z for (y, z) in a]
    for (i, si) in enumerate(surds):
        if (si != 1):
            break
    (_, b1, _) = _split_gcd(*surds[i:])
    a1 = []
    a2 = []
    for (y, z) in a:
        if (z in b1):
            a1.append((y * sqrt(z)))
        else:
            a2.append((y * sqrt(z)))
    p1 = Add(*a1)
    p2 = Add(*a2)
    return (_mexpand((p1 ** 2)) - _mexpand((p2 ** 2)))
Helper function for ``_minimal_polynomial_sq``.

It selects a rational ``g`` such that the polynomial ``p``
consists of a sum of terms whose surds squared have gcd equal to ``g``
and a sum of terms with surds squared prime with ``g``;
then it takes the field norm to eliminate ``sqrt(g)``

See simplify.simplify.split_surds and polytools.sqf_norm.

Examples
========

>>> p = -x + sqrt(2) + sqrt(3) + sqrt(7)
>>> p = _separate_sq(p)
>>> p
-x**2 + 2*sqrt(3)*x + 2*sqrt(7)*x - 2*sqrt(21) - 8
>>> p = _separate_sq(p)
>>> p
-x**4 + 4*sqrt(7)*x**3 - 32*x**2 + 8*sqrt(7)*x + 20
>>> p = _separate_sq(p)
>>> p
-x**8 + 48*x**6 - 536*x**4 + 1728*x**2 - 400
def _separate_sq(p):
    '\n    Helper function for ``_minimal_polynomial_sq``.\n\n    It selects a rational ``g`` such that the polynomial ``p``\n    consists of a sum of terms whose surds squared have gcd equal to ``g``\n    and a sum of terms with surds squared prime with ``g``;\n    then it takes the field norm to eliminate ``sqrt(g)``\n\n    See simplify.simplify.split_surds and polytools.sqf_norm.\n\n    Examples\n    ========\n\n    >>> p = -x + sqrt(2) + sqrt(3) + sqrt(7)\n    >>> p = _separate_sq(p)\n    >>> p\n    -x**2 + 2*sqrt(3)*x + 2*sqrt(7)*x - 2*sqrt(21) - 8\n    >>> p = _separate_sq(p)\n    >>> p\n    -x**4 + 4*sqrt(7)*x**3 - 32*x**2 + 8*sqrt(7)*x + 20\n    >>> p = _separate_sq(p)\n    >>> p\n    -x**8 + 48*x**6 - 536*x**4 + 1728*x**2 - 400\n\n    '

    def is_sqrt(expr):
        return (expr.is_Pow and (expr.exp == Rational(1, 2)))
    p = p.doit()
    a = []
    for y in p.args:
        if (not y.is_Mul):
            if is_sqrt(y):
                a.append((Integer(1), (y ** 2)))
            elif y.is_Atom:
                a.append((y, Integer(1)))
            else:
                raise NotImplementedError
        else:
            sifted = sift(y.args, is_sqrt)
            a.append((Mul(*sifted[False]), (Mul(*sifted[True]) ** 2)))
    a.sort(key=(lambda z: z[1]))
    if (a[(- 1)][1] == 1):
        return p
    surds = [z for (y, z) in a]
    for (i, si) in enumerate(surds):
        if (si != 1):
            break
    (_, b1, _) = _split_gcd(*surds[i:])
    a1 = []
    a2 = []
    for (y, z) in a:
        if (z in b1):
            a1.append((y * sqrt(z)))
        else:
            a2.append((y * sqrt(z)))
    p1 = Add(*a1)
    p2 = Add(*a2)
    return (_mexpand((p1 ** 2)) - _mexpand((p2 ** 2)))