#include <iostream>
#include <math.h>
#include <cmath>
#include <iomanip>
using namespace std;

//=======================================================================
//wektor funkcji

void function_vector(double *f_vector, double x, double y, double z){

    f_vector[0] = x*x + y*y + z*z - 2;
    f_vector[1] = x*x + y*y - 1;
    f_vector[2] = x*x - y;

}


//=======================================================================
//macierz Jakobiego


void Jacobi_Matrix(double **matrix, double x, double y, double z){

    matrix[0][0] = 0;
    matrix[0][1] = 1/(2+(4*x*y));
    matrix[0][2] = y/(x+(2*x*y));
    matrix[1][0] = 0;
    matrix[1][1] = 1/(1+(2*y));
    matrix[1][2] = 1/(-1-(2*y));
    matrix[2][0] = (1+(2*y))/((2*z)+(4*y*z));
    matrix[2][1] = (1+(2*y))/((-2*z)-(4*y*z));
    matrix[2][2] = 0;

}


//=======================================================================


double max3(double x, double y, double z)
{
    if (x > y)
    {
        if (x > z)
            return x;
        return z;
    }
    else if (y > z)
        return y;
    return z;
}


int main()
{
    std::cout << std::setprecision(10);
    double tab[3], TOLX=1.0e-12, TOLF=1.0e-12, estX=1, estF=1, LastWynik[3];
    int max_iterations = 100, j=0;
    double **matrix = new double *[3];
    for(int i=0; i<3; i++) matrix[i] = new double [3];

    double wynik[3]={1.0,1.0,1.0};


//=====================================================
//petla liczaca kolejne przyblizenia trojki (x, y, z)


    while ((estX > TOLX || estF > TOLF) && j < max_iterations)
    {
    Jacobi_Matrix(matrix, wynik[0], wynik[1], wynik[2]);
    function_vector(tab, wynik[0], wynik[1], wynik[2]);

        for(int i=0; i<3; i++)
        {
            LastWynik[i]=wynik[i];
            wynik[i] = wynik[i] - (tab[0]*matrix[i][0] + tab[1]*matrix[i][1] + tab[2]*matrix[i][2]);
            cout<<"wynik["<<i<<"] = "<<wynik[i]<<endl;

        }

        cout<<"Estymatory :\nestX = "<<estX<<"\nestF = "<<estF<<endl;
        cout<<"\n=====================================\n\n";
        estX = fabs(max3(wynik[0], wynik[1], wynik[2]) - max3(LastWynik[0], LastWynik[1], LastWynik[2]));
        estF = fabs(max3(tab[0], tab[1], tab[2]));

        j++;
    }
    cout<<"Liczba iteracji: "<<j<<endl;

        return 0;
}