FFT

namespace
{
    const double PI = std::atan(1.0) * 4;
    int reverse_bits(int n, int log_n)
    {
        int rn = 0;
        for (int i = 0; i < log_n; i++)
        {
            rn <<= 1;
            rn += (n >> i) & 1;
        }
        return rn;
    }
    long long gcd_extended(long long a, long long b, long long &x, long long &y)
    {
        if (a == 0)
        {
            x = 0;
            y = 1;
            return b;
        }
        long long x1, y1;
        long long g = gcd_extended(b % a, a, x1, y1);
        x = y1 - (b / a) * x1;
        y = x1;
        return g;
    }

    long long inverse_elem(long long a, long long m)
    {
        long long x, y;
        gcd_extended(a, m, x, y);
        return (x % m + m) % m;
    }

    const long long mod = 7340033;
    const long long base_prim_root = 5;
    const long long base_prim_root_pow = (1 << 20);
}

void fast_fourier_transform::fft_integer(vector<long long> &data, bool invert)
{
    vector<long long> tmp = data;
    int logn = 0;
    for (size_t i = 1; i < data.size(); i *= 2)
        logn++;
    for (size_t i = 0; i < data.size(); i++)
        tmp[i] = data[reverse_bits(i, logn)];// переставляем числа чтобы избавиться от рекурсии
//  for (size_t i = 0; i < data.size(); i++)
//      tmp[i] %= mod; //??
    data = tmp;

    for (size_t len = 2; len <= data.size(); len *= 2)//перебираем длины подотрезков
    {
        long long cur_prim_root = invert ? inverse_elem(base_prim_root, mod): base_prim_root;
        for (long long i = len; i < base_prim_root_pow; i *= 2)
            cur_prim_root = (cur_prim_root * cur_prim_root) % mod;

        for (size_t i = 0; i < data.size(); i += len) //цикл по подотрезкам
        {
            long long w = 1;
            for (size_t j = 0; j < len / 2; j++) // по подотрезку [i, i + len)
            {
                long long c = data[i + j];
                long long d = (data[i + len / 2 + j] * w) % mod;
                data[i + j] = (c + d) % mod;
                data[i + len / 2 + j] = (c - d >= 0 ? c - d : c - d + mod);
                w = (w * cur_prim_root) % mod;
            }
        }
    }
    if (invert)
    {
        long long n_inv = inverse_elem(data.size(), mod);
        for (size_t i = 0; i < data.size(); i++)
            data[i] = (data[i] * n_inv) % mod;
    }
}