from numpy import *

def qrfact(A):
    m,n=shape(A)
    Q=matrix(zeros((m,n))); R=matrix(zeros((n,n)))
    for k in range(n):
        R[:k,k] = Q[:,:k].T*A[:,k]         
        Q[:,k] = A[:,k] - Q[:,:k]*R[:k,k] 
        R[k,k] = linalg.norm(Q[:,k])
        Q[:,k] = Q[:,k]/R[k,k]
    return Q,R