def _cholesky(self, hermitian=True):
     """Helper function of cholesky.
     Without the error checks.
     To be used privately.
     Implements the Cholesky-Banachiewicz algorithm.
     Returns L such that L*L.H == self if hermitian flag is True,
     or L*L.T == self if hermitian is False.
     """
     L = zeros(self.rows, self.rows)
     if hermitian:
         for i in range(self.rows):
             for j in range(i):
                 L[i, j] = (1 / L[j, j])*expand_mul(self[i, j] -
                     sum(L[i, k]*L[j, k].conjugate() for k in range(j)))
             Lii2 = expand_mul(self[i, i] -
                 sum(L[i, k]*L[i, k].conjugate() for k in range(i)))
             if Lii2.is_positive is False:
                 raise ValueError("Matrix must be positive-definite")
             L[i, i] = sqrt(Lii2)
     else:
         for i in range(self.rows):
             for j in range(i):
                 L[i, j] = (1 / L[j, j])*(self[i, j] -
                     sum(L[i, k]*L[j, k] for k in range(j)))
             L[i, i] = sqrt(self[i, i] -
                 sum(L[i, k]**2 for k in range(i)))
     return self._new(L)

 def _LDLdecomposition(self, hermitian=True):
     """Helper function of LDLdecomposition.
     Without the error checks.
     To be used privately.
     Returns L and D such that L*D*L.H == self if hermitian flag is True,
     or L*D*L.T == self if hermitian is False.
     """
     # https://en.wikipedia.org/wiki/Cholesky_decomposition#LDL_decomposition_2
     D = zeros(self.rows, self.rows)
     L = eye(self.rows)
     if hermitian:
         for i in range(self.rows):
             for j in range(i):
                 L[i, j] = (1 / D[j, j])*expand_mul(self[i, j] - sum(
                     L[i, k]*L[j, k].conjugate()*D[k, k] for k in range(j)))
             D[i, i] = expand_mul(self[i, i] -
                 sum(L[i, k]*L[i, k].conjugate()*D[k, k] for k in range(i)))
             if D[i, i].is_positive is False:
                 raise ValueError("Matrix must be positive-definite")
     else:
         for i in range(self.rows):
             for j in range(i):
                 L[i, j] = (1 / D[j, j])*(self[i, j] - sum(
                     L[i, k]*L[j, k]*D[k, k] for k in range(j)))
             D[i, i] = self[i, i] - sum(L[i, k]**2*D[k, k] for k in range(i))
     return self._new(L), self._new(D)

 def valid(x):
     # this is used to see if gen=x satisfies the
     # relational by substituting it into the
     # expanded form and testing against 0, e.g.
     # if expr = x*(x + 1) < 2 then e = x*(x + 1) - 2
     # and expanded_e = x**2 + x - 2; the test is
     # whether a given value of x satisfies
     # x**2 + x - 2 < 0
     #
     # expanded_e, expr and gen used from enclosing scope
     v = expanded_e.subs(gen, expand_mul(x))
     try:
         r = expr.func(v, 0)
     except TypeError:
         r = S.false
     if r in (S.true, S.false):
         return r
     if v.is_real is False:
         return S.false
     else:
         v = v.n(2)
         if v.is_comparable:
             return expr.func(v, 0)
         # not comparable or couldn't be evaluated
         raise NotImplementedError(
             'relationship did not evaluate: %s' % r)

def test_issue_11230():
    # a specific test that always failed
    a, b, f, k, l, i = symbols('a b f k l i')
    p = [a*b*f*k*l, a*i*k**2*l, f*i*k**2*l]
    R, C = cse(p)
    assert not any(i.is_Mul for a in C for i in a.args)

    # random tests for the issue
    from random import choice
    from sympy.core.function import expand_mul
    s = symbols('a:m')
    # 35 Mul tests, none of which should ever fail
    ex = [Mul(*[choice(s) for i in range(5)]) for i in range(7)]
    for p in subsets(ex, 3):
        p = list(p)
        R, C = cse(p)
        assert not any(i.is_Mul for a in C for i in a.args)
        for ri in reversed(R):
            for i in range(len(C)):
                C[i] = C[i].subs(*ri)
        assert p == C
    # 35 Add tests, none of which should ever fail
    ex = [Add(*[choice(s[:7]) for i in range(5)]) for i in range(7)]
    for p in subsets(ex, 3):
        p = list(p)
        was = R, C = cse(p)
        assert not any(i.is_Add for a in C for i in a.args)
        for ri in reversed(R):
            for i in range(len(C)):
                C[i] = C[i].subs(*ri)
        # use expand_mul to handle cases like this:
        # p = [a + 2*b + 2*e, 2*b + c + 2*e, b + 2*c + 2*g]
        # x0 = 2*(b + e) is identified giving a rebuilt p that
        # is now `[a + 2*(b + e), c + 2*(b + e), b + 2*c + 2*g]`
        assert p == [expand_mul(i) for i in C]

def sqrtdenest(expr, max_iter=3):
    """Denests sqrts in an expression that contain other square roots
    if possible, otherwise returns the expr unchanged. This is based on the
    algorithms of [1].

    Examples
    ========

    >>> from sympy.simplify.sqrtdenest import sqrtdenest
    >>> from sympy import sqrt
    >>> sqrtdenest(sqrt(5 + 2 * sqrt(6)))
    sqrt(2) + sqrt(3)

    See Also
    ========
    sympy.solvers.solvers.unrad

    References
    ==========
    [1] http://researcher.watson.ibm.com/researcher/files/us-fagin/symb85.pdf

    [2] D. J. Jeffrey and A. D. Rich, 'Symplifying Square Roots of Square Roots
    by Denesting' (available at http://www.cybertester.com/data/denest.pdf)

    """
    expr = expand_mul(sympify(expr))
    for i in range(max_iter):
        z = _sqrtdenest0(expr)
        if expr == z:
            return expr
        expr = z
    return expr

    def _eval_as_leading_term(self, x):
        from sympy import expand_mul, factor_terms

        old = self

        expr = expand_mul(self)
        if not expr.is_Add:
            return expr.as_leading_term(x)

        infinite = [t for t in expr.args if t.is_infinite]

        expr = expr.func(*[t.as_leading_term(x) for t in expr.args]).removeO()
        if not expr:
            # simple leading term analysis gave us 0 but we have to send
            # back a term, so compute the leading term (via series)
            return old.compute_leading_term(x)
        elif expr is S.NaN:
            return old.func._from_args(infinite)
        elif not expr.is_Add:
            return expr
        else:
            plain = expr.func(*[s for s, _ in expr.extract_leading_order(x)])
            rv = factor_terms(plain, fraction=False)
            rv_simplify = rv.simplify()
            # if it simplifies to an x-free expression, return that;
            # tests don't fail if we don't but it seems nicer to do this
            if x not in rv_simplify.free_symbols:
                if rv_simplify.is_zero and plain.is_zero is not True:
                    return (expr - plain)._eval_as_leading_term(x)
                return rv_simplify
            return rv

def sqrtdenest(expr, max_iter=3):
    """Denests sqrts in an expression that contain other square roots
    if possible, otherwise returns the expr unchanged. This is based on the
    algorithms of [1].

    Examples
    ========

    >>> from sympy.simplify.sqrtdenest import sqrtdenest
    >>> from sympy import sqrt
    >>> sqrtdenest(sqrt(5 + 2 * sqrt(6)))
    sqrt(2) + sqrt(3)

    See Also
    ========
    sympy.solvers.solvers.unrad

    References
    ==========
    [1] http://www.almaden.ibm.com/cs/people/fagin/symb85.pdf
    """
    expr = expand_mul(sympify(expr))
    for i in range(max_iter):
        z = _sqrtdenest0(expr)
        if expr == z:
            return expr
        expr = z
    return expr

 def _bell_poly(n, prev):
     s = 1
     a = 1
     for k in xrange(2, n + 1):
         a = a * (n - k + 1) // (k - 1)
         s += a * prev[k - 1]
     return expand_mul(_sym * s)

    def _bell_incomplete_poly(n, k, symbols):
        r"""
        The second kind of Bell polynomials (incomplete Bell polynomials).

        Calculated by recurrence formula:

        .. math:: B_{n,k}(x_1, x_2, \dotsc, x_{n-k+1}) =
                \sum_{m=1}^{n-k+1}
                \x_m \binom{n-1}{m-1} B_{n-m,k-1}(x_1, x_2, \dotsc, x_{n-m-k})

        where
            B_{0,0} = 1;
            B_{n,0} = 0; for n>=1
            B_{0,k} = 0; for k>=1

        """
        if (n == 0) and (k == 0):
            return S.One
        elif (n == 0) or (k == 0):
            return S.Zero
        s = S.Zero
        a = S.One
        for m in xrange(1, n - k + 2):
            s += a * bell._bell_incomplete_poly(
                n - m, k - 1, symbols) * symbols[m - 1]
            a = a * (n - m) / m
        return expand_mul(s)

 def _LDLdecomposition(self):
     """Helper function of LDLdecomposition.
     Without the error checks.
     To be used privately.
     """
     # https://en.wikipedia.org/wiki/Cholesky_decomposition#LDL_decomposition_2
     D = zeros(self.rows, self.rows)
     L = eye(self.rows)
     for i in range(self.rows):
         for j in range(i):
             L[i, j] = (1 / D[j, j])*expand_mul(self[i, j] - sum(
                 L[i, k]*L[j, k].conjugate()*D[k, k] for k in range(j)))
         D[i, i] = expand_mul(self[i, i] -
             sum(L[i, k]*L[i, k].conjugate()*D[k, k] for k in range(i)))
         if D[i, i].is_positive is False:
             raise ValueError("Matrix must be positive-definite")
     return self._new(L), self._new(D)

 def _cholesky(self):
     """Helper function of cholesky.
     Without the error checks.
     To be used privately.
     Implements the Cholesky-Banachiewicz algorithm.
     """
     L = zeros(self.rows, self.rows)
     for i in range(self.rows):
         for j in range(i):
             L[i, j] = (1 / L[j, j])*expand_mul(self[i, j] -
                 sum(L[i, k]*L[j, k].conjugate() for k in range(j)))
         Lii2 = expand_mul(self[i, i] -
             sum(L[i, k]*L[i, k].conjugate() for k in range(i)))
         if Lii2.is_positive is False:
             raise ValueError("Matrix must be positive-definite")
         L[i, i] = sqrt(Lii2)
     return self._new(L)

def convolution_fwht(a, b):
    """
    Performs dyadic (*bitwise-XOR*) convolution using Fast Walsh Hadamard
    Transform.

    The convolution is automatically padded to the right with zeros, as the
    *radix-2 FWHT* requires the number of sample points to be a power of 2.

    Parameters
    ==========

    a, b : iterables
        The sequences for which convolution is performed.

    Examples
    ========

    >>> from sympy import symbols, S, I
    >>> from sympy.discrete.convolutions import convolution_fwht

    >>> u, v, x, y = symbols('u v x y')
    >>> convolution_fwht([u, v], [x, y])
    [u*x + v*y, u*y + v*x]

    >>> convolution_fwht([2, 3], [4, 5])
    [23, 22]
    >>> convolution_fwht([2, 5 + 4*I, 7], [6*I, 7, 3 + 4*I])
    [56 + 68*I, -10 + 30*I, 6 + 50*I, 48 + 32*I]

    >>> convolution_fwht([S(33)/7, S(55)/6, S(7)/4], [S(2)/3, 5])
    [2057/42, 1870/63, 7/6, 35/4]

    References
    ==========

    .. [1] https://www.radioeng.cz/fulltexts/2002/02_03_40_42.pdf
    .. [2] https://en.wikipedia.org/wiki/Hadamard_transform

    """

    if not a or not b:
        return []

    a, b = a[:], b[:]
    n = max(len(a), len(b))

    if n&(n - 1): # not a power of 2
        n = 2**n.bit_length()

    # padding with zeros
    a += [S.Zero]*(n - len(a))
    b += [S.Zero]*(n - len(b))

    a, b = fwht(a), fwht(b)
    a = [expand_mul(x*y) for x, y in zip(a, b)]
    a = ifwht(a)

    return a

 def _combine_inverse(lhs, rhs):
     """
     Returns lhs - rhs, but treats oo like a symbol so oo - oo
     returns 0, instead of a nan.
     """
     from sympy.core.function import expand_mul
     from sympy.core.symbol import Dummy
     inf = (S.Infinity, S.NegativeInfinity)
     if lhs.has(*inf) or rhs.has(*inf):
         oo = Dummy('oo')
         reps = {
             S.Infinity: oo,
             S.NegativeInfinity: -oo}
         ireps = dict([(v, k) for k, v in reps.items()])
         eq = expand_mul(lhs.xreplace(reps) - rhs.xreplace(reps))
         if eq.has(oo):
             eq = eq.replace(
                 lambda x: x.is_Pow and x.base == oo,
                 lambda x: x.base)
         return eq.xreplace(ireps)
     else:
         return expand_mul(lhs - rhs)

 def eval(cls, p, q):
     from sympy.simplify.simplify import nsimplify
     if q.is_Number:
         float = not q.is_Rational
         pnew = expand_mul(p)
         if pnew.is_Number:
             float = float or not pnew.is_Rational
             if not float:
                 return pnew % q
             return Float(nsimplify(pnew) % nsimplify(q))
         elif pnew.is_Add and pnew.args[0].is_Number:
             r, p = pnew.as_two_terms()
             p += Mod(r, q)
     return Mod(p, q, evaluate=False)

def convolution_fft(a, b, dps=None):
    """
    Performs linear convolution using Fast Fourier Transform.

    Parameters
    ==========

    a, b : iterables
        The sequences for which convolution is performed.
    dps : Integer
        Specifies the number of decimal digits for precision.

    Examples
    ========

    >>> from sympy import S, I
    >>> from sympy.discrete.convolutions import convolution_fft

    >>> convolution_fft([2, 3], [4, 5])
    [8, 22, 15]
    >>> convolution_fft([2, 5], [6, 7, 3])
    [12, 44, 41, 15]
    >>> convolution_fft([1 + 2*I, 4 + 3*I], [S(5)/4, 6])
    [5/4 + 5*I/2, 11 + 63*I/4, 24 + 18*I]

    References
    ==========

    .. [1] https://en.wikipedia.org/wiki/Convolution_theorem
    .. [2] https://en.wikipedia.org/wiki/Discrete_Fourier_transform_(general%29

    """

    a, b = a[:], b[:]
    n = m = len(a) + len(b) - 1 # convolution size

    if n > 0 and n&(n - 1): # not a power of 2
        n = 2**n.bit_length()

    # padding with zeros
    a += [S.Zero]*(n - len(a))
    b += [S.Zero]*(n - len(b))

    a, b = fft(a, dps), fft(b, dps)
    a = [expand_mul(x*y) for x, y in zip(a, b)]
    a = ifft(a, dps)[:m]

    return a

def _fourier_transform(seq, dps, inverse=False):
    """Utility function for the Discrete Fourier Transform"""

    if not iterable(seq):
        raise TypeError("Expected a sequence of numeric coefficients "
                        "for Fourier Transform")

    a = [sympify(arg) for arg in seq]
    if any(x.has(Symbol) for x in a):
        raise ValueError("Expected non-symbolic coefficients")

    n = len(a)
    if n < 2:
        return a

    b = n.bit_length() - 1
    if n&(n - 1): # not a power of 2
        b += 1
        n = 2**b

    a += [S.Zero]*(n - len(a))
    for i in range(1, n):
        j = int(ibin(i, b, str=True)[::-1], 2)
        if i < j:
            a[i], a[j] = a[j], a[i]

    ang = -2*pi/n if inverse else 2*pi/n

    if dps is not None:
        ang = ang.evalf(dps + 2)

    w = [cos(ang*i) + I*sin(ang*i) for i in range(n // 2)]

    h = 2
    while h <= n:
        hf, ut = h // 2, n // h
        for i in range(0, n, h):
            for j in range(hf):
                u, v = a[i + j], expand_mul(a[i + j + hf]*w[ut * j])
                a[i + j], a[i + j + hf] = u + v, u - v
        h *= 2

    if inverse:
        a = [(x/n).evalf(dps) for x in a] if dps is not None \
                            else [x/n for x in a]

    return a

    def projection(self, pt):
        """Projection of any point on the plane
           Will result in a 3D point.

        Parameters
        ==========

        Point or Point3D

        Returns
        =======

        Point3D

        Notes
        =====

        For the interaction between 2D and 3D points, you should convert the 2D
        point to 3D by using this method. For example for finding the distance
        between Point(1, 2) and Point3D(1, 2, 3) you can convert the 2D point
        to Point3D(1, 2, 0) by projecting it to plane which is parallel to xy
        plane.

        Examples
        ========

        >>> from sympy import Plane, Point, Point3D
        >>> a = Plane(Point3D(1, 1, 1), normal_vector=(1, 1, 1))
        >>> b = Point(1, 2)
        >>> a.projection(b)
        Point3D(1, 2, 0)
        >>> c = Point3D(1, 1, 1)
        >>> a.projection(c)
        Point3D(4/3, 4/3, 4/3)

        """
        x, y, z = map(Dummy, 'xyz')
        k = expand_mul(self.equation(x, y, z))
        a, b, c = [k.coeff(i) for i in (x, y, z)]
        d = k.xreplace({x: pt.args[0], y: pt.args[1], z:0})
        t = -d/(a**2 + b**2 + c**2)
        if isinstance(pt, Point):
            return Point3D(pt.x + t*a, pt.y + t*b, t*c)
        if isinstance(pt, Point3D):
            return Point3D(pt.x + t*a, pt.y + t*b, pt.z + t*c)

def _set_function(f, self):
    from sympy.core.function import expand_mul
    if not self:
        return S.EmptySet
    if not isinstance(f.expr, Expr):
        return
    if self.size == 1:
        return FiniteSet(f(self[0]))
    if f is S.IdentityFunction:
        return self

    x = f.variables[0]
    expr = f.expr
    # handle f that is linear in f's variable
    if x not in expr.free_symbols or x in expr.diff(x).free_symbols:
        return
    if self.start.is_finite:
        F = f(self.step*x + self.start)  # for i in range(len(self))
    else:
        F = f(-self.step*x + self[-1])
    F = expand_mul(F)
    if F != expr:
        return imageset(x, F, Range(self.size))

def intersecting_product(a, b):
    """
    Returns the intersecting product of given sequences.

    The indices of each argument, considered as bit strings, correspond to
    subsets of a finite set.

    The intersecting product of given sequences is the sequence which
    contains the sum of products of the elements of the given sequences
    grouped by the *bitwise-AND* of the corresponding indices.

    The sequence is automatically padded to the right with zeros, as the
    definition of subset based on bitmasks (indices) requires the size of
    sequence to be a power of 2.

    Parameters
    ==========

    a, b : iterables
        The sequences for which intersecting product is to be obtained.

    Examples
    ========

    >>> from sympy import symbols, S, I, intersecting_product
    >>> u, v, x, y, z = symbols('u v x y z')

    >>> intersecting_product([u, v], [x, y])
    [u*x + u*y + v*x, v*y]
    >>> intersecting_product([u, v, x], [y, z])
    [u*y + u*z + v*y + x*y + x*z, v*z, 0, 0]

    >>> intersecting_product([1, S(2)/3], [3, 4 + 5*I])
    [9 + 5*I, 8/3 + 10*I/3]
    >>> intersecting_product([1, 3, S(5)/7], [7, 8])
    [327/7, 24, 0, 0]

    References
    ==========

    .. [1] https://people.csail.mit.edu/rrw/presentations/subset-conv.pdf

    """

    if not a or not b:
        return []

    a, b = a[:], b[:]
    n = max(len(a), len(b))

    if n&(n - 1): # not a power of 2
        n = 2**n.bit_length()

    # padding with zeros
    a += [S.Zero]*(n - len(a))
    b += [S.Zero]*(n - len(b))

    a, b = mobius_transform(a, subset=False), mobius_transform(b, subset=False)
    a = [expand_mul(x*y) for x, y in zip(a, b)]
    a = inverse_mobius_transform(a, subset=False)

    return a

    def eval(cls, arg):
        from sympy.simplify.simplify import signsimp
        from sympy.core.function import expand_mul

        if hasattr(arg, '_eval_Abs'):
            obj = arg._eval_Abs()
            if obj is not None:
                return obj
        if not isinstance(arg, Expr):
            raise TypeError("Bad argument type for Abs(): %s" % type(arg))
        # handle what we can
        arg = signsimp(arg, evaluate=False)
        if arg.is_Mul:
            known = []
            unk = []
            for t in arg.args:
                tnew = cls(t)
                if isinstance(tnew, cls):
                    unk.append(tnew.args[0])
                else:
                    known.append(tnew)
            known = Mul(*known)
            unk = cls(Mul(*unk), evaluate=False) if unk else S.One
            return known*unk
        if arg is S.NaN:
            return S.NaN
        if arg is S.ComplexInfinity:
            return S.Infinity
        if arg.is_Pow:
            base, exponent = arg.as_base_exp()
            if base.is_real:
                if exponent.is_integer:
                    if exponent.is_even:
                        return arg
                    if base is S.NegativeOne:
                        return S.One
                    if isinstance(base, cls) and exponent is S.NegativeOne:
                        return arg
                    return Abs(base)**exponent
                if base.is_nonnegative:
                    return base**re(exponent)
                if base.is_negative:
                    return (-base)**re(exponent)*exp(-S.Pi*im(exponent))
                return
            elif not base.has(Symbol): # complex base
                # express base**exponent as exp(exponent*log(base))
                a, b = log(base).as_real_imag()
                z = a + I*b
                return exp(re(exponent*z))

        if isinstance(arg, exp):
            return exp(re(arg.args[0]))
        if isinstance(arg, AppliedUndef):
            return
        if arg.is_Add and arg.has(S.Infinity, S.NegativeInfinity):
            if any(a.is_infinite for a in arg.as_real_imag()):
                return S.Infinity
        if arg.is_zero:
            return S.Zero
        if arg.is_nonnegative:
            return arg
        if arg.is_nonpositive:
            return -arg
        if arg.is_imaginary:
            arg2 = -S.ImaginaryUnit * arg
            if arg2.is_nonnegative:
                return arg2
        # reject result if all new conjugates are just wrappers around
        # an expression that was already in the arg
        conj = signsimp(arg.conjugate(), evaluate=False)
        new_conj = conj.atoms(conjugate) - arg.atoms(conjugate)
        if new_conj and all(arg.has(i.args[0]) for i in new_conj):
            return
        if arg != conj and arg != -conj:
            ignore = arg.atoms(Abs)
            abs_free_arg = arg.xreplace({i: Dummy(real=True) for i in ignore})
            unk = [a for a in abs_free_arg.free_symbols if a.is_real is None]
            if not unk or not all(conj.has(conjugate(u)) for u in unk):
                return sqrt(expand_mul(arg*conj))