"""Original exact synchronous Jacobi, with SPD divergence counterexample."""
from fractions import Fraction as F
import json

def run(a,b,start,cap=8,tol=F(1,10**8)):
 if type(a) is not list or not 1<=len(a)<=3 or any(type(r) is not list or len(r)!=len(a) for r in a) or any(type(x) not in [int,F] or abs(x)>100 for row in a for x in row) or any(type(v) is not list or len(v)!=len(a) or any(type(x) not in [int,F] or abs(x)>100 for x in v) for v in [b,start]) or type(cap) is not int or not 1<=cap<=16 or type(tol) not in [int,F] or tol<=0:raise ValueError('bounded exact square/vector and positive tolerance/cap')
 a=[list(map(F,r)) for r in a];b=list(map(F,b));x=list(map(F,start));n=len(a)
 if any(a[i][i]==0 for i in range(n)):raise ValueError('nonzero diagonal required')
 rows=[]
 for k in range(cap):
  old=x[:];x=[(b[i]-sum(a[i][j]*old[j] for j in range(n) if j!=i))/a[i][i] for i in range(n)];res=[b[i]-sum(a[i][j]*x[j] for j in range(n)) for i in range(n)];norm=max(map(abs,res));rows.append({'step':k+1,'old':[str(v) for v in old],'new':[str(v) for v in x],'residual_b_minus_Ax':[str(v) for v in res],'residual_infinity':str(norm)})
  if norm<=tol:return rows,'residual_tolerance'
 return rows,'budget_exhausted'

def demo():
 a=[[4,1],[2,3]];b=[6,8];rows,status=run(a,b,[0,0]);direct=[F(b[0]*a[1][1]-a[0][1]*b[1],10),F(a[0][0]*b[1]-b[0]*a[1][0],10)];assert direct==[1,2] and rows[0]['new']==['3/2','8/3']
 # For this B, B^2=(1/6)I; independent closed-form even iterates.
 for k in [2,4,6,8]:assert rows[k-1]['new']==[str(v*(1-F(1,6)**(k//2))) for v in direct]
 inplace=[F(0),F(0)]
 for i in range(2):inplace[i]=(F(b[i])-sum(F(a[i][j])*inplace[j] for j in range(2) if j!=i))/a[i][i]
 assert inplace==[F(3,2),F(5,3)] and inplace!=list(map(F,rows[0]['new']))
 spd=[[F(1) if i==j else F(9,10) for j in range(3)] for i in range(3)];rhs=[F(14,5)]*3;bad,bs=run(spd,rhs,[0,0,0]);assert bs=='budget_exhausted'
 for k,row in enumerate(bad,1):assert row['new']==[str(1-F(-9,5)**k)]*3
 eig=[]
 for vector,lam in [([1,1,1],F(14,5)),([1,-1,0],F(1,10)),([1,0,-1],F(1,10))]:
  actual=[sum(spd[i][j]*vector[j] for j in range(3)) for i in range(3)];assert actual==[lam*v for v in vector];eig.append({'vector':vector,'eigenvalue':str(lam)})
 iteration=[[(F(0) if i==j else -spd[i][j]/spd[i][i]) for j in range(3)] for i in range(3)];iteration_eigen=[]
 for vector,lam in [([1,1,1],F(-9,5)),([1,-1,0],F(9,10)),([1,0,-1],F(9,10))]:
  assert [sum(iteration[i][j]*vector[j] for j in range(3)) for i in range(3)]==[lam*v for v in vector];iteration_eigen.append({'vector':vector,'eigenvalue':str(lam)})
 assert max(abs(F(row['eigenvalue'])) for row in iteration_eigen)==F(9,5)
 rejected=[]
 for label,fn in [('zero_diagonal',lambda:run([[0]],[1],[0])),('ragged',lambda:run([[1,2]],[1],[0])),('nonfinite',lambda:run([[float('inf')]],[1],[0])),('bool_cap',lambda:run([[1]],[1],[0],True)),('negative_tol',lambda:run([[1]],[1],[0],tol=-1))]:
  try:fn()
  except ValueError:rejected.append(label)
  else:raise AssertionError('expected rejection')
 return {'A':a,'b':b,'direct_exact_solution':[str(v) for v in direct],'synchronous_eight_steps':rows,'status':status,'in_place_first_step_different':[str(v) for v in inplace],'SPD_does_not_guarantee_Jacobi_convergence':{'A':[[str(v) for v in row] for row in spd],'b':[str(v) for v in rhs],'three_independent_eigenvectors':eig,'iteration_matrix':[[str(v) for v in row] for row in iteration],'independent_iteration_eigenvectors':iteration_eigen,'iteration_spectral_radius':'9/5','iteration_on_equal_components_factor':'-9/5','steps':bad,'status':bs,'closed_form_checked_every_step':'s_k=1-(-9/5)^k'},'invalid_cases_rejected':rejected,'scope':'small exact synchronous old/new vector iteration, not CG/Gauss-Seidel or direct solver; SPD eigenvalues positive but iteration spectral radius9/5; no speed/performance or real model claim'}
if __name__=='__main__':print(json.dumps(demo(),ensure_ascii=True,indent=2))
