FFT(n, array)
    if (n == 1)
        return array[0]
    FOR (i = 0 ... n - 1)
        b[i & 1].push_back(array[i])
    f0 = FFT(n / 2, b[0])
    f1 = FFT(n / 2, b[1])
    FOR (j = 0 ... n - 1)
        f[j] = f0[j % (n/2)] + exp(2 * Pi * i * j / n) * f1[j % (n/2)] 
    return f


MUL(n, a, m, b)
    // ^ is not xor here
    k : 2^k >= n + m
    // append zeros to a and b to do they long enough
    f_a = FFT(2^k, a)
    f_b = FFT(2^k, b)
    FOR (i = 0 ... 2^k-1)
        f_c[i] = f_a[i] * f_b[i]
    c = FFT(2^k, f_c)
    reverse(c + 1, c + 2^k)
    FOR (i = 0 ... 2^k-1)
        c[i] = round(c[i] / 2^k)
    return c