#include <math.h>
#include "real.h"

int rotate(int p, int q, real**A, int n, real** V){
// Jacobi rotation eliminating A_pq.
// Only upper triangle of A is updated.
// The matrix of eigenvectors V is also updated.
	if(q<p) {int tmp=p; p=q; q=tmp;}
	real app = A[p][p];
	real aqq = A[q][q];
	real apq = A[q][p];
	real phi=0.5*atan2(2*A[q][p],A[q][q]-A[p][p]);
	real c=cos(phi), s=sin(phi);
	A[p][p] = c * c * app + s * s * aqq - 2 * s * c * apq;
	A[q][q] = s * s * app + c * c * aqq + 2 * s * c * apq;
//	A[q][p] = ( c * c - s * s ) * apq + s * c * ( app - aqq );
	A[q][p]=0;

	for(int i=0;i<p;i++){
		double aip=A[p][i],aiq=A[q][i];
		A[p][i] = c*aip-s*aiq;
		A[q][i] = c*aiq+s*aip;
	}

	for(int i=p+1;i<q;i++){
		double api=A[i][p],aiq=A[q][i];
		A[i][p] = c*api-s*aiq;
		A[q][i] = c*aiq+s*api;
	}

	for(int i=q+1;i<n;i++){
		double api=A[i][p],aqi=A[i][q];
		A[i][p] = c*api-s*aqi;
		A[i][q] = c*aqi+s*api;
	}

	for(int i=0;i<n;i++){
		double vip=V[p][i],viq=V[q][i];
		V[p][i] = c*vip-s*viq;
		V[q][i] = c*viq+s*vip;
	}
	return 0;
}