#include <stdio.h>
#include <math.h>

double fdet ( float d, float c, int n );
int fau ( float d, float c, int n, double ua[] );
void fld ( float d, float c, int n, double ld[] );
void fv ( int n, double vld[] );

void main()
{
 int i, i1, j, j1, n = 5 ;
 float d = 2., c = 1., det;
 double  ld[ 10 ], ua[100], vld[100];

 printf("\n     НАЧАЛЬНЫЕ ДАННЫЕ  \n");
 printf("       d = %5.2f  c = %5.2f  n = %5d  \n", d, c, n );

 det = fdet ( d, c, n );
 printf("\n   Детерминант   \n" );
 printf("  det = %5.2f  ", det );

 fau ( d, c, n, ua );
  printf("\n   Обратная матрица   \n" );
 for ( i=0; i < n; i++)
 {
  for ( j=0; j < n; j++)
  {
   i1 = i + 1;
   j1 = j + 1;
   printf("ua(%1d,%1d)=%2.4f ", i1, j1, ua[j + i*n] );
  }
  printf("\n " );
 }

 fld (d, c, n, ld );
 printf("\n   Собственные  значения   \n" );
 for ( i=0; i < n; i++)
 {
  i1 = i + 1;
  printf("  ld(%1d) = %2.2f  ", i1, ld[i] );
 }

 fv ( n, vld );
 printf("\n   Собственные  векторы   \n" );
 for ( i=0; i < n; i++)
 {
  for ( j=0; j < n; j++)
  {
   i1 = i + 1;
   j1 = j + 1;
   printf("  vld(%1d,%1d) = %5.2f  ", i1, j1, vld[j + i*n] );
  }
  printf("\n " );
 }
} // end_of_main

 double fdet ( float d, float c, int n )
// Вычисление определителя
// симметричной якобиевой матрицы
// с одинаковыми диагональными элементами.
//
// Параметры программы:
// d - диагональный элемент
// c - внедиагональный элемент
// n - порядок матрицы
// функция возвращает значение определителя
//
 {
  int i;
  double det, eps = 1.e-10, fi, p, r, r1, r2, r3, s1, s2;
  if ( fabs (d*d - 4*c*c) <= eps ) // совпадающие корни
                     // характеристического многочлена
  {
   for ( i=0, r = 1.; i < n; i++)
   r *= d/2.;
   det = r*( 1.+ n );
   return det;
  }
  if ( d*d > 4*c*c) // различные вещественные корни
                    // характеристического многочлена
  {
   r1 = sqrt( d*d - 4.*c*c );
   r2 = ( d + r1 );
   r3 = ( d - r1 );
   for ( i=0, s1 = 1., s2 = 1., p = 1.; i < n + 1; i++)
   {
    s1 *= r2;
    s2 *= r3;
    p *= 2.;
   }
   det = ( s1 - s2 )/( p*r1);
   return det;
  }
  if ( d*d < 4*c*c)  // мнимые корни
                     // характеристического многочлена
  for ( i=0, r = 1.; i < n; i++)
  r *= c;
  r1 = d/( 2 * c );
  fi = acos( r1 );
  det = r * sin( ( n + 1.)*fi)/sin(fi);
  return det;
 } // end_of fdet

 int fau ( float d, float c, int n, double ua[] )
// Вычисление обратной матрицы для
// симметричной якобиевой матрицы
// с одинаковыми диагональными элементами.
//
// Параметры программы:
// d - диагональный элемент
// c - внедиагональный элемент
// n - порядок матрицы
// ua - массив элементов обратной матрицы
// Функция возвращает 1,если матрица вырождена и
// 0 в противном случае.
 {
  int i, j, k, l ;
  double c1, c2, q1, q2, eps = 1.e-10, fi, r, r1, r2, r3, s1, s2, xjk;
  double sfi1, sfi2, sfi3;
  if ( fabs( d*d - 4*c*c) <= eps )  // совпадающие корни
                     // характеристического многочлена
  {
   q1 = -d/(2.*c);
   for ( j = 1 ; j < n + 1; j++ )
   for ( k = 1 ; k < n + 1; k++ )
   {
    if ( k + 1 >= j )
    {
     for ( l = 0, r = 1.; l < k - j + 1; l++ )
     r *= q1;
    }
    if ( k + 1 < j )
    {
     for ( l = 0, r = 1.; l < j - k - 1; l++ )
     r /= q1;
    }
    r1 = c*( j*q1*q1 - j + n + 1 );
    xjk = r/r1;
    if ( k < j )
    xjk *= ( j - n - 1 )*k;
    else
    xjk *= ( k - n - 1 )*j;
    ua[k-1+n*(j-1)] = xjk;
   }
   return 0;
  }
  if ( d*d > 4*c*c )    // различные вещественные корни
                        // характеристического многочлена
  {
   r1 = sqrt( d*d - 4.*c*c );
   q1 = ( -d + r1 )/(2.*c);
   q2 = ( -d - r1 )/(2.*c);

   for ( l = 0, r1 = 1., r2 = 1.; l < n + 1 ; l++ )
   {
    r1 *= q1;
    r2 *= q2;
   }

   c1 = -1./( c*( q1 - q2 ) * ( r1 - r2 ));

   for ( j = 1 ; j < n + 1; j++ )
   {
    for ( l = 0, r1 = 1., r2 = 1.; l < n + 1 - j ; l++ )
    {
     r1 *= q1;
     r2 *= q2;
    }
    c2 = c1*( r1 - r2 );

    for ( k = 1 ; k < n + 1; k++ )
    {
     for ( l = 0, r1 = 1., r2 = 1.; l < k ; l++ )
     {
      r1 *= q1;
      r2 *= q2;
     }
     xjk = c2*( r1 - r2 );
     ua[ k - 1 + n*(j - 1) ] = xjk;
     if ( k > j )
     {
      for ( l = 0, r1 = 1., r2 = 1.; l < k - j ; l++ )
      {
       r1 *= q1;
       r2 *= q2;
      }
      xjk += ( r1 - r2 )/( c * ( q1 - q2 ) );
      ua[ k - 1 + n*(j - 1) ] = xjk;
     }
    }
   }
   return 0;
  }

  if ( d*d < 4*c*c )   // мнимые корни
                       // характеристического многочлена
  {
   r1 = -d/ ( 2 * c );
   fi = acos( r1 );
   for ( j = 1 ; j < n + 1; j++ )
    {
    for ( k = 1 ; k < n + 1; k++ )
    {
     r1 = sin( fi );
     r2 = sin( (n + 1)*fi );
     r1 *= c*r2;
     if ( fabs(r1) <= eps ) return 1;
     r2 = sin(( n + 1 - j )*fi);
     r3 = sin( k*fi );
     xjk = -r2*r3/r1;
     if ( k > j  )
     {
      r1 = sin( (k - j )*fi );
      r2 = c*sin( fi );
      if ( fabs(r2) <= eps ) return 1;
      xjk += r1/r2;
     }
     ua[k - 1 + n*(j - 1)] = xjk;
    }
   }
  }
  return 0;
 }

 void fld ( float d, float c, int n, double ld[] )
// Вычисление собственных значений для
// симметричной якобиевой матрицы
// с одинаковыми диагональными элементами.
//
// Параметры программы:
// d - диагональный элемент
// c - внедиагональный элемент
// n - порядок матрицы
// ld - массив собственных значений
//
 {
  int i;
  float k ;
  double a, pi, pin1;
  pi = atan(1.)*4.;
  pin1=pi/(n+1);
  if(c>0.)
  {
   for ( k=1; k <= n; k++)
   {
   ld[k-1] = d - 2*c*cos(pin1*k);
   } return;
  }
   for ( k=n; k >= 1; k--)
   {
   ld[n-k] = d - 2*c*cos(pin1*k);
   }return;
 }

 void fv ( int n, double vld[] )
// Вычисление собственных векторов
// симметричной якобиевой матрицы
// с одинаковыми диагональными элементами.
//
// Параметры программы:
// n - порядок матрицы
// vld - массив собственных векторов
 {
      int j, k;
      double  a, fik, fik_j, a1, pi;
      {
      pi = atan(1.)*4.;
      for(k=1; k<=n; k++)
        {
         vld[(k-1)*n+0]=1.;
         fik = pi*k/(n + 1.);
         a1 = sin(fik);
          for(j=2; j<=n; j++)
           {
           fik_j = fik*j;
           a = sin(fik_j);
           vld[(k-1)*n+(j-1)] = pow(-1,j+1)*a/a1;
           }
        }
      }
 }