#include <stdio.h>
#include <math.h>
typedef double R;
void eigen(int n, R m[n][n], R l[n], R vc[n][n]){
R m2[n][n]; {int i=n; while(i--) {int j=n; while(j--) m2[i][j] = m[i][j];}}
{int i=n; while(i--) {int j=n; while(j--) vc[i][j] = i==j;}}
while(1){R mod = 0; int i=0, j=0;
{int k=n; while(k--){int m=n; while((--m)>k){R q = fabs(m2[k][m]);
  if(q > mod) {mod=q; i=k; j=m;}}}}
if(mod < 0.00000000001) break;
R th = 0.5*atan(2*m2[i][j]/(m2[i][i] - m2[j][j]));
{R c = cos(th), s = sin(th);
void twst(R m[n][n]){int k=n; while(k--){R t = c*m[i][k]+s*m[j][k];
   m[j][k] = -s*m[i][k]+c*m[j][k]; m[i][k] = t;}}
{int k=n; while(k--){R t = c*m2[k][i] + s*m2[k][j];
   m2[k][j] = -s*m2[k][i] + c*m2[k][j]; m2[k][i]=t;}}
twst(m2); twst(vc);}}
{int j=n; while(j--) l[j] = m2[j][j];}}

#define n 4
//R x[][n] = {{1, .01, .4, 0}, {.01, -1.6, 0, 0}, {.4, 0, 2, 2}, {0, 0, 2, -0.7}};
R x[][n] = {{1, .6, -1.1, -1.3}, {.6, 1, 2.1, .31}, {-1.1, 2.1, 2, 2}, {-1.3, .31, 2, 2}};
R vec[n][n], vl[n];
int main(){
{int i=n; while(i--) {int j=n; while(j--) if(x[i][j] != x[j][i]) printf("uns %d %d\n", i, j);}}
eigen(n, x, vl, vec);
int i; for(i=0; i<n; ++i){ printf("%11.7f: ", vl[i]);
   {int j; for(j=0; j<n; ++j) printf("%11.7f ", vec[i][j]); printf("\n");}}
{R me = 0; int j = n; while(j--){ // We examine the jth eigen vec and val.
    {int k = n; while(k--){ // We consider nth component of jth vec.
    R c = 0; {int y = n; while(y--) c += x[y][k]*vec[j][y];}
    {int e = fabs(vl[j]*vec[j][k] - c); if(e > me) me = e;}}}}
  printf("max error = %e\n", me);}
return 0;}