#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "task3.h"

double *inv(const int N, const double *M) {
    int i,j, k;
    double xil, xjl;
    double sinfi, cosfi, somebodythatiusedtoknow;
    double *invm;
    double *adjointm;
//    double *Tij;
    invm = (double*) malloc(N*N*sizeof(double));
    adjointm = (double*) malloc(2*N*N*sizeof(double));
//    Tij = (double*) malloc(N*N*sizeof(double));
    for (i=0; i<N; ++i){
        for (j=0; j<N; ++j){
            adjointm[i*2*N+j] = M[i*N+j];
        }
        j=N;
        for (j=N; j<N+i; ++j){
            adjointm[i*2*N+j] = 0;
        }
        adjointm[i*2*N+N+i] = 1;
        j = N+i+1;
        for (j=N+i+1; j<2*N; ++j){
            adjointm[i*2*N+j] = 0;
        }
    }
    i=0;
    j=0;
    for (i=0; i <N-1; ++i){
        for (j=i+1; j<N; ++j){
            //находим матрицу Tij
            if(((adjointm[i*2*N+i])<0) || ((adjointm[i*2*N+i])>0) || ((adjointm[j*2*N+i])<0) || ((adjointm[j*2*N+i])>0)){
            cosfi = adjointm[i*2*N+i]/(sqrt((adjointm[i*2*N+i])*(adjointm[i*2*N+i])+(adjointm[j*2*N+i])*(adjointm[j*2*N+i])));
            sinfi = adjointm[j*2*N+i]/(sqrt((adjointm[i*2*N+i])*(adjointm[i*2*N+i])+(adjointm[j*2*N+i])*(adjointm[j*2*N+i])));
/*            printf("cosfi = %f, sinfi = %f\n",cosfi, sinfi);*/
//            Tij[i*N+i] = M[i*N+i]/sqrt((M[i*N+i])*(M[i*N+i])+(M[j*N+i])*(M[j*N+i]));
//            printf("%f ", adjointm[i*2*N+j]);
//            Tij[j*N+j] = M[i*N+i]/sqrt((M[i*N+i])*(M[i*N+i])+(M[j*N+i])*(M[j*N+i]));
//            Tij[i*N+j] = -M[j*N+i]/sqrt((M[i*N+i])*(M[i*N+i])+(M[j*N+i])*(M[j*N+i]));
//            Tij[j*N+i] = M[j*N+i]/sqrt((M[i*N+i])*(M[i*N+i])+(M[j*N+i])*(M[j*N+i]));
//            for(k=0; k<N; ++k){
//                for(l=0; l<N; ++l){
//                    if (((k*N+l)!=(i*N+i)) && ((k*N+l)!=(i*N+j)) && ((k*N+l)!=(j*N+j)) && ((k*N+l)!=(j*N+i))){
//                        if(k==l){
//                            Tij[k*N+l] = 1;
//                        }
//                        else{
//                            Tij[k*N+l] = 0;
//                        }
//                    }
//                }
//            }
//            k=0;
//            l=0;
            //приводим M к треугольному виду
            for(k=0; k<(2*N); ++k){
                    xil = adjointm[i*2*N+k];
                    xjl = adjointm[j*2*N+k];
                    adjointm[i*2*N+k] = xil*cosfi+xjl*sinfi;
                    adjointm[j*2*N+k] = -xil*sinfi+xjl*cosfi;
            }
            k=0;
            }
        }
    }
/*    for (i=0; i<N; ++i){
        for (j=0; j<2*N; ++j){
            printf("%f ", adjointm[i*2*N+j]);
        }
        printf("\n");
    }
    printf("\n");*/
    i=0;
    j=0;
    k=0;
    //ищем обратную
    if((adjointm[(N-1)*2*N+N-1]<=0) && (adjointm[(N-1)*2*N+N-1]>=0)){
        for (i=0; i<N; ++i){
            for (j=0; j<N; ++j){
                invm[i*N+j] = 0;
//                printf("%f ", adjointm[i*2*N+j]);
            }
        }
        return invm;
    }
    for (i=N-2; i>=0; --i){
        for (k=i; k>=0; --k){
            for (j=0; j!= 2*N; ++j){
                if((adjointm[k*2*N+i+1]>0) || (adjointm[k*2*N+i+1]<0)){
                adjointm[k*2*N+j]=adjointm[k*2*N+j]-adjointm[(i+1)*2*N+j]*adjointm[k*2*N+i+1]/adjointm[(i+1)*2*N+i+1];
                }
            }
        }
        /*for(k=0;k<N;++k){
            for(j=0;j<2*N;++j){
                printf("%f ", adjointm[i*2*N+j]);
            }
            printf("\n");
        }
        printf("\n\n");*/
    }
    for(i=0;i<N;++i){
        if((adjointm[i*2*N+i]<1)||(adjointm[i*2*N+i]>1)){
                 somebodythatiusedtoknow = adjointm[i*2*N+i];
                 for(k=i;k!=2*N;++k){
                     adjointm[i*2*N+k]=adjointm[i*2*N+k]/somebodythatiusedtoknow;
                 }
        }
    }
    i=0;
    j=0;
    for (i=0; i<N; ++i){
        for (j=0; j<N; ++j){
            invm[i*N+j] = adjointm[i*2*N+j+N];
        }
    }
    return invm;
}