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

/* matrix_and_vector.c
   CSC 111 - Fall 2023 - Assignment 4

   Function implementations of matrix
   and vecotor operations based on the
   matrix_and_vector.h header file.

   Elvin Shoolbraid, V01029195 - 2023-10-28
*/

void set_vector(int n, double V[n], double s){
    
    for(int i = 0; i < n; i++){
        V[i] = s;
    }
    
}

void mul_vector_by_scalar(int n, double Vin[n], double Vout[n], double s){
    
    for(int i = 0; i < n; i++){
        Vout[i] = s*Vin[i];
    }    
    
}

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

double norm(int n, double V[n]){
    
    return sqrt(dot_product(n, V, V));
    
}

double angle_between_vectors(int n, double V1[n], double V2[n]){
    
    if (norm(n, V1) == 0 || norm(n, V2) == 0){
        return -1;
    }
    
    return acos((dot_product(n, V1, V2))/(norm(n, V1)*norm(n, V2)));
    
}

void add_vectors(int n, double V1[n], double V2[n], double Vout[n]){
    
    for(int i = 0; i < n; i++){
        Vout[i] = V1[i] + V2[i];
    }   
    
}

/****************************** MATRIX FUNCTIONS ******************************/

void identity(int n, double M[n][n]){
    
    for(int i = 0; i < n; i++){
        for(int j = 0; j < n; j++){
            if(i == j){
                M[i][j] = 1;
            } else {
                M[i][j] = 0;
            }
        }   
    }   
    
}

int matrices_equal(int rows, int cols, double A[rows][cols], double B[rows][cols]){
    
    for(int i = 0; i < rows; i++){
        for(int j = 0; j < cols; j++){
           
            if(A[i][j] != B[i][j]){
                return 0;
            }
            
        }   
    }   
    
    return 1;
    
}

void add_matrices(int rows, int cols, double A[rows][cols], double B[rows][cols], double C[rows][cols]){
    
    for(int i = 0; i < rows; i++){
        for(int j = 0; j < cols; j++){
           
            C[i][j] = A[i][j] + B[i][j];
            
        }   
    }   
    
}

double trace(int n, double A[n][n]){
    
    double sum = 0;
    for(int i = 0; i < n; i++){
        for(int j = 0; j < n; j++){
            if(i == j){
                sum += A[i][j];
            } 
        }   
    }   
    
    return sum;
    
}

void transpose(int n, int k, double A[n][k], double T[k][n]){
    
    for(int i = 0; i < n; i++){
        for(int j = 0; j < k; j++){
            
            T[j][i] = A[i][j];
            
        }   
    }   
}

void omit_row(int n, int k, double A[n][k], double B[n-1][k], int omit_idx){
    
    for(int i = 0; i < omit_idx; i++){
        
            for(int j = 0; j < k; j++){
                 
                B[i][j] = A[i][j];
                
            }
    }     
    
    for(int i = omit_idx; i < n-1; i++){
        
            for(int j = 0; j < k; j++){
                 
                B[i][j] = A[i+1][j];
                
            }
    }     
}


void omit_column(int n, int k, double A[n][k], double B[n][k-1], int omit_idx){
    
    for(int j = 0; j < omit_idx; j++){
        
            for(int i = 0; i < n; i++){
                 
                B[i][j] = A[i][j];
                
            }
    }     
    
    for(int j = omit_idx; j < k-1; j++){
        
            for(int i = 0; i < n; i++){
                 
                B[i][j] = A[i][j+1];
                
            }
    }     
}

void outer_product(int n, int m, double V1[n], double V2[m], double M[n][m]){
    
    for(int i = 0; i < n; i++){
        for(int j = 0; j < m; j++){
            
            M[i][j] = V1[i]*V2[j];
            
        }   
    }   
    
}

void circulant(int n, double A[n][n], double V[n]){
    
    for(int i = 0; i < n; i++){
        for(int j = 0; j < n; j++){
            
            if(j+i < n){
            
                A[i][j] = V[j+i];
                
            } else {
                A[i][j] = V[j+i-n];
            }
            
        }   
    }   
    
}

void matrix_vector_multiply(int n, int k, double A[n][k], double V[k], double Vout[n]){
    
    for(int i = 0; i < n; i++){
    
        double Ai[k];
        
        for(int j = 0; j < k; j++){
            Ai[j] = A[i][j];
        }
        
        Vout[i] = dot_product(k, Ai, V);
    
    }    
    
}

void matrix_multiply(int m, int n, int k, double A[m][n], double B[n][k], double C[m][k]){
    
    for(int i = 0; i < m; i++){
        for(int j = 0; j < k; j++){
            
        double Ai[n];
        double Bj[n];
        
        for(int k = 0; k < n; k++){
            Ai[k] = A[i][k];
        }
            
        for(int k = 0; k < n; k++){
            Bj[k] = B[k][j];
        }
        
        C[i][j] = dot_product(n, Ai, Bj);
            
        }
    
    }    
    
}