#pragma once
#include <iostream>
#include <vector>
#include <stdio.h>

void extended_euclid(long a, long b, long *x, long *y, long *d) {/* calculates a * *x + b * *y = gcd(a, b) = *d */
    long q, r, x1, x2, y1, y2;

    if (b == 0) {
        *d = a, *x = 1, *y = 0;
        return;
    }
    x2 = 1, x1 = 0, y2 = 0, y1 = 1;

    while (b > 0) {
        q = a / b, r = a - q * b;
        *x = x2 - q * x1, *y = y2 - q * y1;
        a = b, b = r;
        x2 = x1, x1 = *x, y2 = y1, y1 = *y;
    }
    *d = a, *x = x2, *y = y2;
}


long inverse(long a, long n) { /* computes the inverse of a modulo n */
    long d, x, y;
    extended_euclid(a, n, &x, &y, &d);
    if (d == 1) return x;
    return 0;
}



class ZpnVector {
public:
    int p, n;
    std::vector <int> X;
    ZpnVector() : n(1), p(1) {
        X.clear();
        X.push_back(0);
    };

    ZpnVector(std::vector <int> _v, int _p) : X(_v), p(_p) {
        n = _v.size();
        for each (int x in X) {
            x = x % p;
        }
    };

    ZpnVector(const ZpnVector& _v) {
        p = _v.p;
        n = _v.n;
        X.clear();
        for (int i = 0; i < _v.n; i++) {
            X.push_back(_v.X[i]);
        }
    };

    ZpnVector& operator= (const ZpnVector& _v) {
        X = _v.X; n = _v.n; p = _v.p;
        return *this;
    };

    ZpnVector& operator+ (const ZpnVector& _v) const {
        std::vector<int> cache;
        for (int i = 0; i < n; i++) {
            cache[i] = X[i] + _v.X[i];
        }
        return ZpnVector(cache, p);
    };

    ZpnVector& operator- (const ZpnVector& _v) const {
        std::vector<int> cache;
        for (int i = 0; i < n; i++) {
            cache.push_back(X[i] - _v.X[i]);
        }
        return ZpnVector(cache, p);
    };

    ZpnVector& operator* (const int& k) const {
        std::vector<int> cache;
        for (int i = 0; i < n; i++) {
            cache.push_back(X[i] * k);
        }
        return ZpnVector(cache, p);
    };

    ZpnVector& operator/ (const int& k) const {
        std::vector<int> cache;
        for (int i = 0; i < n; i++) {
            cache.push_back( X[i] * inverse(k, p) );
        }
        return ZpnVector(cache, p);
    };

    int dotProduct(const ZpnVector& _v) const { //Scalar
        int dP = 0;
        for (int i = 0; i < n; i++) {
            dP += X[i] * _v.X[i];
        }
        return (dP % p);
    };

    double magnitude() const {
        int m = 0;
        for (int i = 0; i < n; i++) {
            m += X[i] * X[i];
        }
        return sqrt(m);
    }
    ~ZpnVector() {};
};

void GaussianElimination(std::vector< ZpnVector > matrix) {
    int rows = matrix.size();
    int cols = matrix[0].X.size();

    if (rows > cols) { //если кол-во векторов больше чем их размерность
        std::cout << "Linear Independent!\n";
        return;
    }

    int curRow = 0; //строка
    int curCol = 0; //колонка

    while (curRow < rows) { // i - vector number
        int nonZeroRow = curRow;
        bool zeroCol = false;
        while (matrix[nonZeroRow].X[curCol] == 0) { //проходимся по curCol элементам векторов начиная с curRow-го
            nonZeroRow++;
            if (curRow == rows) {
                zeroCol = true;
                break;
            }
        }

        if (!zeroCol) { //нашли вектор с ненулевым curCol элементом
            if (nonZeroRow != curRow)
                std::swap(matrix[curRow], matrix[nonZeroRow]);
        }
        else {
            curCol++;
            continue;
        }

        ZpnVector cacheV = matrix[curRow] / matrix[curRow].X[curCol];
        matrix[curRow] = cacheV; //нормализуем вектор по curCol элемету <- PROBLEM!!!

        for (int j = curRow + 1; j < rows; j++) {
            matrix[j] = matrix[j] - matrix[curRow] * matrix[j].X[curCol];
        }
        curCol++;
        curRow++;
    }

    int zeroVectors = 0;
    for (int i = 0; i < rows; i++) {
        if (matrix[i].magnitude() == 0)
            zeroVectors++;
    }

    if (zeroVectors > 0 && rows - zeroVectors < cols)
        std::cout << "Linear Independent!\n";
}