BlueStein 算法

BlueStein 算法是处理循环卷积的利器，需要对单位根要有比较清晰的认识。

1. 主要思想

用于循环卷积，即 $c_k = \sum_i \sum_j [i + j\bmod n = k] a_ib_j$。你肯定要说了，这个有什么用？直接计算长度为 $2n$ 的正常卷积就行了？别急，我们待会看道例题就知道了。

2. 算法流程

仍然考虑计算 $a’(k) = \sum_{i = 0} ^ {n - 1} \omega_n ^ {ik} a(i)$，这本质和 FFT 是相同的。注意一般有模数时需要满足循环卷积长度 $n$ 整除 $p - 1$，这样才能使用单位根。

没有比较优秀的性质，怎么办呢？BlueStein 算法给出了一种化积为和的方案。

观察以下两个式子：

$$\begin{aligned} (x + y) ^ 2 - x ^ 2 - y ^ 2 =& 2xy\\ \binom{x + y}2 - \binom x2 - \binom y2 =& xy \end{aligned}$$

正确性展开即可得到。现在我们发现我们已经巧妙的把积放在了和之间了！

一般来说，第二个式子更常用，因为第一个式子有 $2xy$，可能需要用到二次剩余，但是一个数可能不存在二次剩余。

于是我们可以把上面的式子变形一下：

$$\begin{aligned} a'(k) =& \sum_{i = 0} ^ {n - 1} \omega_n ^ {ik} a(i)\\ =& \omega_n ^ {-\binom k2} \sum_{i = 0} ^ {n - 1} \omega_n ^ {\binom{i + k}2}\omega_n ^ {\binom i2}a(i) \end{aligned}$$

容易发现这其实是 $f(x) = \sum_{i = 0} ^ {n - 1} \omega_n ^ {\binom i2}$ 和 $g(x) = \sum_{i = 0} ^ {n - 1} \omega_n ^ {-\binom i2}$ 的差卷积，于是 NTT 加速，复杂度 $O(n\log n)$。这个办法再很多场合比较常见，这也是 BlueStein 算法的一个精妙之处。似乎也叫 Chirp Z-Transform？

然后类似于 FFT 的，我们只需要把 $c’(k) = a’(k) * b’(k)$，然后逆变换回去即可。逆变换只需要在单位根的位置乘上 -1 即可，最后每一项除以 $n$ 即可。

然后考虑证明这个就是循环卷积。同样，由于 BlueStein 算法与 FFT 的相似性，我们同样可以得到 FFT 其实也是循环卷积。

考虑 $a(x), b(y)$ 对于 $c(k)$ 的贡献为：

$$\sum_{i = 0} ^ {n - 1} a(x) \omega_n ^ {ix} b(y) \omega_n ^ {iy} \omega_n ^ {-ik}$$

通过表示 $a(x), b(y)$ 变换之后在哪一位，然后又逆变换回了 $k$ 这一位，就是上面的式子。

考虑推导：

$$\begin{aligned} &\sum_{i = 0} ^ {n - 1} a(x) \omega_n ^ {ix} b(y) \omega_n ^ {iy} \omega_n ^ {-ik}\\ =& a(x) b(y) \sum_{i = 0} ^ {n - 1} \omega_n ^ {ix} \omega_n ^ {iy} \omega_n ^ {-ik}\\ =& a(x) b(y) \sum_{i = 0} ^ {n - 1} \omega_n ^ {(x + y - k)i} \end{aligned}$$

后面一坨式子可以通过 得到即为 $n\times [n|(x + y - k)]$。（关系好像不太对，怎么感觉应该是单位根反演从这个得到的啊

于是我们就得到当 $(x + y)\bmod n = k$ 时，贡献为 $n\times a(x)b(y)$，否则为 0。最后除以 $n$ 即可得到 $c(k)$ 了。

至此，我们就可以在 $O(n\log n)$ 的之间内求出 $a$ 和 $b$ 的循环卷积了。

3. 例题

T1：BZOJ1919 [CTSC2010]性能优化

题目传送门 Luogu

题目传送门 Hydro

题意：做 $A\times B ^ k$ 长度为 $n$ 的循环卷积，对 $n + 1$ 取模。保证 $n + 1$ 是质数。$n\leq 5\times 10 ^ 5$，6s。

提示：为确保你不被卡常，请使用比较快速的任意模数多项式乘法。（但好像可以不用任意模数多项式乘法

这时 BlueStien 的优势体现出来了：我做循环卷积的话，我可以对变换后的点值直接做 $a’(i) \times b’(i) ^ k$，而普通的变换不行，只有倍增快速幂。证明的话，可以考虑 BlueStein 算法逆变换后就是答案，而他和下一个相乘时又要变换回来，这个过程可以省略。这样可以把所有点值一起乘起来，再做逆变换，而普通的变换不行。

注意需要使用 MTT，但是可以使用大模数 NTT 代替，因为本题值域只有 $(5\times 10 ^ 5) ^ 3$，没实现过不知道常数如何（好像比较小）。

void BlueStein(poly &a, int inv)
{
	bool flag = inv == 1;
	if (inv == -1) inv = 1;
	else inv = Mod - 2;
	int n = a.size();
	poly b(2 * n - 1);
	for (int i = 0; i < 2 * n - 1; ++ i)
		b[i] = qpow(W, (i - 1LL) * i / 2 % (Mod - 1) * (Mod - 1 - inv) % (Mod - 1));
	for (int i = 0; i < n; ++ i)
		a[i] = (LL) a[i] * qpow(W, (i - 1LL) * i / 2 % (Mod - 1) * inv % (Mod - 1)) % Mod;
	std::reverse(a.begin(), a.end());
	a = a * b, a = std::vector<int>(a.data() + n - 1, a.data() + 2 * n - 1);
	for (int i = 0; i < n; ++ i)
		a[i] = (LL) a[i] * qpow(W, (i - 1LL) * i / 2 % (Mod - 1) * inv % (Mod - 1)) % Mod;
	if (flag) return;
	inv = qpow(n);
	for (int &x : a) x = (LL) x * inv % Mod;
}

int main()
{
	int n, C;
	std::cin >> n >> C, W = findrt(Mod = n + 1);
	poly a(n), b(n);
	for (int &x : a) scanf("%d", &x);
	for (int &x : b) scanf("%d", &x);
	BlueStein(a, 1), BlueStein(b, 1);
	for (int i = 0; i < n; ++ i) a[i] = (LL) a[i] * qpow(b[i], C) % Mod;
	BlueStein(a, -1);
	for (int &x : a) printf("%d\n", x);
	puts("");
	return 0;
}