def pr(A):
	try:
		m,n=len(A[0]),len(A)
		for i in range(m):
			print '| ',
			for j in range(n):
				if type(A[j][i])==complex: print A[j][i],		
				else: print "%.5f"%A[j][i],
			print ' |'
	except:
		n=len(A)
		print '| ',
		for j in range(n):
			if type(A[j])==complex: print A[j],		
			else: print "%.5f"%A[j],
		print ' |'

		
def norm(A):
  try:
    m,n=len(A[0]),len(A)
    l=[]
    for j in range(n):
      a=0
      for i in range(m):
        a+=abs(A[j][i])
      l.append(a)
    return max(l)
  except:
    try:
      l=map(abs,A)
      f=lambda x,y:x+y
      return reduce(f,l)
    except:
      return abs(A)

def sum(A,B):
   ma,na,mb,nb=len(A[0]),len(A),len(B[0]),len(B)
   C=[]
   if (ma-mb)**2+(na-nb)**2:
      print 'Operation imposible'
   else:
      for i in range(na):
         c=[]
         for j in range(ma):
            c.append(A[i][j]+B[i][j])
         C.append(c)
   return C

def sub(A,B):
   ma,na,mb,nb=len(A[0]),len(A),len(B[0]),len(B)
   C=[]
   if (ma-mb)**2+(na-nb)**2:
      print 'Operation imposible'
   else:
      for i in range(na):
         c=[]
         for j in range(ma):
            c.append(A[i][j]-B[i][j])
         C.append(c)
   return C


def mult(A,B):
   m,n=len(B[0]),len(A)
   C=[]
   if m!=n:
    print 'Operation imposible'
   else:
      for i in range(len(B)):
         c=[]
         for j in range(len(A[0])):
            t=0
            for k in range(len(A)):
               t+=A[k][j]*B[i][k] 
            c.append(t)
         C.append(c)
   return C 


def kA(k,A):
  C=[]
  if type(A[0])==list:
    for i in range(len(A)):
      c=[]
      for j in range(len(A[0])):
        c.append(k*A[i][j])
      C.append(c)
    return C
  else:
    for j in range(len(A)):   
      C.append(k*A[j])
    return C


def tran(A):
   m,n=len(A[0]),len(A)
   C=[]
   for i in range(m):
      c=[]
      for j in range(n):
         c.append(A[j][i])
      C.append(c)
   print
   return M(C)


def mnr(x,y,A):
   m,n=len(A[0]),len(A)
   if x in range(m+1) and y in range(n+1):
      C=[]
      for j in range(y-1)+range(y,n):
         c=[]
         for i in range(x-1)+range(x,m):
            c.append(A[j][i])
         C.append(c)
      return C
   else:
      print 'Operation imposible'


def inv(A):
   m,n=len(A[0]),len(A)
   if m!=n or det(A)==0:
      print 'Operation imposible'
   elif m==1 and n==1:
      return [[1./A[0][0]]]
   else:
      k,C=1./det(A),[]
      for j in range(n):
         c=[]
         for i in range(m):
            a=((-1)**(i+j))*det(mnr(i+1,j+1,A))
            c.append(a)
         C.append(c)       
      return tran(kA(k,C))  


def det(A):
   m,n=len(A[0]),len(A)
   if m!=n:
    print 'Operation imposible'
   else:
    if m==1:
     d=A[0][0]
    else:
     d=0    
     for k in range(m):
       C=mnr(k+1,1,A)
       d+=((-1)**k)*A[0][k]*det(C)
   return d

 
def replace(x,y,A):
  m,n=len(A[0]),len(A)
  if x in range(1,m+1) and y in range(1,m+1) and x!=y:
   C=[]
   for j in range(n):
      c=[]
      for i in range(m):
       if i!=x-1 and i!=y-1:
         c.append(A[j][i])
       if i==x-1:
         c.append(A[j][y-1])
       if i==y-1:
         c.append(A[j][x-1])
      C.append(c)
   return C
  elif x==y:
   return A
  else:
   print 'Operation imposible'
  

def mix(k,x,y,A):
   m,n=len(A[0]),len(A)
   C=[]
   if x in range(1,m+1) and y in range(1,m+1):
     for j in range(n):
       c=[]
       for i in range(m):
         if i!=y-1:
          c.append(A[j][i])
         else:
          c.append(A[j][y-1]-k*A[j][x-1])
       C.append(c)
     return C
   else:
     print 'Operation imposible'


def gd(k,m,n,A,p):
  M=abs(A[k-1][k-1])
  for j in range(k-1,n):
    for i in range(k-1,m):
      M=max(M,abs(A[j][i]))
      if M==abs(A[j][i]):
        I,J=i,j
  C=replace(k,I+1,A)
  C=replace(k,J+1,tran(C))
  u,v=p[k-1],p[J]
  p[k-1],p[J]=v,u
  return tran(C),p

def gt(m,n,A,p):
  N=min(m,n)
  C,p=gdc(1,m,n,A,p)
  for i in range(2,N+1):
    C,p=gdc(i,m,n,C,p)
  return C,p

def gdc(k,m,n,A,p):
  C,p=gd(k,m,n,A,p)
  if type(C[k-1][k-1])==complex and round(abs(C[k-1][k-1]),5)==0:
    C[k-1][k-1]=0
    return C,p
  elif type(C[k-1][k-1])!=complex and round(C[k-1][k-1],13)==0:
    C[k-1][k-1]=0
    return C,p
  else:
    for i in range(k,m):
      if type(C[k-1][i])!=complex:
        C=mix(float(C[k-1][i])/C[k-1][k-1],k,i+1,C)
      else:
        C=mix(C[k-1][i]/C[k-1][k-1],k,i+1,C)
    C=mix((C[k-1][k-1]-1.)/C[k-1][k-1],k,k,C)
    return C,p    


def mx(m,k,x,a,b):
  c=[]
  for i in range(m):
    c.append(b[i]-x*a[k-1][i])
  return c


def st(m,a,b):
  X=[]
  x=b[m-1]
  X.insert(0,x)
  for i in range(m-1):
    b=mx(m,m-i,X[0],a,b)
    x=b[m-2-i]
    X.insert(0,x)
  return X

def srt(a,p):
  n=len(p)
  A=[]
  for k in range(1,n+1):
    for j in range(n):
      if p[j]==k:
        A.append(a[j])
      else:
        pass
  return A

def cond(A):
  return norm(A)*norm(inv(A))

def solver(A,B):
  m,n,nr=len(A[0]),len(A),len(A)+1
  R=[]
  for j in range(n):
    R.append(A[j])
  R.append(B)
  p,pr=range(1,n+1),range(1,nr+1)
  OBJ=gt(m,n,R,p)
  At,p=OBJ[0],OBJ[1]
  Rt=gt(m,nr,R,pr)[0]
  ra=rr=0
  for i in range(min(m,n)):
    ra+=int(round(abs(At[i][i])))
  for i in range(min(m,nr)):
    rr+=int(round(abs(Rt[i][i])))
  flag=1
  while flag:
    if int(rr)>int(ra):
      print 'Rang A=',int(ra),'\nRang A|B=',int(rr)
      print 'No solution'
      flag=0
    else:
      a,b0,b=[],[],At[nr-1]
      for j in range(int(ra)):
        a.append(At[j])
        b0.append(0.)
      S=[]
      X=st(int(ra),a,b)
      for j in range(int(ra),n):
        X.append(0.)
      S.append(X)
      for k in range(int(ra)+1,n+1):
        b0=mx(int(ra),k,1,At,b0)
        X=st(int(ra),a,b0)
        for j in range(int(ra)+1,k):
          X.append(0.)
        X.append(1.)
        for j in range(k+1,n+1):
          X.append(0.)
        S.append(X)
        b0=[]
        for j in range(int(ra)):
          b0.append(0.)
      S=tran(srt(tran(S),p))
      print 'Rang A=',int(ra),'\nRang A|B=',int(rr)
      for i in range(n-int(ra)+1):
        print u'X%d=%ls'%(i,S[i])
      flag=0
      print
      return S

def eigenvv(A):
  if len(A)>4 or len(A)<2:
    print 'error'
    print 'The Order of Matrix must be  1<n<5 '
    return
  B,e,r,R=[],[],[],[]
  for k in range(len(A)):
    e0=[]
    for l in range(len(A)):
        if l==k: e0.append(1)
        else: e0.append(0)
    e.append(e0)
  e,a=M(e),M(A)
  lst='KV(a0,a1,a2)','TR(a0,a1,a2,a3)','KW()' 

  def SV(R):
    n,rr,b,a,S=len(A),0,[],[],[]
    p=range(1,n+1)
    OBJ=gt(n,n,R,p)
    At,p=OBJ[0],OBJ[1]
    for i in range(n):
      rr+=int(round(abs(At[i][i])))
    for j in range(rr):
      b.append(0)
      a.append(At[j])
    for k in range(rr+1,n+1):
      b=mx(rr,k,1,At,b)
      X=st(rr,a,b)
      for j in range(rr+1,k):
        X.append(0.)
      X.append(1.)
      for j in range(k+1,n+1):
        X.append(0.)
      S.append(X)
      b=[]
      for j in range(rr):
        b.append(0)
    S=tran(srt(tran(S),p))
    return S,rr

  def KV(a0,a1,a2): 
    d=(a1**2-4*a0*a2)
    if type(d)==complex or d<0:
        d+=0j     
    l1=(-a1-d**.5)/2.*a2
    l2=(-a1+d**.5)/2.*a2  
    return l1,l2

  def TR(a0,a1,a2,a3):
    V=[]
    t=(a2*a2-3*a3*a1)/(9*a3*a3)
    q=2*a2**3/(27*a3**3)-a2*a1/(3*a3*a3)+a0/a3
    u0=KV(t**3,q,1)[0]
    v0=KV(t**3,q,1)[1]
    if type(u0)==complex or u0<0: u0+=0j
    if type(v0)==complex or v0<0: v0+=0j
    u,v=u0**(1./3),v0**(1./3)
    e=complex(-.5,.5*3**.5)
    while round(abs(u*v-t),7)!=0:           u*=e
    for s in range(1,4):
      l0=e**(s-1)*u+e**(4-s)*v-a2/(3.*a3) 
      if type(l0)==complex:
        V.append(round(l0.real,7)+round(l0.imag,7)*1j)
      else:
        V.append(round(l0,7))
    return V

  def KW():
   x0=a3/4.
   r=x0**4-a3*x0**3+a2*x0**2-a1*x0+a0
   q=-4*x0**3+3*a3*x0**2-2*a2*x0+a1
   p=6*x0**2-3*a3*x0+a2  
   z=TR(-q*q,p*p-4*r,2*p,1.)[0]
   v1,v2=z**.5,(p+z)/2.
   v3=q/(2.*v1)
   return KV(v2+v3,-v1,1)[0]-x0,KV(v2+v3,-v1,1)[1]-x0,KV(v2-v3,v1,1)[0]-x0,KV(v2-v3,v1,1)[1]-x0
 
  for i in range(len(A)+1):
    b=[]; r.append(det(a-i*e))
    for j in range(len(A)+1):
        b.append(j**i)
    B.append(b)
  R.append(r)
  K=mult(inv(B),R)

  if len(A)==2:
    a0,a1,a2=K[0][0],K[0][1],K[0][2]
  elif len(A)==3:
    a0,a1,a2=K[0][0],K[0][1],K[0][2]
    a3=K[0][3]
  elif len(A)==4:
    a0,a1,a2=K[0][0],K[0][1],K[0][2]
    a3,a4=K[0][3],K[0][4]

  sval,svec,svec0=[],[],[]
  if a.x==(a[0][0]*e).x:
    for l in e.x:
      print
      print 'eigenvalue=',a[0][0]
      print 'eigenvector=',l
      sval.append(a[0][0])
      svec.append(l)
  else:
    for l in eval(lst[len(A)-2]):
      rr=SV((a-l*e).x)[1]
      print
      print 'eigenvalue=',l
      sval.append(l)
      for i in range(len(A)-rr):
        print 'eigenvector=',SV((a-l*e).x)[0][i];
        svec0.append(SV((a-l*e).x)[0][i])
      print
      svec.append(svec0[0])
      svec0=[]
  for pip in range(len(sval)):
      print "eigenvalue: ",
      print sval[pip]
      print "eigenvector: "
      pr(svec[pip])
  return sval,svec

class M:
  def __init__(self,x):
    self.x=x
    pr(x)
  def __getitem__(self,i):
    try:
      if self.x[i]==list:
        return M(self.x[i])
      else:
        return self.x[i]
    except:
      return self.x

  def __len__(self):
    try:
      return len(self.x)
    except:
      return 0

  def __add__(self,y):
    try:
      print
      pr(self.x)
      print '+'
      pr(y)
      print '='
      return M(sum(self.x,y))
    except:
      if type(y)==list:
        print 'operation cancel' 
      elif type(self.x)==list:
        print 'operation cancel'
      else:
        return self.x+y

  def __radd__(self,y):
    try:
      print
      pr(self.x)
      print '+'
      pr(y)
      print '='
      return M(sum(y,self.x))
    except:
      if type(y)==list:
        print 'operation cancel' 
      elif type(self.x)==list:
        print 'operation cancel'
      else:
        return self.x+y

  def __sub__(self,y):
    try:
      print
      pr(self.x)
      print '-'
      pr(y)
      print '='
      return M(sub(self.x,y))
    except:
      if type(y)==list:
        print 'operation cancel' 
      elif type(self.x)==list:
        print 'operation cancel'
      else:
        return self.x-y

  def __rsub__(self,y):
    try:
      print
      pr(self.x)
      print '-'
      pr(y)
      print '='
      return M(sub(y,self.x))
    except:
      if type(y)==list:
        print 'operation cancel' 
      elif type(self.x)==list:
        print 'operation cancel'
      else:
        return y-self.x

  def __mul__(self,y):
    try:
      print
      pr(self.x)
      print '*'
      pr(y)
      print '='
      return M(mult(self.x,y))
    except:
      if type(y)==list:
        print self.x
        return M(kA(self.x,y))
      elif type(self.x)==list:
        print y
        return M(kA(y,self.x))  
      else:
        return y*self.x

  def __rmul__(self,y):
    try:
      print
      pr(self.x)
      print '*'
      pr(y)
      print '='
      return M(mult(y,self.x))
    except:
      if type(y)==list:
        print self.x
        return M(kA(self.x,y))
      elif type(self.x)==list:
        print y
        return M(kA(y,self.x))  
      else:
        return y*self.x

  def __pow__(self,n):
    a=self.x
    try:
      print
      for i in range(n-1):
        a=mult(a,self.x)
      return M(a)
    except:
      return a**n

  def __setitem__(self,j,b):
    try:
      self.x[j]=b
    except:
        print 'operation cancel'      
  def __abs__(self):
    return norm(self.x)  

  def __repr__(self):
    try:
      return '%ls'%self.x
    except:
      return str(self.x)

