//
//  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;
            }
        }
    }
    cout<<"max ne diag = "<<matrix[l][m]<<". Where l = "<<l<<" and m = "<<m<<endl;
}

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 =0;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;
}

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 sin_cos(double** arr){
    double Plm = 2*arr[l][m]/(arr[l][l]-arr[m][m]);
    printf("P(l,m) =  2 * %f / (%f - %f) = %f\n",arr[l][m],arr[l][l],arr[m][m],2*arr[l][m]/(arr[l][l]-arr[m][m]));
    double _cos =pow(0.5*(1+1/sqrt(pow(Plm,2)+1)),0.5);
    printf("cos = sqrt(0,5 * (1 + 1/ sqrt((%f)^2+1))) = %f\n",Plm,_cos);
     double _sin = Plm/abs(Plm) * pow(0.5*(1-1/sqrt(pow(Plm,2)+1)),0.5);
    printf("sin = %f *sqrt(0,5 * (1 - 1/ sqrt((%f)^2+1))) = %f\n", Plm/abs(Plm),Plm,_sin);
}

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++;
        printf("~~~~ итерация № %d ~~~~\n",iteration);
        Max_ne_diag(arr);
        sin_cos(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);
        cout<<"~~~~~~~~~~"<<endl;
    }
}
 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;
    printf("~~~~ итерация № %d ~~~~\n",iteration);
    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 ++;
        printf("~~~~ итерация № %d ~~~~\n",iteration);
        pred_lambd = lamd;
        bool bool_norm;
        cout<<"Будем делать нормировку? {0,1}"<<endl;
        cin>>bool_norm;
        if (bool_norm){
            double k = first_Norm(buff);
                   Normirovka(buff, k);
        }
        Norm_prisvaivanie(buff_pred, buff);
        buff = matrix_multiplication(arr, buff);
        print_all_I_need(buff_pred,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;
    printf("~~~~ итерация № %d ~~~~\n",iteration);
    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 ++;
           printf("~~~~ итерация № %d ~~~~\n",iteration);
           pred_lambd = lamd;
           bool bool_norm;
           cout<<"Будем делать нормировку? {0,1}"<<endl;
           cin>>bool_norm;
           if (bool_norm){
               double k = first_Norm(buff_pred);
                      Normirovka(buff_pred, k);
           }
           Norm_prisvaivanie(buff_ppred, buff_pred);
           buff_pred =matrix_multiplication(arr, buff_ppred);
           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;
       }
       return lamd;
}


int main() {
    srand(time(NULL));
    double ** matrix=new double *[size];
    for (int i = 0; i< size; ++i)
         matrix[i] = new double[size];
    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;
    }
}