#include <stdio.h>
#include <math.h>
#include <iostream>
#include <string.h>
#include <iomanip>
#define epsilon 0.000001
#define n 3
#define debug

using namespace std;

long double przedzialy[3][3],J[3][3], F[3], F_transponowane[3], J_transponowane[3][3], X[3], J_F[3], X_nowe[3];


//___________________funkcje jakobianu:__________________________

long double Jacobian_x0(long double y, long double z)
{
    return (y*y*y*y + z*z);

}

long double Jacobian_x1(long double x)
{

    return (1.0);
}

long double Jacobian_x2(long double x)
{

    return (0.0);
}

long double Jacobian_y0(long double x,long double y)
{

    return  (4*x*y*y*y);

}

long double Jacobian_y1(long double y)
{
    if(y==-4.0)
    {
        cout<<"Dzielenie przez 0! -1/z+5. STOP"<<endl;
        exit(0);
    }
    else
    {
    return (-1/(y+4));
    }
}

long double Jacobian_y2(long double y)
{

    return (3*y*y);
}

long double Jacobian_z0(long double z,long double x)
{
    //cout<<"X to "<<x<<"Z to: "<<z<<endl;
    long double wynik = 2*x*z;
    return wynik;
}

long double Jacobian_z1(long double z)
{
    if(z==-5.0)
    {
        cout<<"Dzielenie przez 0! -1/z+5. STOP"<<endl;
        exit(0);
    }
    else
    {
    return (-1/(z+5));
    }
}

long double Jacobian_z2(long double z)
{

    return (-1.0);
}


//_____________________Koniec funkcji jakobianu_____________________


//_____________________Funkcje macierzy F_________________________
long double funkcja_0(long double x, long double y, long double z)
{

    return (x*(y*y*y*y+z*z)-5);
}

long double funkcja_1(long double x, long double y, long double z)
{/*
    if(y<=-4 || z<=-5)
    {
        cout<<"Ujemna liczba logarytmowana. STOP"<<endl;
        cout<<"Liczby to:"<<endl;
        cout<<"x    =      "<<x<<endl;
        cout<<"y    =      "<<y<<endl;
        cout<<"z    =      "<<z<<endl;
        exit(0);
    }
    else
    {*/
    return (x-log10(fabs(y+4))-log10(fabs(z+5)));
    //}
}

long double funkcja_2(long double x, long double y, long double z)
{

    return (y*y*y-z);
}

//_______________________Koniec funkcji F___________________________



//_______________________Odwracanie macierzy_______________________

long double odwracanie_macierzy(long double X[n][n], long double wynik[n][n])
{
    int i, j, k;
    long double wsp_1;
    long double tempor2;
    long double Temp[3][2*3];

    //wczytanie macierzy do pomocniczej tablicy

    for (i = 0; i < n; ++i)
    {
        for (j = 0; j < n; ++j)
        {
            Temp[i][j] = X[i][j];
        }
    }

    long double sarrus =Temp[0][0]*Temp[1][1]*Temp[2][2]+Temp[0][1]*Temp[1][2]*Temp[2][0]+Temp[0][2]*Temp[1][0]*Temp[2][1]-Temp[0][2]*Temp[1][1]*Temp[2][0]-Temp[0][1]*Temp[1][0]*Temp[2][2]-Temp[0][0]*Temp[1][2]*Temp[2][1];
    if (sarrus == 0)
    {
        printf("Macierzy nie mozna odwrocic.\n");
        return 0;

    }

    // utworzenie macierzy jednostkowej

    for (i = 0; i < n; ++i)
    {
        for (j = 0; j < n; ++j)
        {
            if (i == j) Temp[i][n+j] = 1.0;
            else Temp[i][n+j] = 0.0;
        }
    }


    for (i = 0; i < n; ++i)
    {
        if (Temp[i][i] == 0)
        {
            j = i+1;
            while ( Temp[j][i] == 0 && j < n) ++j;
            if (j < n && Temp[j][i] != 0)
                for (k = 0; k < n+n; ++k)
                {
                    tempor2 = Temp[j][k];
                    Temp[j][k] = Temp[i][k];
                    Temp[i][k] = tempor2;
                }

        }
        for (k = n+n-1; k >= 0; --k) Temp[i][k] /= Temp[i][i];
        for (j = 0; j < n; ++j)
        {
            wsp_1 = -Temp[j][i];
            for (k = i; k < n+n; ++k)
                if (j != i)
                    Temp[j][k] += wsp_1 * Temp[i][k];
        }
    }

    // przepisanie odwroconej macierzy do wejsciowej tablicy

    for (i = 0; i < n; ++i)
    {
        for (j = 0; j < n; ++j)
        {
            wynik[i][j] = Temp[i][n+j];
        }
    }

    return 1;
}
//______________________________koniec odwracania macierzy__________________________________



long double norma_1(long double X[n][n])    // sumowanie wyrazow w kolumnach
{
    long double max_norma=0;
    long double temp;

    temp = 0.0;

    for (int i=0; i<n; i++)
    {
        for (int j=0; j<n; j++)
        {
            temp+=fabs(X[j][i]);

        }


        if (temp>max_norma) max_norma = temp;

        temp=0;

    }

    return max_norma;  // zwroc max. sume wyrazow
}

long double norma_1_J_F(long double X[n])    // sumowanie wyrazow w kolumnach
{
    long double max_norma=0;
    long double temp;

    temp = 0.0;

    for (int i=0; i<n; i++)
    {

        temp+=fabs(X[i]);


        if (temp>max_norma) max_norma = temp;

        temp=0;

    }

    return max_norma;  // zwroc max. sume wyrazow
}




void wypelnij_F(long double X[n])
{
    F[0]=funkcja_0(X[0],X[1],X[2]);
    F[1]=funkcja_1(X[0],X[1],X[2]);
    F[2]=funkcja_2(X[0],X[1],X[2]);
}

void wypelnij_J(long double X[n])
{
    J[0][0]=Jacobian_x0(X[1],X[2]);
    J[0][1]=Jacobian_y0(X[0],X[2]);
    J[0][2]=Jacobian_z0(X[2],X[0]);
    J[1][0]=1.0;//Jacobian_x1(X[0]);
    J[1][1]=Jacobian_y1(X[1]);
    J[1][2]=Jacobian_z1(X[2]);
    J[2][0]=0.0;//Jacobian_x2(X[0]);
    J[2][1]=Jacobian_y2(X[1]);
    J[2][2]=-1.0;//Jacobian_z2(X[2]);
}

void mnozenie_J_F(void)
{

    for (int i=0; i<3; i++)
    {
        for (int j=0; j<3; j++)
        {
            for (int k=0; k<3; k++)
            {
                J_F[i]+=J[i][k]*F[k];
            }

        }
    }


}

void pomnozMacierz33(long double A[n][n], long double B[n][n],long double Wyn[n][n])
{ //Wyn = A*B   Macierze 3x3 * 3x3
        int i,j,k;
        for(i=0;i<n;i++)
          for(j=0;j<n;j++)
          {
                Wyn[i][j] = 0;
                for(k=0;k<n;k++)
                  Wyn[i][j] += A[i][k] * B[k][j];
          }
}


void pomnozMacierz31(long double A[n][n], long double B[n],long double Wyn[n])
{ //Wyn = A*B   Macierze 3x3 * 3x1
        int i,j,k,m=1;
        for(i=0;i<n;i++)
          for(j=0;j<m;j++)
          {
                Wyn[i] = 0;
                for(k=0;k<n;k++)
                  Wyn[i] += A[i][k] * B[k];
          }
}


void sprawdz_warunek_zbieznosci(long double X[n])
{


    cout<<"Sprawdzenie warunku zbieżności:"<<endl;
    //wypelnij_J();
    cout<<"MACIERZ J:"<<endl;
    for (int i=0;i<3;i++)
    {
        for (int j=0;j<3;j++)
        {
            cout<<J[i][j]<<"    ";
        }
        cout<<endl;
    }

    //odwracanie_macierzy(3);
    cout<<"MACIERZ J po odwróceniu:"<<endl;
    for (int i=0;i<3;i++)
    {
        for (int j=0;j<3;j++)
        {
            cout<<J[i][j]<<"    ";
        }
        cout<<endl;
    }
    cout<<endl;
//    cout<<"norma to: "<<norma_1(3)<<endl;

//    wypelnij_F();
//    mnozenie_J_F();

    cout<<endl;
//    cout<<"norma J*F to: "<<norma_1_J_F(3)<<endl;

}

int warunek(long double X[n],long double Xn[n])
{   // gdy obie tablice maj&#261; wartosci rowne (z dok&#322;. do eps)
    // funkcja zwraca 0, w przeciwnym wypadku 1
        int i;
        int t[3]={0,0,0};
        long double roznica[3],max=0;
        for(i=0;i<n;i++)
        {
                roznica[i] = fabsl(X[i]-Xn[i]);
                if(roznica[i] < epsilon)t[i]=1;
        }
        for(i=0;i<n;i++)if(roznica[i]>max)max = roznica[i];
        printf("\n\n|X(k)-X(k+1)|: ");
        cout<< max << "\n\n";

        for(i=0;i<n;i++)if(t[i]==0)return 1;
        return 0;
}


void odejmij11(long double X1[n],long double X2[n],long double wynik[n])
{
  for(int i=0;i<n;i++)
  {
      wynik[i]=X1[i]-X2[i];
  }

}

void przepisz(long double in[n], long double out[n])
{