"""Original exact LU partial pivoting with earlier L columns swapped."""
from fractions import Fraction as F
from itertools import permutations
import json

def matrix(a):
 if type(a) is not list or not 1<=len(a)<=4 or any(type(row) is not list or len(row)!=len(a) for row in a) or any(type(x) not in [int,F] or abs(x)>100 for row in a for x in row):raise ValueError('bounded exact square1..4')
 return [list(map(F,row)) for row in a]

def multiply(a,b):return [[sum((a[i][k]*b[k][j] for k in range(len(b))),F(0)) for j in range(len(b[0]))] for i in range(len(a))]

def det(a):
 # Independent Leibniz permutations, no elimination.
 n=len(a);total=F(0)
 for perm in permutations(range(n)):
  sign=(-1)**sum(perm[i]>perm[j] for i in range(n) for j in range(i+1,n));v=F(sign)
  for i,j in enumerate(perm):v*=a[i][j]
  total+=v
 return total

def factor(a,swap_previous_L=True):
 u=matrix(a);n=len(u);l=[[F(i==j) for j in range(n)] for i in range(n)];order=list(range(n));trace=[]
 for k in range(n):
  pivot=max(range(k,n),key=lambda i:(abs(u[i][k]),-i))
  if u[pivot][k]==0:raise ValueError('singular exact pivot')
  if pivot!=k:
   u[k],u[pivot]=u[pivot],u[k];order[k],order[pivot]=order[pivot],order[k]
   if swap_previous_L:
    for j in range(k):l[k][j],l[pivot][j]=l[pivot][j],l[k][j]
  eliminated=[]
  for i in range(k+1,n):
   l[i][k]=u[i][k]/u[k][k]
   for j in range(k,n):u[i][j]-=l[i][k]*u[k][j]
   eliminated.append({'row':i,'multiplier':str(l[i][k])})
  trace.append({'column':k,'pivot_row_before_swap':pivot,'row_order_after':order[:],'swapped_previous_L_columns':list(range(k)) if swap_previous_L and pivot!=k else [],'eliminated':eliminated})
 return l,u,order,trace

def solve(a,b):
 aa=matrix(a);n=len(aa)
 if type(b) is not list or len(b)!=n or any(type(x) not in [int,F] or abs(x)>100 for x in b):raise ValueError('bounded exact right-hand side')
 l,u,p,tr=factor(aa);pb=[F(b[i]) for i in p];y=[F(0)]*n;x=y[:]
 for i in range(n):y[i]=pb[i]-sum(l[i][j]*y[j] for j in range(i))
 for i in reversed(range(n)):x[i]=(y[i]-sum(u[i][j]*x[j] for j in range(i+1,n)))/u[i][i]
 assert multiply(l,u)==[aa[i] for i in p] and [sum(aa[i][j]*x[j] for j in range(n)) for i in range(n)]==list(map(F,b))
 return l,u,p,tr,pb,y,x

def demo():
 a=[[0,2,1],[1,1,1],[2,3,4]];b=[7,6,20];l,u,p,tr,pb,y,x=solve(a,b);assert x==[1,2,3] and [r['pivot_row_before_swap'] for r in tr]==[2,2,2]
 d=det(matrix(a));assert d==-3;cramer=[]
 swap_count=sum(row['pivot_row_before_swap']!=row['column'] for row in tr);upper_diagonal_product=u[0][0]*u[1][1]*u[2][2];assert swap_count==2 and (-1)**swap_count*upper_diagonal_product==d
 for j in range(3):
  c=[row[:] for row in a]
  for i in range(3):c[i][j]=b[i]
  cramer.append(det(matrix(c))/d)
 assert cramer==x
 wrong_l,wrong_u,wrong_p,_=factor(a,False);assert multiply(wrong_l,wrong_u)!=[matrix(a)[i] for i in wrong_p]
 inverse_order=[p.index(i) for i in range(3)];assert [multiply(l,u)[i] for i in inverse_order]==matrix(a)
 rejected=[]
 for label,bad in [('singular',[[1,2],[1,2]]),('ragged',[[1,2],[3]]),('empty',[]),('nonfinite',[[float('inf')]]),('bool',[[True]])]:
  try:factor(bad)
  except ValueError:rejected.append(label)
  else:raise AssertionError('expected rejection')
 enc=lambda a:[[str(x) for x in row] for row in a]
 return {'A':a,'b':b,'L':enc(l),'U':enc(u),'row_permutation_PA_order':p,'inverse_permutation_for_A_equals_PLU':inverse_order,'pivot_trace':tr,'Pb':[str(v) for v in pb],'forward_y':[str(v) for v in y],'backward_x':[str(v) for v in x],'independent_Leibniz_det':str(d),'swap_count':swap_count,'upper_diagonal_product':str(upper_diagonal_product),'determinant_permutation_sign_checked':True,'independent_Cramer_solution':[str(v) for v in cramer],'wrong_without_previous_L_swap':{'L':enc(wrong_l),'U':enc(wrong_u),'PA_equals_LU':False},'invalid_inputs_rejected':rejected,'scope':'own exact PA=LU; SciPy uses inverse permutation A=P L U; no SciPy execution/floating stability/sparse solver/benchmark'}
if __name__=='__main__':print(json.dumps(demo(),ensure_ascii=True,indent=2))