【Luogu3803】模板：多项式乘法

Link:
Luogu https://www.luogu.com.cn/problem/P3803

---

Solution

日常复习一波 FFT

比较重要的就是，DFT 是代入单位根求点值
然后奇偶拆分下放

[omega_n^k=cosfrac{2kpi}{n}+isinfrac{2kpi}{n} ]

[omega_{n}^{k+frac{n}{2}}=-omega_n^k ]

然后 IDFT 的时候

[c_k=na_k ]

就是要除以 (n) 注意一下。

首先一点，模板题（多项式乘法）的确做的是卷积
但是为什么不是 f[i] * g[N-i] 呢？
因为 N = n + m
（

[x^ncdot x^m=x^{n+m} ]

好吧。

提前丢人：今天又把 n/2 + i/2 负优化成了 n + i >> 1
complex step 写错了几次，一次是写成 /lim
一次是没写 cos 和 sin
另外数组一开始开小了（）
开了 1<<20 然后两个点 RE

总而言之下面是递归FFT
非递归等一等哦。

#include<cstdio>
#include<cstdlib>
#include<iostream>
#include<algorithm>
#include<cctype>
#include<cmath>

using namespace std;

const int MAXN = 1 << 22;

struct complex
{
    double real, imag;
    complex(double aug1 = 0, double aug2 = 0)
    {
        real = aug1;
        imag = aug2;
    }
    void operator = (const complex &b)
    {
        real = b.real;
        imag = b.imag;
    }
    complex operator + (const complex &b)
    {
        return complex(real + b.real, imag + b.imag);
    }
    complex operator - (const complex &b)
    {
        return complex(real - b.real, imag - b.imag);
    }
    complex operator * (const complex &b)
    {
        return complex(real * b.real - imag * b.imag, real * b.imag + imag * b.real);
    }
    complex operator * (const double &x)
    {
        return complex(real * x, imag * x);
    }
    //下面不能写 const& 不然自乘可能会炸
    void operator *= (complex b)
    {
        *this = *this * b;
    }
} temp[MAXN], cf[MAXN], cg[MAXN], cur;

const double tpi = acos(-1.0);
int n, m, lim = 1;
double modpi;

// rev == +1 :  DFT
// rev == -1 : IDFT
void FFT(complex *f, int n, int rev)
{
    if (n <= 1)
    {
        return;
    }

    int half_n = n >> 1;

    for (int i = 0; i < n; ++i)
    {
        temp[i] = f[i];
    }

    for (int i = 0; i < n; ++i)
    {
        if (i & 1) 
        {
            f[(n >> 1) + (i >> 1)] = temp[i];
        }
        else
        {
            f[i >> 1] = temp[i];
        }
    }

    complex *g = f;
    complex *h = f + half_n;

    FFT(g, half_n, rev);
    FFT(h, half_n, rev);

    cur = complex(1, 0);
    complex step(cos(tpi / half_n), sin(tpi / half_n * rev));

    for (int i = 0; i < half_n; ++i)
    {
        temp[i] = g[i] + h[i] * cur;
        temp[i + half_n] = g[i] - h[i] * cur;
        cur *= step;
    }

    for (int i = 0; i < n; ++i)
    {
        f[i] = temp[i];
    }
}

int main()
{

    ios::sync_with_stdio(false);

    cin >> n >> m;

    while (lim <= n + m)
    {
        lim <<= 1;
    }

    modpi = 2.0 * tpi / lim;

    for (int i = 0; i <= n; ++i)
    {
        cin >> cf[i].real;
    }   

    for (int i = 0; i <= m; ++i)
    {
        cin >> cg[i].real;
    }

    FFT(cf, lim, +1);
    FFT(cg, lim, +1);

    for (int i = 0; i < lim; ++i)
    {
        cf[i] *= cg[i];
    }

    FFT(cf, lim, -1);

    n += m;
    for (int i = 0; i <= n; ++i)
    {
        cout << (int)(cf[i].real/lim+0.05) << ' ';
    }

    // system("pause");

    return 0;
}

---

继续。
讲到递归FFT，就不得不提一下蝴蝶变换
那么究竟为什么可以二进制翻转呢？我想要提供一个比较直观的理解角度

然后我也不知道为什么突然就想通了

比如说第一次就是按照二进制最后一位是 0 还是 1 分成左右两边
然后左边的，新的下标最高位都会是 0；右边的都是 1
就完成了原下标最低位 → 新下标最高位的翻转。

比如这样：
000 001 010 011 | 100 101 110 111
变成
000 010 100 110 | 001 011 101 111 这是原下标
0xx 0xx 0xx 0xx | 1xx 1xx 1xx 1xx 这是新下标
接下来只考虑未确定的部分
考虑其中一边就可以了，另外一边类比
比如说右边，
001 011 | 101 111 原下标
然后：
$00 $01 | $10 $11 暂时下标，等于原下标除以 2 的商，也就是原下标第一位 ~ 倒数第二位
所以这个暂时下标的二进制末位是原下标二进制倒数第二位
$00 $10 | $01 $11 按照暂时下标奇偶分开之后的下标
所以，现在这半边确定的新下标是
00x 00x | 01x 01x

至于
Rev[i] = (Rev[i>>1]>>1) | ((i&1)<<(Log - 1));
这个我就不解释了。。懂的都懂
注意一下里面利用了之前翻转的结果 Rev[i>>1] 应该马上就明白了

非递归：

void FFT(Complex* A, const bool& Type = 0)
{
	static Complex t;
	for (R int i = 0; i < Limit; ++i) if(Rev[i] > i) swap(A[i], A[Rev[i]]);
	for (R int Mid = 1; Mid < Limit; Mid <<= 1)
	{
		for (R int Len = Mid << 1, Pos = 0; Pos < Limit; Pos += Len)
		{
			for (R int Delta = 0; Delta < Mid; ++Delta)
			{
				t = Wn[Limit/Len*Delta][Type] * A[Pos + Mid + Delta];
				A[Pos + Mid + Delta] = A[Pos + Delta] - t;
				A[Pos + Delta] += t;
			}
		}
	}
}