//
//  main.cpp
//  KR_PO_VMA
//
//  Created by Berkovich Pavel on 12/8/19.
//  Copyright © 2019 Berkovich Pavel. All rights reserved.
//



#include <iostream>
#include<time.h>
#include <cmath>
#include <iomanip>
#include <float.h>

using namespace std;
const int size =3;
const double eps = pow(10, -3);
int iteration = 1;
int m,l;
double fi(0);




double** Build_U_matrix(){
    double ** matrix=new double *[size];
    for (int i = 0; i< size; ++i)
        matrix[i] = new double[size]{0};
    for(int i = 0; i<size;++i){
        matrix[i][i]=1;
    }
    matrix[l][l] = cos(fi);
    matrix[m][m] = cos(fi);
    matrix[m][l] = sin(fi);
    matrix[l][m] = -sin(fi);
    return matrix;
}

void Max_ne_diag (double**matrix){
    double max = DBL_MIN;
    for(int i=0;i<size;++i){
        for (int j=0; j<size; ++j) {
            if (abs(max)<abs(matrix[i][j]) && (i<j)){
                max =abs( matrix[i][j]);
                l = i;
                m=j;
            }
        }
    }
}

void counting_fi (double**matrix){
    fi = 0.5*atan((2*matrix[l][m])/(matrix[l][l] - matrix[m][m]));
}


double scalar(double* V1,double* V2,int n)
{
    double rez(0);
    for(int i = 0; i<size;++i)
        rez+= V1[i]*V2[i];
    return rez;
}

void Print_Vector (double *arr){
    for (int i = 0; i < size; ++i){
        cout<<setw(7);
        printf("%.15f \n",arr[i]);
    }
    cout<<endl;
}

void Print_Matrix (double** arr){
    for (int i=0;i<size;++i){
        for(int j=0;j<size;j++){
            printf("%15.5f ",arr[i][j]);
        }
    cout<<endl;
    }
    cout<<endl;
}

void Norm_prisvaivanie (double** arr1,double** arr2){
    for (int i=0;i<size;++i){
        for(int j=0;j<size;j++){
            arr1[i][j] = arr2[i][j];
        }
    }
}

void Norm_prisvaivanie (double* arr1,double* arr2){
    for (int i=0;i<size;++i){
            arr1[i] = arr2[i];
        }
    }

void Transpouse_matrix (double **arr){
    for(int i = 0; i < size; ++i){
            for(int j = i; j < size; ++j){
                swap(arr[i][j] , arr[j][i]);
            }
        }
}

void matrix_multiplication(double **arr1, double **arr2, double ** answer){
    for (int i=0; i<size; ++i){
        for(int j=0;j<size;++j){
            for(int k = 0; k < size; ++k){
                 answer[i][j] += arr1[i][k] * arr2[k][j];
            }
        }
    }
}

double * matrix_multiplication(double **arr1, double *arr2){
    double *buff_vector = new double[size] {0};
    for (int i=0; i<size; ++i){
        for(int j=0;j<size;++j){
            buff_vector[i]+=arr1[i][j]*arr2[j];
        }
    }
return buff_vector;
}


double first_Norm(double*x){
    double max = DBL_MIN;
           for(int i =1;i<size;++i){
               (abs(x[i])> abs(max))?(max = x[i]):(max=max);
           }
           return abs(max) ;
       }

double mark_for_Test_1(double** matrix){
    double summ(0);
    for(int i=0;i<size;++i){
        for (int j=0; j<size; ++j) {
            if (i ==j )
                continue;
            summ +=pow(matrix[i][j],2);
        }
    }
    return summ;
}

int mark_for_aprior(double** matrix){
    double summ(0),p(0);
    int answ(0);
    summ =mark_for_Test_1(matrix);
    p = 1.0 - 2.0/(size*(size-1));
    answ = log(eps/summ)/log(p)+1;
    return answ;
}


void Normirovka (double* arr1,double k){
for (int i=0;i<size;++i){
        arr1[i] /=k;
    }
}

void Obnulenie(double** matrix){
    for(int i = 0;i<size;++i)
        for(int j = 0;j<size;++j)
            matrix[i][j]=0;
}

void Task_1(double** arr){
    iteration=0;
    double ** buf2=new double *[size];
    for (int i = 0; i< size; ++i)
        buf2[i] = new double[size]{0};
    double ** buf=new double *[size];
    for (int i = 0; i< size; ++i)
        buf[i] = new double[size]{0};
    while (mark_for_Test_1(arr)>eps){
        iteration++;
        Max_ne_diag(arr);
        counting_fi(arr);
        Norm_prisvaivanie(buf2, Build_U_matrix());
        cout<<"U matrix:"<<endl;
        Print_Matrix(buf2);
        Transpouse_matrix(buf2);
        Obnulenie(buf);
        matrix_multiplication(buf2, arr, buf);
        Obnulenie(arr);
        matrix_multiplication(buf,Build_U_matrix(), arr);
        cout<<"Next matrix:"<<endl;
        Print_Matrix(arr);
    }
}
 void print_all_I_need (double *buff_pred,double * buff){
     cout<<"Наши значения y(k),y(k+1):"<<endl;
     Print_Vector(buff_pred);
     cout<<"~~~~~~~~~~"<<endl;
     Print_Vector(buff);

 }


double Task_2_1_and_2(double **arr,double* y){
    iteration = 1;
    double lamd(0),pred_lambd(0);
    double * buff_pred = y;
    double * buff =matrix_multiplication(arr, y);
    print_all_I_need(buff_pred,buff);
    lamd =  scalar(buff, buff_pred, size)/scalar(buff_pred, buff_pred, size);
    cout<<lamd<<endl;

    while (abs(lamd-pred_lambd)>eps) {
        iteration ++;
        pred_lambd = lamd;
        double k = first_Norm(buff);
        Normirovka(buff, k);
        Norm_prisvaivanie(buff_pred, buff);
        buff = matrix_multiplication(arr, buff);
        lamd =  scalar(buff, buff_pred, size) / scalar(buff_pred, buff_pred, size) ;
        cout<<lamd<<endl;
    }
    return lamd;
}


void print_all_I_need (double * buff_ppred,double * buff_pred,double * buff){
    cout<<"Наши значения y(k-1),y(k),y(k+1):"<<endl;
    Print_Vector(buff_ppred);
    cout<<"~~~~~~~~~~"<<endl;
    Print_Vector(buff_pred);
    cout<<"~~~~~~~~~~"<<endl;
    Print_Vector(buff);
}

double Task_2_3(double **arr,double* y){
     iteration = 1;
    double lamd(0), pred_lambd(0);
    double * buff_ppred = y;
    double * buff =matrix_multiplication(arr, y);
    double * buff_pred =new double[size];
    Norm_prisvaivanie(buff_pred, buff);
    buff =matrix_multiplication(arr, buff_pred);
    print_all_I_need(buff_ppred,buff_pred,buff);
    lamd =  sqrt(scalar(buff, buff_pred, size)/scalar(buff_pred, buff_ppred, size));
    cout<<lamd<<endl;
       while (abs(lamd-pred_lambd)>eps) {
           iteration ++;
           pred_lambd = lamd;
           //double k = first_Norm(buff);
           //Normirovka(buff, k);
           //Norm_prisvaivanie(buff_pred, buff);
           buff_ppred = buff_pred;
           buff_pred = buff;
           buff = matrix_multiplication(arr, buff);
           print_all_I_need(buff_ppred,buff_pred,buff);
           lamd =  sqrt(scalar(buff, buff_pred, size)/scalar(buff_pred, buff_ppred, size));
           cout<<lamd<<endl;
       }
       return lamd;
}


int main() {
    srand(time(NULL));
    double ** matrix=new double *[size];
    for (int i = 0; i< size; ++i)
         matrix[i] = new double[size];
    //Generate_matrix_2_0(matrix);
    //Getting_the_right_matrix(matrix);
    for (int i = 0; i< size; ++i)
        for (int j =0;j<size;++j){
            cin>>matrix[i][j];
        }
    cout<<"Матрица c которой будем работать:"<<endl;
    Print_Matrix(matrix);
    //Start
    double ** buf=new double *[size];
    for (int i = 0; i< size; ++i)
        buf[i] = new double[size]{0};
    Norm_prisvaivanie(buf, matrix);
    Task_1(buf);
    cout<<"Количество итераций для вычисления = "<<iteration<<"."<<endl;
    //Task_2
    double * y = new double [size]{0};
    y[0] = 1;
    int choise;
    cout<<"У вас случай 1 и 2 или случай 3:"<<endl;
    cin>>choise;
    switch (choise) {
        case 1:
        {
            double p = Task_2_1_and_2(matrix, y);
            printf("Максимальное по модулю собственное значение, полученное итерационно-степенным методом %10.15f.\n",p);
            break;
        }
        case 2:
        {
            double p = Task_2_1_and_2(matrix, y);
            printf("Максимальное по модулю собственное значение, полученное итерационно-степенным методом %10.15f.\n",p);
            break;

        }
        case 3:
        {
            double p = Task_2_3(matrix, y);
            printf("Максимальное по модулю собственное значение, полученное итерационно-степенным методом %10.15f.\n",p);
            break;
        }
        default:
            break;
    }
}