import numpy as np

# Calculates the QR decomposition of A using Givens rotations.
def QRGivens(A):
    m, n = A.shape
    R = A.copy()
    Q = np.eye(m)

    for j in range(n):
        for i in range(m-1, j, -1):
            a, b = R[i-1, j], R[i, j]
            G = np.array([[a, b], [-b, a]]) / np.sqrt(a**2 + b**2)
            R[i-1:i+1, :] = G @ R[i-1:i+1, :]
            Q[i-1:i+1, :] = G @ Q[i-1:i+1, :]

    return Q.T, R

if __name__ == "__main__":
    A = np.array([[2, 1, 1], [-1, -1, -3], [2, 3, 3], [1, 1, -1]]).astype(float)
    Q, R = QRGivens(A)

    print("Approximating example matrix with QR decomposition gives error: {}".format(np.linalg.norm(Q @ R - A)))

    # To generate random matricies one can use functions from the np.random
    # module.

    print("Approximation errors for some randomly generated matricies:")
    for i in range(50):
        M = np.random.standard_normal((i, i))
        Q, R = QRGivens(M)
        print(np.linalg.norm(Q @ R - M))