void normalize_submatrix(double **data, size_t k, size_t n)
{
    for (size_t i = k; i < n; ++i)
    {
        double temp = data[i][k];
        if (!(fabs(temp) < EPS))
        {
            for (size_t j = 0; j < n + 1; ++j)
                data[i][j] /= temp;
            if (i != k)
                for (size_t j = 0; j < n + 1; ++j)
                    data[i][j] -= data[k][j];
        }
    }
}

double find_max_in_submatrix(double **data, size_t n, size_t k, size_t *index_i, size_t *index_j)
{
    double max = fabs(data[k][k]);
    *index_i = k;
    *index_j = k;
    for (size_t i = k; i < n; ++i)
        for (size_t j = k; j < n; ++j)
            if (fabs(data[i][j]) - max > EPS)
            {
                max = fabs(data[i][j]);
                *index_i = i;
                *index_j = j;
            }

    return max;
}

void swap_pointers(double **a, double **b)
{
    double *temp = *a;
    *a = *b;
    *b = temp;
}

void swap_columns(double **data, size_t n, size_t k, size_t index, size_t *lead_columns)
{
    for (size_t i = 0; i < n; ++i)
    {
        double temp = data[i][k];
        data[i][k] = data[i][index];
        data[i][index] = temp;
    }
    size_t temp = lead_columns[k];
    lead_columns[k] = lead_columns[index];
    lead_columns[index] = temp;
}

int gauss(double **data, double *answers, size_t n)
{
    double max;
    size_t *lead_columns = malloc(n * sizeof(size_t));
    if (!lead_columns)
        return ERROR_MEMORY;
    for (size_t i = 0; i < n; ++i)
        lead_columns[i] = i;
    for (size_t k = 0; k < n; ++k)
    {
        size_t index_i, index_j;
        max = find_max_in_submatrix(data, n, k, &index_i, &index_j);
        if (max < EPS)
        {
            free(lead_columns);
            return ERROR_NOT_SOLVABLE;
        }

        swap_columns(data, n, k, index_j, lead_columns);
        swap_pointers(data + k, data + index_i);
        normalize_submatrix(data, k, n);
    }

    for (size_t k = n; k > 0; --k)
    {
        answers[lead_columns[k - 1]] = data[k - 1][n];
        for (size_t i = 0; i < k - 1; ++i)
            data[i][n] -= data[i][k - 1] * answers[lead_columns[k - 1]];
    }
    free(lead_columns);

    return OK;
}