r"""
    Computes the quadratic sum

    .. math::

        V = \sum_{j=0}^{\infty} A^j B A^{j'}

    V is computed by solving the corresponding discrete lyapunov
    equation using the doubling algorithm.  See the documentation of
    `util.solve_discrete_lyapunov` for more information.

    Parameters
    ----------
    A : array_like(float, ndim=2)
        An n x n matrix as described above.  We assume in order for
        convergence that the eigenvalues of :math:`A` have moduli bounded by
        unity
    B : array_like(float, ndim=2)
        An n x n matrix as described above.  We assume in order for
        convergence that the eigenvalues of :math:`A` have moduli bounded by
        unity
    max_it : scalar(int), optional(default=50)
        The maximum number of iterations

    Returns
    ========
    gamma1: array_like(float, ndim=2)
        Represents the value :math:`V`

def m_quadratic_sum(A, B, max_it=50):
    r"""
    Computes the quadratic sum

    .. math::

        V = \sum_{j=0}^{\infty} A^j B A^{j'}

    V is computed by solving the corresponding discrete lyapunov
    equation using the doubling algorithm.  See the documentation of
    `util.solve_discrete_lyapunov` for more information.

    Parameters
    ----------
    A : array_like(float, ndim=2)
        An n x n matrix as described above.  We assume in order for
        convergence that the eigenvalues of :math:`A` have moduli bounded by
        unity
    B : array_like(float, ndim=2)
        An n x n matrix as described above.  We assume in order for
        convergence that the eigenvalues of :math:`A` have moduli bounded by
        unity
    max_it : scalar(int), optional(default=50)
        The maximum number of iterations

    Returns
    ========
    gamma1: array_like(float, ndim=2)
        Represents the value :math:`V`

    """

    gamma1 = solve_discrete_lyapunov(A, B, max_it)

    return gamma1