#include <iostream>
#include <vector>
#include <cmath>
#include <stdexcept>
#include <iomanip>
#define Eps 0.00000001

bool istiBrojevi(double x, double y){
    return std::fabs(x - y) <=Eps * (std::fabs(x)+ std::fabs(y));
}

bool isteMatrice(std::vector<std::vector<double>> a, std::vector<std::vector<double>> b){
    for(int i=0; i<a.size(); i++){
        for(int j=0; j<a.at(0).size(); j++) if(!istiBrojevi(a.at(i).at(j), b.at(i).at(j))) return false;
    }
    return true;
}

std::vector<std::vector<double>> MnozenjeMatrica(std::vector<std::vector<double>> m1,std::vector<std::vector<double>> m2) {
    if(m1.size() == 0 && m2.size()==0) return std::vector<std::vector<double>>();
    if(m1.size()==0 or m2.size()==0) throw std::domain_error("Nije moguce mnozenje matrica!");
    for(int i=0; i<m1.size(); i++){
        if(m1.at(i).size() != m1.at(0).size()) throw std::domain_error("Neispravna matrica!");
    }
    for(int i=0; i<m2.size(); i++){
        if(m2.at(i).size() != m2.at(0).size()) throw std::domain_error("Neispravna matrica!");
    }

    if(m1.at(0).size() != m2.size()) throw std::domain_error("Nije moguce mnozenje matrica!");


    std::vector<std::vector<double>> m3(m1.size(), std::vector<double>(m2.at(0).size()));
    for(int i = 0; i < m1.size(); i++)
    for(int j = 0; j < m2.at(0).size(); j++) {
        double suma =0;
        for(int k = 0; k < m2.size(); k++)
        suma += m1.at(i).at(k) * m2.at(k).at(j);
        m3.at(i).at(j) = suma;
    }
    return m3;
}


std::vector<std::vector<std::vector<double>>> DoolitleLU(std::vector<std::vector<double>> H){
     std::vector<std::vector<std::vector<double>>> vrati;

     if(H.size()==0) return vrati;
     for(int i=0; i<H.size(); i++) if(H.at(i).size()!=H.at(0).size()) throw std::domain_error("Neispravna matrica!");
     if(H.size() != H.at(0).size()) throw std::domain_error("Matrica mora biti kvadratnog oblika!");

    vrati.resize(2);
    vrati.at(0).resize(H.size(), std::vector<double>(H.size()));
    vrati.at(1).resize(H.size(), std::vector<double>(H.size()));

    for (int i = 0; i < H.size(); i++) {
        for (int k = i; k < H.size(); k++) {
            double suma = 0;
            for (int j = 0; j < i; j++)
                suma+= (vrati.at(0).at(i).at(j) * vrati.at(1).at(j).at(k));
            vrati.at(1).at(i).at(k) = H.at(i).at(k) - suma;
        }
        for (int k = i; k < H.size(); k++) {
            if (i == k) vrati.at(0).at(i).at(i) = 1;
            else {
                double suma = 0;
                for (int j = 0; j < i; j++)
                    suma+= (vrati.at(0).at(k).at(j) * vrati.at(1).at(j).at(i));
                vrati.at(0).at(k).at(i) = (H.at(k).at(i) - suma) / vrati.at(1).at(i).at(i);
            }
        }
    }
    return vrati;

}
bool ProvjeriFaktorizaciju(std::vector<std::vector<double>> H, std::vector<std::vector<std::vector<double>>> K){

    if(H.size()==0 and K.size() == 0) return true;
    if(K.size()!=2) return false;
    for(int i=0; i<K.at(0).size(); i++) if(K.at(0).at(i).size() != K.at(0).at(0).size()) return false;
    for(int i=0; i<K.at(1).size(); i++) if(K.at(1).at(i).size() != K.at(1).at(0).size()) return false;
    for(int i=0; i<H.size(); i++) if(H.at(i).size()!=H.at(0).size()) return false;
    if(H.size()!=H.at(0).size()) return false;
    if(K.at(0).size() != K.at(0).at(0).size()) return false;
    if(K.at(1).size() != K.at(1).at(0).size()) return false;
    if(H.size()!=K.at(0).size() || H.size()!=K.at(1).size() || K.at(0).size()!=K.at(1).size()) return false;
    std::vector<std::vector<double>> m;
    try{
        m=MnozenjeMatrica(K.at(0), K.at(1));
        return isteMatrice(m, H);
    }
    catch (std::domain_error e){
        return false;
    }

    return true;


}

int main(){


    return 0;
}