# -*- coding: utf-8 -*-
from __future__ import division
import math

def rotate(p, q, A, V):
	if q < p:
		p, q = q, p
	else:
		pass 
	n = len(A)
	app = A[p][p]
	aqq = A[q][q]
	apq = A[q][p]
	phi = 0.5*math.atan2(2*apq,(aqq-app))
	c = math.cos(phi)
	s = math.sin(phi)
	A[p][p] = c**2*app + s**2*aqq - 2*s*c*apq
	A[q][q] = s**2*app + c**2*aqq + 2*s*c*apq
	A[q][p] = 0

	for i in range(p): 
		aip = A[p][i]
		aiq = A[q][i]
		A[p][i] = c*aip-s*aiq	
		A[q][i] = c*aiq+s*aip	
	for i in range(p+1, q): 
		api = A[i][p]
		aiq = A[q][i]
		A[i][p] = c*api-s*aiq	
		A[q][i] = c*aiq+s*api
	for i in range(q+1, n):
		api = A[i][p]
		aqi = A[i][q]
		A[i][p] = c*api-s*aqi	
		A[i][q] = c*aqi+s*api

	for i in range(n):
		vip = V[p][i]
		viq = V[q][i]
		V[p][i] = c*vip-s*viq
		V[q][i] = c*viq+s*vip

def jacobi(M):
	n = len(M)
	m = len(M[0])
	if n != m:
		return 'not square matrix'
	else:
		pass
	V = [ [ 0 if i != j else 1 for i in range(n)] for j in range(n) ]
	A = M
	eps = 1e-12
	rotated = 1
	sweeps = 0
	while rotated == 1:
		rotated = 0
		for r in range(n):
			for c in range(r+1, n):
				if abs(A[c][r]) > eps*(abs(A[c][c])+abs(A[r][r])):
					rotated = 1
					rotate(r, c, A, V)
					sweeps += 1
	E = [A[i][i] for i in range(n)]
	return E, V, sweeps