Метод Гаусса по модулю

Ребят вот написал код решение систем методом Гаусса по модулю:

#include <iostream>
using namespace std;
int gcdex(int a, int b, int &x, int &y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    int x1, y1;
    int d1 = gcdex(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return d1;
}
int ReverseElement(int a, int N) {
    int x, y, d;
    d = gcdex(a, N, x, y);
    if (d != 1) {
        return 0;
    }
    else
    {
        if (x < 0)
            x += N;
        return x;
    }
}
// Вывод системы уравнений
void sysout(double **a, double *y, int n)
{
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            cout << a[i][j] << "*x" << j;
            if (j < n - 1)
                cout << " + ";
        }
        cout << " = " << y[i] << endl;
    }
    return;
}
double * gauss(double **a, double *y, int n, int p)
{
    double *x, max;
    int k, index;
    const double eps = 0.00001;  // точность
    x = new double[n];
    k = 0;
    while (k < n)
    {
        // Поиск строки с максимальным a[i][k]
        max = abs(a[k][k]);
        index = k;
        for (int i = k + 1; i < n; i++)
        {
            if (abs(a[i][k]) > max)
            {
                max = abs(a[i][k]);
                index = i;
            }
        }
        // Перестановка строк
        if (max < eps)
        {
            // нет ненулевых диагональных элементов
            cout << "Решение получить невозможно из-за нулевого столбца ";
            cout << index << " матрицы A" << endl;
            return 0;
        }
        for (int j = 0; j < n; j++)
        {
            double temp = a[k][j];
            a[k][j] = a[index][j];
            a[index][j] = temp;
        }
        double temp = y[k];
        y[k] = y[index];
        y[index] = temp;
        // Нормализация уравнений
        for (int i = k; i < n; i++)
        {
            double temp = ReverseElement(a[i][k],p);
            if (abs(temp) < eps) continue; // для нулевого коэффициента пропустить
            for (int j = 0; j < n; j++)
                a[i][j] = a[i][j] * temp;
            y[i] = y[i] * temp;
            if (i == k)  continue; // уравнение не вычитать само из себя
            for (int j = 0; j < n; j++)
                a[i][j] = a[i][j] - a[k][j];
            y[i] = y[i] - y[k];
        }
        k++;
    }
    // обратная подстановка
    for (k = n - 1; k >= 0; k--)
    {
        x[k] = int(y[k])%p;
        if (x[k] < 0)
            x[k] += p;
        for (int i = 0; i < k; i++)
            y[i] = y[i] - a[i][k] * x[k];
    }
    return x;
}
int main()
{
    int p = 46;
    double **a, *y, *x;
    int n;
    system("chcp 1251");
    system("cls");
    cout << "Введите количество уравнений: ";
    cin >> n;
    a = new double*[n];
    y = new double[n];
    for (int i = 0; i < n; i++)
    {
        a[i] = new double[n];
        for (int j = 0; j < n; j++)
        {
            cout << "a[" << i << "][" << j << "]= ";
            cin >> a[i][j];
        }
    }
    for (int i = 0; i < n; i++)
    {
        cout << "y[" << i << "]= ";
        cin >> y[i];
    }
//  sysout(a, y, n);
    x = gauss(a, y, n,p);
    for (int i = 0; i < n; i++)
        cout << "x[" << i << "]=" << x[i] << endl;
    cin.get(); cin.get();
    return 0;
}

Он почти такой же как и обычный, только тут надо не делить на наибольший элемент в строке, а умножать на обратный элемент по алгоритму Эвклида. Отлично работает на матрице:

101 1
030 8
500 12 mod46

Проблема в том что, возьмем например матрицу:

101 1
110 2
021 7 mod46

Если будем решать по методу Гаусса при первом занулении столбца все проходит хорошо. Получается так:

101 1
01-1 1
021 7

Потом на втором круге, мы делаем перестановку 3 строки со 2, так как коэффициент больше:

101 1
021 7
01-1 1

И здесь мы берем обратный элемент 2 по модулю 46. Такого элемента не существует и я не знаю как дальше быть. Что делать дальше?