ENH: Added discrete Lyapunov transformation solve

The existing direct solution method to the discrete Lyapunov equation `solve_discrete_lyapunov` is not good for large matrices since if the original matrix is MxM, the current method requires solving a linear system of order M^2.

This request adds a solution method that uses a bilinear transformation from the discrete to the continuous Lyapunov equation, which can then be solved using the existing continuous solver which can take advantage of the LAPACK *trsyl routines.

It allows selecting between the two methods ('direct', 'bilinear'), and defaults to using the direct method when the input matrix dimension is less than 10 and the bilinear method otherwise.
This is now duplicated between here and the two private functions. I suggest to expand this docstring with both references and notes on math, and simplify the private ones as much as possible.