"""Owned two-column classical Gram-Schmidt plus triangular least-squares solve."""
import math,json
def dot(x,y):return sum(a*b for a,b in zip(x,y))
def qr(A):
    if len(A)<2 or any(len(r)!=2 for r in A):raise ValueError('m-by-2 only')
    if not all(math.isfinite(x) for r in A for x in r):raise ValueError('finite only')
    a,b=zip(*A);r11=math.sqrt(dot(a,a))
    if r11<=1e-12:raise ValueError('first column has no usable norm')
    q1=[x/r11 for x in a];r12=dot(q1,b);v=[x-r12*q for x,q in zip(b,q1)];r22=math.sqrt(dot(v,v))
    if r22<=1e-12:raise ValueError('dependent columns in this subset')
    return [[u,x/r22] for u,x in zip(q1,v)],[[r11,r12],[0.,r22]]
def demo():
    A=[[1.,1.],[1.,0.],[0.,1.]];b=[1.,2.,3.];Q,R=qr(A)
    qcols=list(zip(*Q));qtb=[dot(q,b) for q in qcols]
    x2=qtb[1]/R[1][1];x1=(qtb[0]-R[0][1]*x2)/R[0][0];x=[x1,x2]
    recon=[[sum(Q[i][k]*R[k][j] for k in range(2)) for j in range(2)] for i in range(3)]
    gram=[[dot(a,c) for c in qcols] for a in qcols];pred=[dot(row,x) for row in A];res=[z-y for z,y in zip(b,pred)];atr=[dot(c,res) for c in zip(*A)]
    assert max(abs(gram[i][j]-(i==j)) for i in range(2) for j in range(2))<1e-12
    assert max(abs(recon[i][j]-A[i][j]) for i in range(3) for j in range(2))<1e-12
    assert abs(x1-2/3)<1e-12 and abs(x2-5/3)<1e-12 and max(map(abs,atr))<1e-12
    rejected=[]
    for name,B in [('dependent',[[1.,2.],[2.,4.]]),('zero_first_column',[[0.,1.],[0.,2.]])]:
        try:qr(B)
        except ValueError:rejected.append(name)
        else:raise AssertionError(name)
    return {'A':A,'b':b,'reduced_Q':Q,'upper_R':R,'Q_transpose_Q':gram,'QR_reconstruction':recon,'Q_transpose_b':qtb,'least_squares_x':x,'predictions':pred,'residual_b_minus_Ax':res,'A_transpose_residual':atr,'residual_squared_norm':dot(res,res),'own_subset_rejected':rejected,'algorithm':'two-column classical Gram-Schmidt, not LAPACK','NumPy_executed':False}
if __name__=='__main__':print(json.dumps(demo(),indent=2))
