zoukankan      html  css  js  c++  java
  • 多项式简单操作

    多项式求导

    原理

    如果有

    [A(x) = sum_{i = 0}^{n} c_i x^i ]

    那么由于求导的加法原则,每一项可以单独考虑,则有

    [A'(x) = sum_{i = 1}^{n} ic_ix^{i - 1} ]

    实现

    void Der(int *f, int *g, int lf) {
    	For (i, 1, lf) g[i - 1] = 1ll * i * f[i] % Mod; g[lf] = 0;
    }
    

    多项式积分

    原理

    如果有

    [A(x) = sum_{i = 0}^{n} c_i x^i ]

    同样由于积分的加法原则,每一项可以单独考虑,则有(常数项可以忽略)

    [int A(x) mathrm{d} x = sum_{i = 1}^{n + 1} frac{c_{i - 1}}{i} x^{i} ]

    实现

    void Int(int *f, int *g, int lf) {
        g[0] = 0;
        For (i, 1, lf + 1)
            g[i] = 1ll * f[i - 1] * inv[i] % Mod;
    }
    

    多项式乘法

    原理

    定义两个多项式 (A(x) = sum_{i = 0}^{n} a_i x^i, B(x) = sum_{i = 0}^m b_i x^i) 。那么乘积为 (C(x) = sum_{i = 0}^{n} sum_{j = 0}^m a_ib_jx^{i + j})

    我们考虑利用快速傅里叶变换解决,我们求出 (omega_n^k) 代入的所有的点值,然后用逆变换就可以求出所有系数了。

    对于模意义下的我们利用 (g^{frac{p - 1}n}) 代替 (omega) (其中 (g) 为原根)就可以了。

    实现

    我们接下来都默认用费马质数来作为模数进行 ( ext{NTT})

    int rev[Maxn], W[Maxn], len;
    void NTT(int *P, int opt) {
    	Rep (i, len) if (i < rev[i]) swap(P[i], P[rev[i]]);
    	for (int i = 2, p; p = i >> 1, i <= len; i <<= 1) {
    		W[0] = 1; W[1] = fpm(~opt ? g : invg, (Mod - 1) / i);
    		For (k, 2, p - 1) W[k] = 1ll * W[k - 1] * W[1] % Mod;
    		for (int j = 0; j < len; j += i) Rep (k, p) {
    			int u = P[j + k], v = 1ll * P[j + k + p] * W[k] % Mod;
    			P[j + k] = (u + v) % Mod; P[j + k + p] = (u - v + Mod) % Mod;
    		}
    	}
    	if (!~opt) {
    		int invn = fpm(len, Mod - 2);
    		Rep (i, len) P[i] = 1ll * P[i] * invn % Mod;
    	}
    }
    
    void Prepare(int lc) {
    	int cnt = -1; for (len = 1; len <= lc; len <<= 1) ++ cnt;
    	Rep (i, len) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << cnt);
    }
    
    int A[Maxn], B[Maxn], C[Maxn];
    void Mult(int *a, int *b, int *c, int la, int lb) {
    	int lc = la + lb; Prepare(lc);
    	Rep (i, len) A[i] = i <= la ? a[i] : 0; NTT(A, 1);
    	Rep (i, len) B[i] = i <= lb ? b[i] : 0; NTT(B, 1);
    	Rep (i, len) C[i] = 1ll * A[i] * B[i] % Mod; NTT(C, -1);
    	For (i, 0, lc) c[i] = C[i];
    }
    

    多项式牛顿迭代

    原理

    给出函数 (f) ,找到一个多项式 (A),使得 (f(A(x)) equiv 0 pmod {x^n})

    (pmod x) 时,可以求出 (A(x)) 的常数项。

    已经求出 (f(A(x)) equiv 0 mod {x^n}) ,现在求 (f(B(x)) equiv 0 mod {x^{2n}})(B(x))(A(x))(n) 项系数相同。

    (mod {x^{2n}}) 下,在 (A(x)) 处展开 (f(B(x)))

    [f(B(x)) = sum_{i ge 0} frac{f^{(i)}(A(x))}{i!}(B(x) - A(x))^i ]

    不难发现在 (i ge 2) 的时候,最低项已经超过了 (x^{2n}) ,所以只有 (i = 0, 1) 要考虑,也就是

    [egin{aligned} f(B(x)) &equiv f(A(x)) + f'(A(x))(B(x) - A(x)) &pmod {x^{2n}}\ B(x) &equiv A(x) - frac{f(A(x))}{f'(A(x))} &pmod {x^{2n}} end{aligned} ]

    大多数情况下,多项式复合是一个比较难求解的一个东西,但对于一些特殊的 (f) 可以方便求出。

    多项式求逆

    原理

    好像这个不好套牛迭,我们现推算了。

    (F(x))(P(x))(pmod {x^n}) 下的逆元,那么表示出来就是

    [F(x)P(x) equiv 1 pmod {x^n} ]

    (n = 1) 求逆元即可,假设我们求出在 (mod x^{frac n 2}) 意义下的逆元 (G(x)) ,那么有

    [G(x)P(x) equiv 1 pmod {x^frac n 2}\ ]

    我们有

    [egin{aligned} F(x) - G(x) &equiv 0 pmod {x^frac n 2}\ F^2(x) - 2F(x)G(x) + G^2(x) &equiv 0 pmod {x^n}\ F(x) &equiv 2G(x) - P(x)G^2(x) pmod {x^n} end{aligned} ]

    那么递归求解即可,复杂度 (displaystyle mathcal T(n) = mathcal T(frac n2) + mathcal O(n log n) = mathcal O(n log n))

    实现

    由于这个是大部分多项式操作的一个基础操作(加减乘除)之一,所以要尽量加快速度。

    可以先 ( ext{NTT}) 然后做点值乘法,能少一半常数。(注意一开始需要是 (2^k ge len) 形式)

    void Inv(int *f, int *g, int lf) {
    	if (lf == 1) return void(g[0] = fpm(f[0], Mod - 2));
    	Inv(f, g, lf >> 1); Prepare(lf << 1);
    	Rep (i, len) A[i] = i < lf ? f[i] : 0; NTT(A, 1); 
    	Rep (i, len) B[i] = i < lf ? g[i] : 0; NTT(B, 1);
    	Rep (i, len) C[i] = 1ll * A[i] * B[i] % Mod * B[i] % Mod; NTT(C, -1);
    	Rep (i, lf) g[i] = (g[i] * 2ll + Mod - C[i]) % Mod;
    }
    

    多项式开根

    原理

    (F^2(x) equiv P(x) pmod {x^n}) ,此时的 (f(A(x)) = A^2 - P) 我们要求它的零点。

    那么我们套用牛迭的式子就得到了

    [egin{aligned} B(x) &= A(x) - frac{A^2(x) - P(x)}{2A(x)}\ &= frac{A(x)}2 + frac{P(x)}{2A(x)} end{aligned} ]

    复杂度 (displaystyle mathcal T(n) = mathcal T(frac n2) + mathcal O(n log n) = mathcal O(n log n)) 。 。。主定理是真的强大

    注意常数项非 (1) 的时候,通常都会保证存在模意义下的二次剩余,然后用 ( ext{BSGS}) 求出它是 (g^k) ,然后开根就是 (g^{k/2})

    实现

    void Sqrt(int *f, int *g, int lf) {
    	if (lf == 1) return void(g[0] = getsqrt(f[0]));
    	Sqrt(f, g, lf >> 1); static int tmp[Maxn], tinv[Maxn];
    	Rep (i, lf) tmp[i] = 2 * g[i] % Mod; Inv(tmp, tinv, lf);
    	Mult(f, tinv, tmp, lf, lf);
    	const int inv2 = fpm(2, Mod - 2);
    	Rep (i, lf) g[i] = (1ll * g[i] * inv2 % Mod + tmp[i]) % Mod;
    }
    

    多项式ln

    原理

    (G(x) = ln F(x)) ,我们有

    [G(x) = int frac{F'(x)}{F(x)} mathrm{d}x ]

    那么我们求导然后再求逆即可,复杂度 (mathcal O(n log n))

    常数项需要是 (1) ,其他形如 (ln x) 的无理数在模意义下无法表示。

    实现

    void Ln(int *f, int *g, int lf) {
    	static int der[Maxn], tmp[Maxn];
    	Der(f, der, lf); Inv(f, tmp, lf);
    	Mult(der, tmp, tmp, lf, lf); Int(tmp, g, lf);
    }
    

    多项式exp

    原理

    可以套用牛顿迭代和多项式 (ln) 解决。

    分治的做法虽然多了个 (mathcal O(log)) 但好写许多(而且通常由于常数优势会比较快,除非你会论文哥那个牛迭优化),原理如下

    [egin{aligned} B(x) &= exp(A(x))\ B'(x) &= A'(x) exp(A(x))\ B'(x) &= A'(x) B(x)\ B(x) &= int A'(x)B(x) mathrm{d}x end{aligned} ]

    我们考虑把 (A(x)) 先求导成 (A'(x)) ,令 (mid = lfloor frac{l + r}2 floor) 然后考虑分治求出 ([l, mid])(B(x)) ,然后考虑对于 ((mid, r]) 的贡献,也就是 (A_{0 sim r - l})(B_{l sim mid}) 卷然后贡献上去。

    注意 (A(0) = 0) 才在模意义下有含义,此时 (B(0) = 1) ,以及积分需要平移一位并乘上 (frac{1}{n}) 。复杂度 (mathcal O(n log^2 n))

    实现

    void Exp(int *a, int *b, int l, int r) {
    	if (l == r) return void(b[l] = (l ? 1ll * b[l] * inv[l] % Mod : 1));
    	int mid = (l + r) >> 1;
    	Exp(a, b, l, mid); static int tmp[Maxn];
    	Mult(a, b + l, tmp, r - l, mid - l);
    	For (i, mid + 1, r) b[i] = (b[i] + tmp[i - l - 1]) % Mod;
    	Exp(a, b, mid + 1, r);
    }
    
    void Exp(int *f, int *g, int lf) {
    	static int der[Maxn];
    	Der(f, der, lf); der[lf] = 0;
        memset(g, 0, sizeof(int) * (lf + 1));
    	Exp(der, g, 0, lf);
    }
    

    多项式快速幂

    原理

    (G(x) = F^k(x)) ,我们取对数有 (ln G(x) = k ln F(x)) 也就是 (G(x) = exp(k ln F(x)))

    但是 (F(0)) 可能不为 (1) ,那么我们把 (F(x)) 除掉 (F(0)) 即可。那么还可能为 (0) 我们平移即可,也就是除掉 (x^a)

    然后 (exp) 结束记得乘回来它的 (k) 次幂即可。

    代码

    不想判 (0) 了。。

    void Pow(int *f, int *g, int lf, int power) {
    	static int tf[Maxn], tmp[Maxn], invf; assert(f[0]); 
    
    	invf = fpm(f[0], Mod - 2);
    	Rep (i, lf) tf[i] = 1ll * f[i] * invf % Mod;
    
    	Ln(tf, tmp, lf); Rep (i, lf) tmp[i] = 1ll * tmp[i] * power % Mod; Exp(tmp, g, lf);
    
    	invf = fpm(fpm(invf, Mod - 2), power);
    	Rep (i, lf) g[i] = 1ll * g[i] * invf % Mod;
    }
    

    多项式复合逆

    原理

    (G(F(x)) = x)(F(x)) 的一项。此时 (G(x)) 可能是和 (F) 同阶的一个多项式,套用牛顿迭代项数会特别多,根本无法求解。

    但要求单项即 ([x^n]F(x)) 有一个可以快速计算的式子

    [[x^n] F(x) = frac 1n [x^{n - 1}] (frac{x}{G(x)})^n ]

    至于证明可以参看这篇博客

    实现

    int Lagrange(int *f, int lf, int x) {
    	static int tmp[Maxn], res[Maxn];
    	Rep (i, lf) res[i] = f[i + 1]; Inv(res, tmp, lf); Power(tmp, res, n);
    	return 1ll * res[x - 1] * fpm(x, Mod - 2) % Mod;
    }
    
  • 相关阅读:
    Nginx负载均衡+代理+ssl+压力测试
    Nginx配置文件详解
    HDU ACM 1690 Bus System (SPFA)
    HDU ACM 1224 Free DIY Tour (SPFA)
    HDU ACM 1869 六度分离(Floyd)
    HDU ACM 2066 一个人的旅行
    HDU ACM 3790 最短路径问题
    HDU ACM 1879 继续畅通工程
    HDU ACM 1856 More is better(并查集)
    HDU ACM 1325 / POJ 1308 Is It A Tree?
  • 原文地址:https://www.cnblogs.com/zjp-shadow/p/10963718.html
Copyright © 2011-2022 走看看