拉插、FFT 和 NTT

吃老本

拉傅数

拉格朗日插值法

一个次数为 $n-1$ 的多项式可以被 $n$ 个点确定，具体地，

$$F(k)=\sum_{i=1}^{n}y_i\prod_{j\neq i}\frac{k-x_j}{x_i-x_j} $$

代码

int n;
ll x[N], y[N], k;
ll ans;

ll qpow(ll a, ll b) {
	ll ans = 1;
	for (; b; b >>= 1, a = a * a % mod) 
		if (b & 1) ans = ans * a % mod;
	return ans;
} 

int main() {
	scanf ("%d%lld", &n, &k);
	for (int i = 1; i <= n; i++)
		scanf ("%lld%lld", &x[i], &y[i]);
	for (int i = 1; i <= n; i++)
	{
		ll tmp = y[i];
		for (int j = 1; j <= n; j++)
			if(i != j)
				tmp = tmp * ((k - x[j] + mod) % mod) % mod * qpow((x[i] - x[j] + mod) % mod, mod - 2) % mod;
		ans = (ans + tmp) % mod;
	}
	printf ("%lld\n", ans);
	return 0;
}

快速傅里叶变换 FFT

求卷积。

从系数表示转成点值表示（DFT）来乘，然后再转回系数表示（IDFT）。

对于一个 $n-1$ 次的多项式，其点值表示法：

$$F(x)=\{(x_0,F(x_0)),(x_1,F(x_1)),\cdots,(x_{n-1},F(x_{n-1}))\} $$

注意次数，点值表示的原理就是，一个次数为 $n-1$ 的多项式可以被 $n$ 个点确定。对于其中的 $x_i$ ，我们从一个单位圆中等分 $n$ 个点，以它们表示的复数 $\omega_n^i$ 代入，并称之为单位根。DFT 过后数组里存储的值 f[i] 就是 $F(\omega_n^i)$ 。

单位根有很显然的性质：

  1. $\omega_{n}^{k}=\omega_{2n}^{2k}$ ，可以用在分治；
  2. $\omega_{n}^{k + \frac{n}{2}} = -\omega_{n}^{k}$ ，可以用于讨论大于半圆的单位根情况；
  3. 点值表示的值乘上单位根的逆元能回到系数表示乘 $n$ ，结合等比数列求和手推可得之。这也是为什么选单位根。

分治优化 DFT

设多项式 $F(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}$ ，然后按照奇偶性分为两份：

$$\begin{aligned}F(x)&=(a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})+(a_1x+a_3x^3+\cdots+a_{n-1}x^{n-1})\\ &=(a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})+x(a_1+a_3x^2+\cdots+a_{n-1}x^{n-2})\end{aligned}$$

设

$$F_1(x)=a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{\frac{n}{2}-1}\\ F_2(x)=a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{\frac{n}{2}-1}$$

那么，

$$F(x)=F_1(x^2)+x\cdot F_2(x^2) $$

将 $x=\omega_n^k\quad(k<\frac{n}{2})$ 代入得

$$\begin{aligned}F(\omega_n^k)&=F_1(\omega_n^{2k})+\omega_n^k\cdot F_2(\omega_n^{2k})\\ &=F_1(\omega_{\frac{n}{2}}^{k})+\omega_n^k\cdot F_2(\omega_{\frac{n}{2}}^{k}) \end{aligned}$$

将 $x=\omega_n^{k+\frac{n}{2}}$ 代入得

$$\begin{aligned}F(\omega_n^{k+\frac{n}{2}})&=F_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}\cdot F_2(\omega_n^{2k+n})\\ &=F_1(\omega_n^{2k}\cdot\omega_n^n)+\omega_n^{k+\frac{n}{2}}\cdot F_2(\omega_n^{2k}\cdot\omega_n^n)\\ &=F_1(\omega_{\frac{n}{2}}^{k})-\omega_n^{k}\cdot F_2(\omega_{\frac{n}{2}}^{k})\\ \end{aligned}$$

然后我们就可以通过 $\mathcal{O}(n\log n)$ 的时间复杂度递归了。

IDFT 就如同单位根性质 3 所述。

蝴蝶变换优化

递归 FFT 常数大，所以尝试把它变成递推形式。

按奇偶性划分可以发现一个规律：

十进制 二进制
$0~1~2~3~4~5~6~7$ $000~001~010~011~100~101~110~111$
$0~2~4~6,1~3~5~7$ $000~010~100~110,001~011~101~111$
$0~4,2~6,1~5,3~7$ $000~100,010~110,001~101,011~111$
$0,4,2,6,1,5,3,7$ $000,100,010,110,001,101,011,111$

每一个数最终的位置，就是把它的二进制翻转过来。

非常好考验位运算能力，令 $n$ 为 $2$ 的幂，设 tr[i] 表示 $i$ 的最终位置（ $i<n$ ），则有 tr[i] = (tr[i >> 1] >> 1) | (i & 1? n >> 1: 0)。其意义（或曰原理）就是，对于 $i$ ，考虑把末位插到前 $|i|-1$ 位的值的前面（其中 $|i|=\lfloor\log_2i\rfloor$ ）。

所以每个数的位置就能够 $\mathcal{O}(n)$ 预处理出来。

代码

const int N = 1500010;

const double PI = acos(-1.0);

struct Complex {
	double x, y;
	Complex (double ix, double iy) {x = ix, y = iy;}
	Complex () {x = y = 0;}
	Complex operator + (Complex &b) {
		return Complex(x + b.x, y + b.y);
	}
	Complex operator - (Complex &b) {
		return Complex(x - b.x, y - b.y);
	}
	Complex operator * (Complex &b) {
		return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
	}
}f[N << 1], g[N << 1];

int n, m;
int tr[N << 1];

void FFT (Complex *f, bool isDFT) {
	for (int i = 1; i <= n; i++)
		if (i < tr[i]) swap (f[i], f[tr[i]]);
	
	for (int p = 2; p <= n; p <<= 1) {  // p 表示这次循环结束后，合并的块长
		int len = p >> 1;  // 上次循环合并的块长
		Complex omega(cos(2 * PI / p), sin(2 * PI / p));
		if (!isDFT) omega.y *= -1;
		for (int k = 0; k < n; k += p) {
			Complex tmp(1, 0);
			for (int i = k; i < k + len; i ++) {
				Complex y = tmp * f[i + len];
				f[i + len] = f[i] - y;
				f[i] = f[i] + y;
				tmp = tmp * omega;
			}
		}
	}
    if(!isDFT) for (int i = 0; i <= n; i++) f[i].x /= n;
}

int main() {
	scanf ("%d%d", &n, &m);
	for (int i = 0; i <= n; i++) scanf ("%lf", &f[i].x);
	for (int i = 0; i <= m; i++) scanf ("%lf", &g[i].x);
		
	for (m += n, n = 1; n <= m; n <<= 1);
	for (int i = 1; i <= n; i++)
		tr[i] = (tr[i >> 1] >> 1) | (i & 1? n >> 1: 0);
	
	FFT(f, 1), FFT(g, 1);
	
	for (int i = 0; i <= n; i++) f[i] = f[i] * g[i];
		
	FFT(f, 0);
	for (int i = 0; i <= m; i++) printf ("%d ", (int)(f[i].x + 0.49));
	return 0;
}

快速数论变换 NTT

由于 FFT 会有很大的精度误差，所以用数论代替，就是 NTT 了，NTT 也比 FFT 要快很多。

我超，原（根）

若整数 $g$ 是奇素数 $p$ 的一个原根，则满足 $g,g^2,g^3,\cdots,g^{p-1}$ 在模 $p$ 意义下互不相同。

在模 $998244353$ 意义下最小原根 $g=3$ !

NTT 的基本流程

若 $p$ 满足 $2^k\cdot r+1$ （比如 $998244353=2^{23}\times199+1$ ），把 $g^{\frac{p-1}{2^k}}$ 作为 $\omega_n^1$ ，发现原本单位根的性质它都能满足，那么就这么跑 FFT 可以了。

但 NTT 也有局限性，只能处理 $n\leq 2^k$ 的情况，遇到模数 $p=10^9+7$ 时，NTT 就做不了。不过有任意模数 NTT。

代码

const int N = 1500010;

const ll mod = 998244353ll, G = 3, invG = 332748118ll;

ll f[N << 1], g[N << 1];

int n, m;
int tr[N << 1];

ll qpow(ll a, ll b) {
	ll ans = 1;
	for (; b; b >>= 1, a = a * a % mod)
		ans = ans * (b & 1? a: 1) % mod;
	return ans;
}

void NTT (ll *f, bool isDFT) {
	for (int i = 1; i <= n; i++)
		if (i < tr[i]) swap (f[i], f[tr[i]]);
	
	for (int p = 2; p <= n; p <<= 1) {
		int len = p >> 1;
		ll omega = qpow(isDFT? G: invG, (mod - 1) / p);
		for (int k = 0; k < n; k += p) {
			ll tmp = 1ll;
			for (int i = k; i < k + len; i ++) {
				ll y = tmp * f[i + len] % mod;
				f[i + len] = (f[i] - y + mod) % mod;
				f[i] = (f[i] + y) % mod;
				tmp = tmp * omega % mod;
			}
		}
	}
    if(!isDFT) {
		ll invn = qpow(n, mod - 2);
		for (int i = 0; i <= n; i++) f[i] = f[i] * invn % mod;
	}
}

int main() {
	scanf ("%d%d", &n, &m);
	for (int i = 0; i <= n; i++) scanf ("%lld", &f[i]);
	for (int i = 0; i <= m; i++) scanf ("%lld", &g[i]);
		
	for (m += n, n = 1; n <= m; n <<= 1);
	for (int i = 1; i <= n; i++)
		tr[i] = (tr[i >> 1] >> 1) | (i & 1? n >> 1: 0);
	
	NTT(f, 1), NTT(g, 1);
	
	for (int i = 0; i <= n; i++) f[i] = f[i] * g[i];
		
	NTT(f, 0);
	for (int i = 0; i <= m; i++) printf ("%lld ", f[i]);
	return 0;
}

更快的 NTT

	#define clr(f, x) memset(f, 0, sizeof(int) * (x))
	#define cpy(f, g, x) memcpy(f, g, sizeof(int) * (x))
	int add (int x, int y) { return x + y > mod? x + y - mod: x + y; }
	int dec (int x, int y) { return x - y < 0? x - y + mod: x - y; }
	int mul (int x, int y) { return 1ll * x * y % mod; }
	int qpow (int a, ll b) {
		int ans = 1;
		for (; b; b >>= 1, a = mul(a, a))
			if (b & 1) ans = mul(ans, a);
		return ans;
	}

	struct poly {
		int tr[pN];

		void init(int n, int siz) {
			for (int i = 0; i < n; ++i)
				tr[i] = (tr[i >> 1] >> 1) | ((i & 1) << (siz - 1));
		}

		void NTT(int *f, int n, int op = 1) {
			for (int i = 0; i < n; i++)
				if (i < tr[i]) swap(f[i], f[tr[i]]);
			for (int p = 1; p < n; p <<= 1) {
				int omega = qpow(op? G: invG, (mod - 1) / (p << 1));
				for (int j = 0; j <	n; j += p << 1) 
					for (int w = 1, k = 0; k < p; k++, w = mul(w, omega)) {
						int x = f[j | k], y = mul(w, f[j | p | k]);
						f[j | k] = add(x, y), f[j | p | k] = dec(x, y);
					}
			}
			if (!op) {
				int invn = qpow(n, mod - 2);
				for (int i = 0; i < n; i++) f[i] = mul(f[i], invn);
			}
		}
		void times (int *f, int *g, int n, int m, int T) {
			int lim = 1, lsiz = 0; for (; lim < n + m; lim <<= 1, lsiz++);
			init(lim, lsiz);
			static int tmp[pN]; clr(f + n, lim - n), cpy (tmp, g, m), clr(tmp + m, lim - m);
			NTT(f, lim), NTT(tmp, lim);
			for (int i = 0; i < lim; i++) f[i] = mul(f[i], tmp[i]);
			NTT(f, lim, 0);
			clr (f + T, lim - 1), clr (tmp, lim);
		}
	} p;
EOF

评论

暂无评论

发表评论

可以用@mike来提到mike这个用户，mike会被高亮显示。如果你真的想打“@”这个字符，请用“@@”。

博客信息