#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define tiny 1e-12
#include "real.h"

int rotate(int p, int q, real**A, int n, real** V);

void make_unity_matrix(real** V,int n){
	for(int col=0;col<n;col++){ V[col][col]=1;
		for(int row=col+1;row<n;row++) V[col][row]=V[row][col]=0; } }

int jacobi(real** A, int n, real** V){
// Jacobi diagonalization.
// Upper triangle of A is destroyed.
// V accumulates eigenvectors.

	make_unity_matrix(V,n);

	int max_rotations=5*n*n;
	int rotated, rotations=0;
	do{
		rotated=0;
		for(int row=0;row<n;row++)for(int col=row+1;col<n;col++){
			if( tiny*fabs(A[row][row])<fabs(A[col][row])
					|| tiny*fabs(A[col][col])<fabs(A[col][row])){
				rotated=1; 
				rotations++;
				if(rotations>max_rotations){
					fprintf(stderr,"max_rotations reached: %d\n",rotations);
					return -rotations;}
				rotate(row,col,A,n,V);
			}
		}
	}while(rotated==1);
	return rotations;
}