docs: 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