"""Original exact rank-deficient pseudoinverse identities with actual NumPy comparisons."""
from fractions import Fraction as F
import math,json
import numpy as np

def transpose(a):return [list(r) for r in zip(*a)]
def mul(a,b):return [[sum(F(x)*y for x,y in zip(row,col)) for col in zip(*b)] for row in a]
def mv(a,b):return [sum(F(x)*y for x,y in zip(row,b)) for row in a]
def validate(a,b,rcond):
 if type(a) is not list or not 1<=len(a)<=4 or any(type(row) is not list or not 1<=len(row)<=4 for row in a) or len({len(row) for row in a})!=1 or any(type(v) not in [int,F] or abs(v)>100 or F(v).denominator>10**12 for row in a for v in row) or type(b) is not list or len(b)!=len(a) or any(type(v) not in [int,F] or abs(v)>100 for v in b) or type(rcond) not in [float,int] or not math.isfinite(rcond) or not 0<=rcond<=1:raise ValueError('bounded exact rectangular matrix/RHS and finite cutoff0..1')

def identities(a,p):
 assert mul(mul(a,p),a)==a and mul(mul(p,a),p)==p and transpose(mul(a,p))==mul(a,p) and transpose(mul(p,a))==mul(p,a)
 return {'A_P_A_equals_A':True,'P_A_P_equals_P':True,'A_P_symmetric':True,'P_A_symmetric':True}

def compare_numpy(a,b,p,rcond):
 validate(a,b,rcond);aa=np.array(a,dtype=float);bb=np.array(b,dtype=float);pp=np.linalg.pinv(aa,rcond=rcond);x,res,rank,s=np.linalg.lstsq(aa,bb,rcond=rcond);xp=pp@bb;expected=np.array(mv(p,b),dtype=float)
 assert np.allclose(pp,np.array(p,dtype=float),rtol=1e-12,atol=1e-12) and np.allclose(x,expected,rtol=1e-12,atol=1e-12) and np.allclose(xp,expected,rtol=1e-12,atol=1e-12)
 return {'rcond':rcond,'pinv':pp.tolist(),'pinv_b':xp.tolist(),'lstsq_x':x.tolist(),'reported_rank':int(rank),'singular_values':s.tolist(),'lstsq_residual_array':res.tolist(),'residual_array_empty_is_not_zero_residual_claim':True}

def demo():
 u=[F(1),F(2)];v=[F(1),F(1)];a=[[x*y for y in v] for x in u];den=sum(x*x for x in u)*sum(y*y for y in v);p=[[y*x/den for x in u] for y in v];b=[F(1),F(3)];ids=identities(a,p);x=mv(p,b);assert x==[F(7,10)]*2
 rows=[]
 for k in range(-8,9):
  t=F(k,4);xx=[x[0]+t,x[1]-t];res=[z-y for z,y in zip(mv(a,xx),b)];sse=sum(z*z for z in res);norm=sum(z*z for z in xx);assert sse==F(1,5) and norm==F(49,50)+2*t*t;rows.append({'t':str(t),'x':[str(z) for z in xx],'residual_squared':str(sse),'solution_norm_squared':str(norm)})
 actual=compare_numpy(a,b,p,1e-12);under_a=[[F(1),F(1),F(0)],[F(0),F(0),F(1)]];under_p=[[F(1,2),F(0)],[F(1,2),F(0)],[F(0),F(1)]];under=compare_numpy(under_a,[F(2),F(3)],under_p,1e-12);identities(under_a,under_p)
 cutoff=[]
 diagonal=np.diag([1.0,1e-8])
 for rc in [1e-6,1e-10]:
  pp=np.linalg.pinv(diagonal,rcond=rc);cutoff.append({'rcond':rc,'pinv':pp.tolist()})
 assert cutoff[0]['pinv']==[[1.0,0.0],[0.0,0.0]] and cutoff[1]['pinv']==[[1.0,0.0],[0.0,1e8]]
 zeros=[[F(0),F(0)],[F(0),F(0)]];identities(zeros,zeros);zero=compare_numpy(zeros,[F(1),F(2)],zeros,1e-12)
 rejected=[]
 for label,aa,bb,rc in [('shape_b',a,[1],1e-12),('ragged',[[1],[1,2]],[1,2],1e-12),('NaN',[[float('nan')]],[1],1e-12),('bool',[[True]],[1],1e-12),('empty',[],[],1e-12),('negative_cutoff',a,b,-1),('infinite_cutoff',a,b,float('inf')),('overbudget',[[0]*5],[1],1e-12)]:
  try:validate(aa,bb,rc)
  except ValueError:rejected.append(label)
  else:raise AssertionError('expected rejection')
 encode=lambda mat:[[str(z) for z in row] for row in mat]
 return {'A':encode(a),'b':[str(z) for z in b],'rank_one_formula_denominator':str(den),'exact_pseudoinverse':encode(p),'four_exact_MP_identities':ids,'minimum_norm_x':[str(z) for z in x],'all17_family_checks':rows,'actual_NumPy_comparison':actual,'underdetermined_example':{'A':encode(under_a),'exact_pinv':encode(under_p),'actual':under},'cutoff_diagonal_1_and_1e_minus8':cutoff,'zero_matrix_actual':zero,'invalid_cases_rejected':rejected,'numpy_version':np.__version__,'scope':'own exact rank-one/underdetermined finite examples versus actual existing NumPy; numeric cutoff is not true-rank discovery, residual-array empty is not zero SSE, no normal-equation inverse/general numerical stability/model accuracy/Opera claim'}
if __name__=='__main__':print(json.dumps(demo(),ensure_ascii=True,indent=2))
