equation

#include <iostream>
#include <cmath>
#include <vector>
#include <stack>
typedef long long i64;
using namespace std;
bool increment(vector<i64>&pair, vector<vector<i64>>&v);
vector<i64> perebor(vector<vector<i64>>& v, vector<i64>& p, vector<i64>& idempotents);
vector<i64>decomposition(i64 n);
vector<i64>find_x(vector<i64>& v, i64 n);
i64 calc_value(vector<i64>& v, i64 x, i64 n);
bool diophant(i64 a, i64 b, i64 c, i64 &x, i64 &y);
vector < i64 > idempotents_from_vector(vector <i64> & v);
vector<i64> find_idempotents(vector<i64> &idempotents, vector<i64> &v_z);
i64 calc_number(vector<i64>&idempotents, vector<i64>& v_z, vector<i64>& v_ost);

int main(){
    i64 a,m,n; vector <i64> v;
    cout<<"Input degree of polynom(m) and base of Zn: ";
    cin>>m>>n;
    cout<<"Input a[0], a[1], ..., a[m] from polynom a[0]*x^m+a[1]*x^(m-1)+...+a[m] = 0: ";
    for (i64 i=0; i<=m; ++i){
        cin>>a;
        v.push_back(a);
    }
    cout<<endl<<"SOLUTION: "<<endl;
    vector<i64> p = decomposition(n);
    vector<i64> mas = idempotents_from_vector(p);
    vector<i64> idempotents = find_idempotents(mas, p);
    vector<vector<i64>> x;
    for (auto &i:p){
        x.push_back(find_x(v,i));
        cout<<"for p = "<<i<<": ";
        for(auto &j : x.back())cout<<j<<" ";
        cout<<endl;
    }
    cout << "ALLRIGHT" << endl;
    vector<i64> answer= perebor(x, p, idempotents);
    cout<<"ANSWER: ";
    for(auto &i : answer) cout<<i<<" ";
    cout<<endl<<"CHECK THIS: "<<endl;
    for(auto &x : answer)
    cout<<"for x="<<x<<": polynom(x) = "<< calc_value(v, x, n)<<endl;
}


bool increment(vector<i64>&pair, vector<vector<i64>>&v){
    for (i64 i = 0; i< v.size(); ++i){
        pair[i]++;
        if (pair[i] == v[i].size()) pair[i]=0; else  return true;
    }
    return false;
}

vector<i64> perebor(vector<vector<i64>>& v, vector<i64>& p, vector<i64>& idempotents){
    vector<i64>pair, w, answer;
    for(i64 i=0; i<v.size(); ++i){ pair.push_back(0); w.push_back(0);}
    bool flag;
    do {
        for (i64 i=0; i<v.size(); ++i) w[i]=v[i][pair[i]];
        answer.push_back(calc_number(idempotents, p, w));
        flag = increment(pair, v);
    } while (flag);
    return answer;
}

vector<i64>decomposition(i64 n){
    vector<i64>v;
    for (i64 i = 2; i <= n; ++i)
        if (n % i == 0)
        {
            v.push_back(i);
            n /= i;
            while (n % i == 0)
            {
                n /= i;
                v.back() *= i;
            }
        }
  return v;
}

vector<i64>find_x(vector<i64>& v, i64 n){
    vector<i64>x;
    for (i64 i=0; i<n; ++i){
        i64 y = calc_value(v, i, n);
        if (y==0)x.push_back(i);
    }
    return x;
}

i64 calc_value(vector<i64>& v, i64 x, i64 n){
    i64 b=0;
    for(auto &a : v){
        b = x*b + a;
        while (b>=n) b-=n;
        while (b<0) b+=n;
    }
    return b;
}

bool diophant(i64 a, i64 b, i64 c, i64 &x, i64 &y){
//  Решение диофантового уравнения a*x + b*y = c методом Евклида
    stack < i64 >q;
    i64 r;
    while (a % b != 0)
    {
        r = a % b;
        q.push(a / b);
        a = b;
        b = r;
    }
    if (b != 1)
    {
        return false;
    }
    else
    {
        i64 temp;
        x = 1;
        y = q.top();
        q.pop();
        while (!q.empty())
        {
            temp = x;
            x = -y;
            y = -temp - q.top() * y;
            q.pop();
        }
        x *= c;
        y *= -c;
        return true;
    }
}

vector < i64 > idempotents_from_vector(vector <i64> & v){
//  Нахождение идемпотентов в кольце вычетов по модулю n, где n - произведение всех элементов вектора v, состоящего из взаимно простых чисел
    i64 n = 1, m;
//  Вычисление n
    for (auto &i : v) n *= i;
//  Двоичный перебор всех представлений n = a * b , где НОД(a,b) = 1;
    int c = v.size();
    i64 limit = pow(2, c - 1), a, b, x, y;
    vector < i64 > answer; // Будет хранить найденные идемпотенты
//  Достаточно перебрать всего 2^(c-1) диофантовых уравнений, откуда мы вытащим один идемпотент
//  И для найденного идемпотента х посчитаем 1-х парный идемпотент
    for (i64 i = 1; i < limit; ++i)
    {
        m = i;
        a = b = 1;
        for (i64 j = 0; j < c; ++j)
        {
            if (m % 2 == 1)
                a *= v[j];
            else
                b *= v[j];
            m /= 2;
        }
        // Решаем полученное диофантово уравнение a*x + b*y = c
        diophant(a, b, 1, x, y);
        // a*x - один идемпотент
        x *= a;
        // Приводим х к нормальному виду: 0 <= x < n
        while (x < 0) x += n;
        while (x >= n) x -= n;
        // Находим парный к нему идемпотент
        y = n + 1 - x;
        // Кидаем оба идемпотента в вектор
        answer.push_back(x);
        answer.push_back(y);
    }
//  Ну и возвращаем вектор
    return answer;
}

vector<i64> find_idempotents(vector<i64> &idempotents, vector<i64> &v_z){
//  Вектор idempotents хранит все идемпотенты кольца вычетов с основанием n
//  Вектор v_z хранит все основания колец вычетов , введенные пользователем
//  Вектор v_ost хранит все остатки от сравнения по модулю с соответсвующим основанием кольца вычетов
    i64 n = 1 ;
//  Так как мы не передаём n , то посчитаем его как произведение всех элементов вектора v_z
    vector<i64> answer;
    for (auto &i : v_z) { n *= i; answer.push_back(0);}
//  Затем сопоставим каждому j-ому кольцу вычетов v_z[j] i-ый идемпотент кольца вычетов по модулю n
    for (auto &a : idempotents){
        int t = -1, k = 0;
        bool flag = true;
        for(auto &b : v_z){
            // Ищем единственную возможную единицу в остатке при делении a[i] на b[j] .
            // Все остальные остатки при делении a[i] на b[k], где k != j, должны быть нулями
            if (a % b == 1)
                if (t != -1){
                    flag = false;
                    break;
                } else
                    t = k;
            ++k;
        }
        // Если мы нашли соответсвующий идемпотент для v_z[t], то сразу умножаем его на v_ost[t] и прибавляем к искомому числу
        if ( flag )
            answer[t] = a;
    }
// Возвращаем искомое число
    return answer;
}


i64 calc_number(vector<i64>&idempotents, vector<i64>& v_z, vector<i64>& v_ost){
    i64 n=1, num =0;
    for (auto &i : v_z) n *= i;
    for (int i=0; i<v_z.size(); ++i){
        num += idempotents[i]*v_ost[i];
        if (num>=n) num%=n;
    }
    return num;
}