#include <perms.h>
int _imp_mainep(int _imp_argc, char **_imp_argv) {
  static const int Nmax = 4;
  auto void Testmatrix(double A, int N);
  auto void Householder(double A, double W, int N, int K);
  int I;
  double A[4 /*1:4*/][4 /*1:4*/];
  double W[4 /*1:4*/];
  Testmatrix(A, Nmax);
  Printstring(_imp_str_literal("\nTEST MATRIX EIGENVALUES\n"));
  Householder(A, W, Nmax, 1);
  for (I = 1; I <= Nmax; I++) {
    Printfl(W[I], 8);
    Newline();
  }
  Pprofile;
  exit(0);
  void Testmatrix(double A, int N) {
    int I;
    int J;
    double T;
    double C;
    double D;
    double F;
    C = N * (N + 1) * (2 * N - 5) / 6;
    D = 1 / C;
    A = -D;
    for (I = 1; I <= N - 1; I++) {
      F = I;
      A = D * F;
      A = A;
      A = D * (C - REXP(F, 2));
      for (J = 1; J <= I - 1; J++) {
        T = J;
        A = -D * F * T;
        A = A;
      }
    }
    for (I = 1; I <= N; I++) {
      for (J = 1; J <= N; J++) Printfl(A, 4);
      Newline();
    }
  }
  void Householder(double A, double W, int N, int K) {
    auto void Householdertridiagonalisation(double A, double B, double C,
                                            int N);
    auto void Tridibisection(double C, double B, double W, int N, double *Norm);
    auto void Tridiinverseiteration(double C, double B, double W, double Z,
                                    int N, double Norm);
    auto void Backtransformation(double A, double B, double Z, int N);
    double Z[N][N];
    double B[N];
    double C[N];
    int I;
    int J;
    double Norm;
    static const double Ten = 10;
    Householdertridiagonalisation(A, B, C, N);
    Tridibisection(C, B, W, N, Norm);
    if (K == 2) return;
    for (I = 1; I <= N; I++) {
      Printfl(B[I], 6);
      Printfl(C[I], 6);
      Printfl(W, 6);
      Newline();
    }
    Tridiinverseiteration(C, B, W, Z, N, Norm);
    Backtransformation(A, B, Z, N);
    for (I = 1; I <= N; I++)
      for (J = 1; J <= N; J++) A = Z[I][J];
    return;
    void Householdertridiagonalisation(double A, double B, double C, int N) {
      int I;
      int J;
      int K;
      double Ai;
      double Sigma;
      double H;
      double Bj;
      double Bigk;
      double Bi;
      double Q0[N - 1];
      for (I = N; I >= 3; I--) {
        Sigma = 0;
        for (K = 1; K <= I - 1; K++) Sigma += REXP(A, 2);
        Ai = A;
        if (Ai >= 0)
          Bi = -Sqrt(Sigma);
        else
          Bi = Sqrt(Sigma);
        B = Bi;
        Printfl(Bi, 5);
        Printfl(Sigma, 5);
        Newline();
        if (!Bi) continue;
        H = Sigma - Ai * Bi;
        A = Ai - Bi;
        for (J = I - 1; J >= 1; J--) {
          Bj = 0;
          for (K = I - 1; K >= J; K--) Bj += A * A;
          for (K = J - 1; K >= 1; K--) Bj += A * A;
          Q0[J] = Bj / H;
        }
        Bigk = 0;
        for (J = I - 1; J >= 1; J--) Bigk += A * Q0[J];
        Bigk = Bigk / (2 * H);
        for (J = I - 1; J >= 1; J--) Q0[J] = Q0[J] - Bigk * A;
        for (J = I - 1; J >= 1; J--)
          for (K = J; K >= 1; K--) A = A - A * Q0[K] - A * Q0[J];
      }
      for (I = N; I >= 1; I--) C = A;
      B = A;
      B = 0;
    }
    void Tridiinverseiteration(double C, double B, double W, double Z, int N,
                               double Norm) {
      int I;
      int J;
      double Bi;
      double Bi1;
      double Z1;
      double Lambda;
      double U;
      double S;
      double V;
      double H;
      double Eps;
      double Eta;
      double M[N];
      double P[N];
      double Q[N];
      double R[N];
      double Int[N];
      double X[N + 2];
      Lambda = Norm;
      Eps = Norm / (REXP(Ten, 11));
      for (J = 1; J <= N; J++) {
        Lambda -= Eps;
        if (W < Lambda) Lambda = W;
        U = C - Lambda;
        V = B;
        if (!V) V = Eps;
        for (I = 1; I <= N - 1; I++) {
          Bi = B;
          if (!Bi) Bi = Eps;
          Bi1 = B;
          if (!Bi1) Bi1 = Eps;
          if (Mod(U) > Mod(Bi)) {
            M[I + 1] = Bi / U;
            P[I] = U;
            Q[I] = V;
            R[I] = 0;
            U = C - Lambda - M[I + 1] * V;
            V = Bi1;
            Int[I + 1] = -1;
          } else {
            M[I + 1] = U / Bi;
            if (M[I + 1] == 0 && Bi <= Eps) M[I + 1] = 1;
            P[I] = Bi;
            Q[I] = C - Lambda;
            R[I] = Bi1;
            U = V - M[I + 1] * Q[I];
            V = -M[I + 1] * R[I];
            Int[I + 1] = 1;
          }
        }
        P[N] = U;
        Q[N] = 0;
        R[N] = 0;
        X[N + 1] = 0;
        X[N + 2] = 0;
        H = 0;
        Eta = 1 / N;
        for (I = N; I >= 1; I--) {
          U = Eta - Q[I] * X[I + 1] - R[I] * X[I + 2];
          if (P[I])
            X[I] = U / P[I];
          else
            X[I] = U / Eps;
          H += Mod(X[I]);
        }
        H = 1 / H;
        for (I = 1; I <= N; I++) X[I] = X[I] * H;
        for (I = 2; I <= N; I++)
          if (Int[I] <= 0)
            X[I] = X[I] - M[I] * X[I - 1];
          else {
            U = X[I - 1];
            X[I - 1] = X[I];
            X[I] = U - M[I] * X[I - 1];
          }
        H = 0;
        for (I = N; I >= 1; I--) {
          U = X[I] - Q[I] * X[I + 1] - R[I] * X[I + 2];
          if (!P[I])
            X[I] = U / Eps;
          else
            X[I] = U / P[I];
          H += REXP(X[I], 2);
        }
        H = 1 / Sqrt(H);
        for (I = 1; I <= N; I++) Z = X[I] * H;
      }
    }
    void Tridibisection(double C, double B, double W, int N, double *Norm) {
      double L;
      double G;
      double H;
      double Lambda;
      double P1;
      double Q1;
      double Y;
      int I;
      int J;
      int K;
      int A1;
      int A2;
      double P[N];
      P[1] = 0;
      *Norm = Mod(C) + Mod(B);
      for (I = 2; I <= N; I++) {
        L = Mod(B) + Mod(C) + Mod(B);
        if (L > *Norm) *Norm = 1;
      }
      for (I = 1; I <= N - 1; I++)
        if (!B)
          P[I + 1] = *Norm * *Norm / REXP(Ten, 23);
        else
          P[I + 1] = REXP(B, 2);
      for (K = 1; K <= N; K++) {
        G = *Norm;
        H = -*Norm;
        for (J = 1; J <= 39; J++) {
          Lambda = (G + H) / 2;
          P1 = 0;
          Q1 = 1;
          A1 = 0;
          for (I = 1; I <= N; I++) {
            Y = (C - Lambda) * Q1 - P[I] * P1;
            Printstring(_imp_str_literal("$"));
            Printfl(P1, 5);
            Printfl(Q1, 5);
            Printfl(Y, 5);
            Newline();
            P1 = Q1;
            Q1 = Y;
            if ((P1 >= 0 && Q1 >= 0) || (P1 < 0 && Q1 < 0)) {
              A1++;
              Printstring(_imp_str_literal("+"));
              Write(A1, 1);
              Printfl(P1, 5);
              Printfl(Q1, 5);
              Printfl(Lambda, 5);
              Newline();
            }
          }
          if (Q1 == 0 && P1 > 0) {
            A1--;
            if (K == 3 && J < 4) {
              Printstring(_imp_str_literal("-"));
              Write(A1, 1);
              Printfl(P1, 5);
              Printfl(Q1, 5);
              Newline();
            }
          }
          if (A1 >= K)
            H = Lambda;
          else
            G = Lambda;
          if (K == 3 && J < 4) {
            Printstring(_imp_str_literal("/"));
            Write(A1, 2);
            Write(K, 2);
            Printfl(H, 4);
            Printfl(G, 4);
            Printfl(Lambda, 4);
            Newline();
          }
        }
        W = (G + H) / 2;
      }
      return;
    }
    void Backtransformation(double A, double B, double Z, int N) {
      int I;
      int J;
      int K;
      double S;
      for (J = 1; J <= N; J++)
        for (K = 3; K <= N; K++) {
          if (!B) continue;
          S = 0;
          for (I = 1; I <= K - 1; I++) S += A * Z;
          S = S / (B * A);
          for (I = 1; I <= K - 1; I++) Z += S * A;
        }
    }
  }
  exit(0);
  return (1);
}