形式幂级数 [学习笔记]

形式幂级数

沉迷多项式，无法自拔...

不具体写了看笔记本，这里稍微记一下。

多项式的各种运算
伯努利数
拉格朗日反演

任意模数卷积

我的三模数ntt跑得好慢，然后拆系数fft跑的好快

设(M = lceil P ceil)，将整数表示成(kcdot M+b)的形式

对(x)和(y)进行卷积，分别表示称(acdot M + b, ccdot M + d)

[(acdot M + b)cdot (ccdot M + d) = ac M^2 + (ad+bc)M + bd ]

对(a,b,c,d)进行dft和idft即可

每个数大小在(10^{14})级别，可以使用复数下fft，共进行7次运算

通常M取(32768=2^{15})

~~根据猜测，系数表示不溢出double点值表示就不会溢出double。这玩意应该只能承受一次点值乘法~~

	void mul_any(int *x, int *y, int lim) {
		for(int i=0; i<lim; i++) {
			a[i].x = x[i] >> 15; b[i].x = x[i] & 32767;
			c[i].x = y[i] >> 15; d[i].x = y[i] & 32767;
		}
		dft(a, 1); dft(b, 1); dft(c, 1); dft(d, 1);
		for(int i=0; i<n; i++) {
			cd _a = a[i], _b = b[i], _c = c[i], _d = d[i];
			a[i] = _a * _c;
			b[i] = _a * _d + _b * _c;
			c[i] = _b * _d;
		}
		dft(a, -1); dft(b, -1); dft(c, -1);
		for(int i=0; i<lim; i++) {
			ll _a = (ll) floor(a[i].x + 0.5) %mo, _b = (ll) floor(b[i].x + 0.5) %mo, _c = (ll) floor(c[i].x + 0.5) %mo;
			printf("%lld ", ((_a << 30) %mo + (_b << 15) %mo + _c) %mo);
		}
	}

以下复杂度均为(T(n) = T(n/2) + O(nlogn) =O(nlogn))

模板在最下方

多项式求逆元

求

[A(x) * B(x) equiv 1 pmod {x^n} ]

注意到这时候(A(x)*B(x))的(1...n-1)次项系数为0

用倍增的思想，已知(mod x^{lceil frac{l}{2} ceil})的逆元(B_0(x))求(mod x^l)下的逆元(B(x))

(l=1)时，(b_0 = a_0^{-1})，可以发现多项式有逆的充要条件是常数项有逆

两式相减，然后平方，同乘(A(x))，得到

[B(x) equiv B_0(x) * (2 - A(x) * B_0(x)) pmod {x^l} ]

处理(mod x^l)时，(l)就是当前的次数界，次数(ge l)的都整除没了。

可以理解为只关心前l项

多项式开根

求

[B^2(x) equiv A(x) pmod {x^n} ]

同样倍增的思想

已知(mod x^{lceil frac{l}{2} ceil})的平方根(B_0(x))求(mod x^l)下的平方根(B(x))

(l=1)时，(b_0 equiv sqrt{a_0} pmod x)，可能需要二次剩余

移项化简后得到

[B(x) equiv 2^{-1}cdot (B_0(x) + A(x) * B_0^{-1}(x)) pmod{x^n} ]

同时还需要求逆...

牛顿迭代法

给出(G(x))，求(F(x))，

[G(F(x)) equiv 0 pmod{x^n} ]

倍增的思想。将(G(F(X)))在(F_0(x))处泰勒展开得到

[F(x) equiv F_0(x) - frac{G(F_0(x))}{G'(F_0(x))} pmod{x^l} ]

也可以用这个式子求逆元和开根，最后的结果式子一样。

多项式求ln和exp就是和对应的麦克劳林级数复合，所以要求常数项为0

多项式求ln

给出(F(x) = 1 + sum_{i ge 1}f_ix^i)

求一下导

[ln F(x) = int frac{F'(x)}{F(x)}dx ]

多项式求exp

给出(A(x) = sum_{i ge 1}a_ix^i)

[e^{A(x)} - B(x) equiv 0 pmod {x^n} ]

取对数后使用牛顿迭代法

[B(x) = B_0(x)(1 - ln B_0(x) + A(x)) ]

多项式k次幂

当(A(x))的常数项为1

[A(x)^k = exp(k ln A(x)) ]

否则提取最低次项(ax^d)

[A(x)^k = (ax^d)^k (frac{A(x)}{ax^d})^k ]

伯努利数

wiki

用来解决等幂求和问题

[egin{align} sum_{i=0}^{n-1} i^m = frac{1}{m+1}sum_{i=0}^m inom{m+1}{i}B_i^- n^{m+1-i} \ S_m(n) = sum_{i=1}^{n} i^m = frac{1}{m+1}sum_{i=0}^m inom{m+1}{i}B_i^+ n^{m+1-i} \ end{align} ]

复杂度与幂次有关

除了(B_1)，其他奇数项都是0。(B_1^+ = frac{1}{2},B_1^-=-frac{1}{2})

(0^0=1)

递推关系

令(n=1, m eq 0),

[sum_{i=0}^m inom{m+1}{i}B_i^- = 0, B_0=1, B_1^-=-frac{1}{2} ]

指数型生成函数

对于(B^-),

[sum_{i=0}^infty B_ifrac{x^i}{i!} =frac{x}{e^x - 1} = (sum_{i=0}^infty frac{x^i}{(i+1)!})^{-1} ]

使用多项式求逆元即可预处理伯努利数. 求(mod x^{n+1})意义下逆元

预处理任意模数伯努利数的模板

拉格朗日反演

复合逆(反函数)：

没有常数项的(f(x), g(x))，(f(g(x))=x)，那么互为复合逆，(f_1g_1=1,g(f(x))=x)

用拉格朗日反演可以(O(nlogn))求复合逆某一项的系数

[[x^n]g(x) = frac{1}{n}[omega^{n-1}](frac{omega}{f(omega)})^n ]

可以配合多项式k次幂使用。

生成函数中出现x之后可以用啦。

模板

namespace ntt {
	int g = 3, rev[N];
	void dft(int *a, int n, int flag) {
		int k = 0; while((1<<k) < n) k++;
		for(int i=0; i<n; i++) {
			rev[i] = (rev[i>>1]>>1) | ((i&1)<<(k-1));
			if(i < rev[i]) swap(a[i], a[rev[i]]);
		}
		for(int l=2; l<=n; l<<=1) {
			int m = l>>1, wn = Pow(g, flag == 1 ? (P-1)/l : P-1-(P-1)/l);
			for(int *p = a; p != a+n; p += l)
				for(int k=0, w=1; k<m; k++, w = (ll)w*wn %P) {
					int t = (ll) w * p[k+m] %P, r = p[k];
					p[k+m] = (r - t + P) %P;
					p[k] = (r + t) %P;
				}
		}
		if(flag == -1) {
			ll inv = Pow(n, P-2);
			for(int i=0; i<n; i++) a[i] = a[i] * inv %P;
		}
	}
	
	void inverse(int *a, int *b, int l) {
		static int t[N];
		if(l == 1) {b[0] = Pow(a[0], P-2); return;}
		inverse(a, b, l>>1);
		int n = l<<1;
		for(int i=0; i<l; i++) t[i] = a[i], t[i+l] = 0;
		dft(t, n, 1); dft(b, n, 1);
		for(int i=0; i<n; i++) b[i] = (ll) b[i] * (2 - (ll) t[i] * b[i] %P + P) %P;
		dft(b, n, -1); for(int i=l; i<n; i++) b[i] = 0;
	}

	void sqrt(int *a, int *b, int l) {
		static int t[N], ib[N];
		if(l == 1) {b[0] = 1; return;}
		sqrt(a, b, l>>1);
		int n = l<<1;
		for(int i=0; i<l; i++) t[i] = a[i], t[i+l] = 0, ib[i] = ib[i+l] = 0;
		inverse(b, ib, l);
		dft(t, n, 1); dft(b, n, 1); dft(ib, n, 1);
		for(int i=0; i<n; i++) b[i] = (ll) inv2 * (b[i] + (ll) t[i] * ib[i] %P) %P;
		dft(b, n, -1); for(int i=l; i<n; i++) b[i] = 0;
	}

	void ln(int *a, int *b, int l) {
		static int da[N], ia[N];
		int n = l<<1;
		for(int i=0; i<n; i++) da[i] = ia[i] = 0;
		for(int i=0; i<l-1; i++) da[i] = (ll) (i+1) * a[i+1] %P;
		inverse(a, ia, l);
		dft(da, n, 1); dft(ia, n, 1);
		for(int i=0; i<n; i++) b[i] = (ll) da[i] * ia[i] %P;
		dft(b, n, -1);
		for(int i=l-1; i>0; i--) b[i] = (ll) inv[i] * b[i-1] %P; b[0] = 0;
		for(int i=l; i<n; i++) b[i] = 0;
	}
	
	void exp(int *a, int *b, int l) {
		static int t[N];
		if(l == 1) {b[0] = 1; return;}
		exp(a, b, l>>1);
		int n = l<<1;
		for(int i=0; i<n; i++) t[i] = 0;
		ln(b, t, l);
		for(int i=0; i<l; i++) t[i] = (a[i] - t[i] + P) %P; t[0] = (t[0] + 1) %P;
		dft(b, n, 1); dft(t, n, 1);
		for(int i=0; i<n; i++) b[i] = (ll) b[i] * t[i] %P;
		dft(b, n, -1); for(int i=l; i<n; i++) b[i] = 0;
	}

	void power(int *a, int k, int l) {
		static int t[N];
		int n = l<<1;
		for(int i=0; i<n; i++) t[i] = 0;
		ln(a, t, l);
		for(int i=0; i<l; i++) t[i] = (ll) k * t[i] %P;
		for(int i=0; i<n; i++) a[i] = 0;
		exp(t, a, l);
	}
}