fix Stewart's algorithm

[araujoms]

Fixes #65

As described in Stewart's paper, one needs a diagonal matrix of signs to make the QR decomposition unique, and the resulting distribution uniform. This PR adds it, in a custom type that keeps the complexity of the algorithm at O(d^2), as opposed to the O(d^3) we would get with a naïve matrix conversion.

The resulting eigenvalue distribution looks pretty uniform:

using LinearAlgebra
using RandomMatrices
using Plots

function distribution(m, N)
    theta = Vector{Float64}(undef, m * N)
    for i = 1:N
        u = Stewart(ComplexF64, m)
        ev = eigvals(Matrix(u))
        theta[((i-1)*m+1):(i*m)] .= angle.(ev)
    end
    histogram(theta; normalize = true)
end

distribution(51, 10000)

[dlfivefifty]

Cool. We can drop Julia v1.0 support if you want

one thing bothers me: every orthogonal matrix is a product of reflections. Doesn’t that mean we can still represent the orthogonal matrix as an QRPackedQ ? We just need to figure out what the right reflections are… though maybe this is hard in practice

[araujoms]

Thanks, but the error I'm getting is for AdjointQ, which was introduced in 1.10. I think dropping 1.9 would be a bit extreme.

You're right, of course, but the problem is that computing the proper reflections is not numerically stable, which is why pretty much every implementation of QR in existence doesn't do it, but instead computes a stable set of reflections that makes a mess of the signs in the diagonal of the R factor.

[araujoms]

Ok, now I got something working from 1.6 to 1.10. It's not beautiful. I dropped 1.0 because it doesn't support Diagonal, and I really need it.

[araujoms]

I only contributed this PR because you asked for it, I'd appreciate if you would review it before I forget everything.