文章目录
- 多项式的逆元
- 理论推导
- 模板
- 例题:[luogu P4238]【模板】多项式乘法逆
- 题目
- code
- 多项式的除法/取模
- 理论推导
- 多项式牛顿迭代法
- 模板
- 例题:[luoguP4512]【模板】多项式除法
- 题目
- code
- 多项式对数ln
- 理论推导
- 模板
- 例题
- 题目
- code
- 多项式开根sqrt
- 理论推导
- 模板
- 例题
- 题目
- code
- 多项式指数exp
- 理论推导
- 模板
- 例题
- 题目
- code
- 多项式快速幂
- 理论推导
- 例题
- 题目
- code
由于本蒟蒻比较差,每一份代码都是交ACACAC了才来的,所以不用怀疑它的真实行,当然代码和模板都不是luogu上面加强版,所以常数项恰好都是111,温馨提醒:不要乱用
多项式的逆元
理论推导
知识储备知识储备知识储备
(a+b)%c=(a%c+b%c)%c(a+b)\%c=(a\%c+b\%c)\%c(a+b)%c=(a%c+b%c)%c
(a−b)%c=(a%c−b%c+c)%c(a-b)\%c=(a\%c-b\%c+c)\%c(a−b)%c=(a%c−b%c+c)%c
(a∗b)%c=(a%c)∗(b%c)%c(a*b)\%c=(a\%c)*(b\%c)\%c(a∗b)%c=(a%c)∗(b%c)%c
也就是说如果你的算法只使用了加法、减法和乘法,你可以在所有的中间步骤取模,这和只对答案取模是等效的
同余的一些运算
a≡b(modm),c≡d(modm)=>a∗c≡b∗d(modm)a\equiv b ( \mathrm{mod\:} m ),c\equiv d ( \mathrm{mod\:} m )=>a*c \equiv b*d ( \mathrm{mod\:} m )a≡b(modm),c≡d(modm)=>a∗c≡b∗d(modm)
a≡b(modm),c≡d(modm)=>a±c≡b±d(modm)a\equiv b ( \mathrm{mod\:} m ),c\equiv d ( \mathrm{mod\:} m )=>a±c \equiv b±d ( \mathrm{mod\:} m )a≡b(modm),c≡d(modm)=>a±c≡b±d(modm)
a≡b(modm)=>ka≡kb(modkm)a\equiv b ( \mathrm{mod\:} m )=>ka \equiv kb ( \mathrm{mod\:} km )a≡b(modm)=>ka≡kb(modkm)
ac≡bc(modm)=>a≡b(modm(c,m))ac\equiv bc ( \mathrm{mod\:} m )=>a \equiv b ( \mathrm{mod\:} \frac{m}{(c,m)} )ac≡bc(modm)=>a≡b(mod(c,m)m)
接下来进入正题:
设给定的多项式为A(x)A(x)A(x),求A−1(x)A^{-1}(x)A−1(x)满足:👇
A(x)∗A−1(x)≡1(modxn)A(x)*A^{-1}(x)\equiv1(mod\ x^n)A(x)∗A−1(x)≡1(mod xn)
modxnmod\ x^nmod xn的意思即为舍去多项式里次数≥n\ge n≥n的项,只保留次数∈[0,n−1]∈[0,n-1]∈[0,n−1]的项
首先如果只有常数项,我们可以直接快速幂算出逆元(利用费马小定理)
假设已知A(x)A(x)A(x)在关于modx⌈n2⌉mod\ x^{\lceil\frac{n}{2}\rceil}mod x⌈2n⌉意义下的逆元B0(x)B_0(x)B0(x),满足
A(x)∗B0(x)≡1(modx⌈n2⌉)A(x)*B_0(x)\equiv1(mod\ x^{\lceil\frac{n}{2}\rceil})A(x)∗B0(x)≡1(mod x⌈2n⌉)
记要求的答案为B(x)B(x)B(x),有
A(x)∗B(x)≡1(modxn)A(x)*B(x)\equiv1(mod\ x^n)A(x)∗B(x)≡1(mod xn)
两式相减,则有
A(x)∗(B(x)−B0(x))≡0(modx⌈n2⌉)A(x)*(B(x)-B_0(x))\equiv0(mod\ x^{\lceil\frac{n}{2}\rceil})A(x)∗(B(x)−B0(x))≡0(mod x⌈2n⌉)
因为A(x)A(x)A(x)肯定不为000,所以上述式子要成立的话只能寄希望于(B(x)−B0(x))(B(x)-B_0(x))(B(x)−B0(x)),所以化简得
B(x)−B0(x)≡0(modx⌈n2⌉)B(x)-B_0(x)\equiv0(mod\ x^{\lceil\frac{n}{2}\rceil})B(x)−B0(x)≡0(mod x⌈2n⌉)
我们怎么把modx⌈n2⌉mod\ x^{\lceil\frac{n}{2}\rceil}mod x⌈2n⌉升级到modxnmod\ x^nmod xn,很简单 将式子平方一下
可以得到
(B(x)−B0(x))2≡0(modxn)(B(x)-B_0(x))^2\equiv0(mod\ x^n)(B(x)−B0(x))2≡0(mod xn)
展开这个式子:
B2(x)+B02(x)−2∗B(x)∗B0(x)≡0(modxn)B^2(x)+B_0^2(x)-2*B(x)*B_0(x)\equiv0(mod\ x^n)B2(x)+B02(x)−2∗B(x)∗B0(x)≡0(mod xn)
我们来证明一下为什么可以modxnmod\ x^nmod xn,思考多项式乘法的过程,ci,ai,bic_i,a_i,b_ici,ai,bi都表示各自多项式中xix^ixi的系数
ci=∑j=0iaj∗bi−j……①c_i=\sum_{j=0}^{i}a_j*b_{i-j}……①ci=j=0∑iaj∗bi−j……①
B(x)−B0(x)≡0(modx⌈n2⌉)……②B(x)-B_0(x)\equiv0(mod\ x^{\lceil\frac{n}{2}\rceil})……②B(x)−B0(x)≡0(mod x⌈2n⌉)……②
首先这两个式子的正确性不用证明了吧!我们把B(x)−B0(x)B(x)-B_0(x)B(x)−B0(x)看成一个整体C(x)C(x)C(x),虽然不能保证B,B0B,B_0B,B0在modx⌈n2⌉mod\ x^{\lceil\frac{n}{2}\rceil}mod x⌈2n⌉意义下的系数都为000,但是因为②②②我们可知C(x)C(x)C(x)在modx⌈n2⌉mod\ x^{\lceil\frac{n}{2}\rceil}mod x⌈2n⌉意义下的系数都为000,那么平方就转换成了C(x)∗C(x)C(x)*C(x)C(x)∗C(x),套用①①①
ansi=∑j=0icj∗ci−jans_i=\sum_{j=0}^ic_j*c_{i-j}ansi=j=0∑icj∗ci−j
发现在<n<n<n的项至少是由其中一个C(x)C(x)C(x)中的<⌈n2⌉<\lceil\frac{n}{2}\rceil<⌈2n⌉项乘以其他项的结果,但是当=n=n=n时,是会遇到c⌈n2⌉∗c⌈n2⌉c_{\lceil\frac{n}{2}\rceil}*c_{\lceil\frac{n}{2}\rceil}c⌈2n⌉∗c⌈2n⌉,结合上面的公式B(x)−B0(x)≡0(modx⌈n2⌉)B(x)-B_0(x)\equiv0(mod\ x^{\lceil\frac{n}{2}\rceil})B(x)−B0(x)≡0(mod x⌈2n⌉),x⌈n2⌉x^{\lceil\frac{n}{2}\rceil}x⌈2n⌉是被舍掉了的,无法保证一定为000,
所以平方后的式子[0,n−1][0,n-1][0,n−1]次项都是000,但是nnn次项不确定,就成为一条分水岭,故把取模的界点设在nnn可以砍掉,而不设在n−1n-1n−1或者n+1n+1n+1等位置
大家最关心的时间复杂度,因为不停的/2/2/2用递归,就是logloglog级别,所以O(nlogn)O(nlogn)O(nlogn)左右皆大欢喜
模板
void solve ( int deg, int *b ) {if ( deg == 1 ) {b[0] = qkpow ( a[0], mod - 2 );return;}solve ( ( deg + 1 ) >> 1, b );len = 1;l = 0;while ( len < deg * 2 - 1 ) {len <<= 1;l ++;}inv = qkpow ( len, mod - 2 );for ( int i = 0;i < len;i ++ )r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );for ( int i = 0;i < deg;i ++ )c[i] = a[i];for ( int i = deg;i < n;i ++ )c[i] = 0;NTT ( c, 1 );NTT ( b, 1 );for ( int i = 0;i < len;i ++ )b[i] = ( 2 - b[i] * c[i] % mod + mod ) % mod * b[i] % mod; NTT ( b, -1 );for ( int i = deg;i < n;i ++ )b[i] = 0;
}
例题:[luogu P4238]【模板】多项式乘法逆
题目
给定一个多项式 F(x)F(x)F(x) ,请求出一个多项式 G(x)G(x)G(x), 满足 F(x)∗G(x)≡1(modxn)F(x) * G(x) \equiv 1 ( \mathrm{mod\:} x^n )F(x)∗G(x)≡1(modxn)系数对998244353998244353998244353 取模。
输入格式
首先输入一个整数 nnn, 表示输入多项式的项数
接着输入nnn 个整数,第 iii 个整数 aia_iai代表 F(x)F(x)F(x)次数为i−1i-1i−1的项的系数。
输出格式
输出 nnn个数字,第iii个整数 bib_ibi,代表 G(x)G(x)G(x)次数为 i−1i-1i−1 的项的系数。 保证有解。
输入输出样例
输入
5
1 6 3 4 9
输出
1 998244347 33 998244169 1020
说明/提示
1≤n≤105,0≤ai≤1091 \leq n \leq 10^5, 0 \leq a_i \leq 10^91≤n≤105,0≤ai≤109
code
998244353998244353998244353的原根是333,所以就直接用了
#include <cstdio>
#include <iostream>
using namespace std;
#define mod 998244353
#define MAXN 300005
#define int long long
int n, l, pi, ni, len, inv;
int a[MAXN], b[MAXN], c[MAXN], r[MAXN];int qkpow ( int x, int y ) {int ans = 1;while ( y ) {if ( y & 1 )ans = ans * x % mod;x = x * x % mod;y >>= 1;}return ans;
}void NTT ( int *c, int f ) {for ( int i = 0;i < len;i ++ )if ( i < r[i] )swap ( c[i], c[r[i]] );for ( int i = 1;i < len;i <<= 1 ) {int omega = qkpow ( f == 1 ? pi : ni, ( mod - 1 ) / ( i << 1 ) );for ( int j = 0;j < len;j += ( i << 1 ) ) {int w = 1;for ( int k = 0;k < i;k ++, w = w * omega % mod ) {int x = c[k + j], y = w * c[k + j + i] % mod;c[k + j] = ( x + y ) % mod;c[k + j + i] = ( x - y + mod ) % mod;}}}if ( f == -1 )for ( int i = 0;i < len;i ++ )c[i] = c[i] * inv % mod;
}void polyinv ( int deg, int *a, int *b ) {//求多项式a的逆元b if ( deg == 1 ) {b[0] = qkpow ( a[0], mod - 2 );return;}polyinv ( ( deg + 1 ) >> 1, a, b );len = 1;l = 0;while ( len < deg * 2 - 1 ) {len <<= 1;l ++;}inv = qkpow ( len, mod - 2 );for ( int i = 0;i < len;i ++ )r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );for ( int i = 0;i < deg;i ++ )c[i] = a[i];for ( int i = deg;i < n;i ++ )c[i] = 0;NTT ( c, 1 );NTT ( b, 1 );for ( int i = 0;i < len;i ++ )b[i] = ( 2 - b[i] * c[i] % mod + mod ) % mod * b[i] % mod; NTT ( b, -1 );for ( int i = deg;i < n;i ++ )b[i] = 0;
}signed main() {pi = 3;ni = mod / pi + 1;scanf ( "%lld", &n );for ( int i = 0;i < n;i ++ )scanf ( "%lld", &a[i] );polyinv ( n, a, b );for ( int i = 0;i < n;i ++ )printf ( "%lld ", b[i] );return 0;
}
多项式的除法/取模
理论推导
给定一个n−1n-1n−1次的多项式A(x)A(x)A(x)和m−1m-1m−1次的多项式B(x)B(x)B(x),求D(x),R(x)D(x),R(x)D(x),R(x)满足
A(x)=D(x)∗B(x)+R(x)……①A(x)=D(x)*B(x)+R(x)……①A(x)=D(x)∗B(x)+R(x)……①
其中自然D(x)D(x)D(x)的最高次为n−mn-mn−m,R(x)R(x)R(x)最高次<m−1<m-1<m−1
或者A(x)≡R(x)(modB(x))A(x)\equiv R(x)(mod\ B(x))A(x)≡R(x)(mod B(x))
因为有余数R(x)R(x)R(x)难以对其进行处理,所以考虑去掉这个玩意儿带来的不良影响
定义反转操作,即将各项系数反转 dalao让我们这么搞,只能跟着
AR(x)=xn−1A(1x)=∑i=0n−1an−i−1xi……②A^R(x)=x^{n-1}A(\frac{1}{x})=\sum_{i=0}^{n-1}a_{n-i-1}x^i……②AR(x)=xn−1A(x1)=i=0∑n−1an−i−1xi……②
证明很简单,把中间的那一步A(1x)A(\frac{1}{x})A(x1)展开再与外面的xn−1x^{n-1}xn−1相乘,你会发现是一样的,好下一步
把1x\frac{1}{x}x1代入①①①,再同乘各多项式的xn−1x^{n-1}xn−1得到
xn−1A(1x)=xn−1∗D(1x)∗B(1x)+xn−1∗R(1x)x^{n-1}A(\frac{1}{x})=x^{n-1}*D(\frac{1}{x})*B(\frac{1}{x})+x^{n-1}*R(\frac{1}{x})xn−1A(x1)=xn−1∗D(x1)∗B(x1)+xn−1∗R(x1)
我们把xn−1x^{n-1}xn−1拆分成xn−m∗xm−1x^{n-m}*x^{m-1}xn−m∗xm−1,因为它们分别对应了D(x),B(x)D(x),B(x)D(x),B(x)的最高次,再把xn−1x^{n-1}xn−1拆分成xn−m+1∗xm−2x^{n-m+1}*x^{m-2}xn−m+1∗xm−2,因为xm−2x^{m-2}xm−2是R(x)R(x)R(x)能达到的最高次,我们就强买强卖要求R(x)R(x)R(x)达到最高次,大不了系数为0,即
xn−1A(1x)=xn−mD(1x)∗xm−1B(1x)+xn−m+1∗xm−2R(1x)x^{n-1}A(\frac{1}{x})=x^{n-m}D(\frac{1}{x})*x^{m-1}B(\frac{1}{x})+x^{n-m+1}*x^{m-2}R(\frac{1}{x})xn−1A(x1)=xn−mD(x1)∗xm−1B(x1)+xn−m+1∗xm−2R(x1)
由②②②可转换为
AR(x)=DR(x)BR(x)+xn−m+1RR(x)A^R(x)=D^R(x)B^R(x)+x^{n-m+1}R^R(x)AR(x)=DR(x)BR(x)+xn−m+1RR(x)
由于DR(x)D^R(x)DR(x)是n−mn-mn−m次的,所以在modxn−m+1mod\ x^{n-m+1}mod xn−m+1意义下,就可以得到AR(x)≡DR(x)∗BR(x)(modxn−m+1)A^R(x)\equiv D^R(x)*B^R(x)\ (mod\ x^{n-m+1})AR(x)≡DR(x)∗BR(x) (mod xn−m+1),猥琐目的达到
所以我们可以用多项式求逆得到DR(x)D^R(x)DR(x),反转即为D(x)D(x)D(x),然后再带入求R(x)R(x)R(x)即可
所以上面的操作真的是妙极了!!!
多项式牛顿迭代法
有一个关于多项式f(x)f(x)f(x)的方程g(f(x))=0g(f(x))=0g(f(x))=0
假设已经求出f(x)f(x)f(x)的前nnn项f0(x)f_0(x)f0(x),则有
f(x)≡f0(x)(modxn)f(x)\equiv f_0(x)\ (mod\ x^n)f(x)≡f0(x) (mod xn)
g(f0(x))≡0(modxn)g(f_0(x))\equiv 0\ (mod\ x^n)g(f0(x))≡0 (mod xn)
对g(f(x))g(f(x))g(f(x))在f(x)f(x)f(x)上进行泰勒展开
g(f(x))=g(f0(x))+g′(f0(x))1!(f(x)−f0(x))1+g′′(f0(x))2!(f(x)−f0(x))2+……g(f(x))=g(f_0(x))+\frac{g'(f_0(x))}{1!}(f(x)-f_0(x))^1+\frac{g''(f_0(x))}{2!}(f(x)-f_0(x))^2+……g(f(x))=g(f0(x))+1!g′(f0(x))(f(x)−f0(x))1+2!g′′(f0(x))(f(x)−f0(x))2+……
与上面内容结合f(x)−f0(x)f(x)-f_0(x)f(x)−f0(x)的前nnn项系数为000,于是g(f(x))≡g(f0(x))+g′(f0(x))(f(x)−f0(x))≡0(modx2n)g(f(x))\equiv g(f_0(x))+g'(f_0(x))(f(x)-f_0(x))\equiv 0\ (mod\ x^{2n})g(f(x))≡g(f0(x))+g′(f0(x))(f(x)−f0(x))≡0 (mod x2n)
为什么从第三项开始就可以省略了呢??里面有清晰的证明
即f(x)≡f0(x)−g(f0(x))g′(f0(x))(modx2n)f(x)\equiv f_0(x)-\frac{g(f_0(x))}{g'(f_0(x))}\ (mod\ x^{2n})f(x)≡f0(x)−g′(f0(x))g(f0(x)) (mod x2n)
g(f(x))≡lnf(x)−A(x)g(f(x))\equiv lnf(x)-A(x)g(f(x))≡lnf(x)−A(x)
这个函数的零点为eA(x)e^A(x)eA(x),即A(x)A(x)A(x)的expexpexp
简单来说,多项式牛顿迭代法就是用来解关于多项式的方程的
——摘自某dalao
用这种方法也可以推多项式求逆,但是更多的是后面的一些基础操作,先卖个关子
模板
void polydiv ( int n, int m,int *d, int *a, int *b, int *r ) {for ( int i = 0;i < n;i ++ )//进行反转操作 A[i] = a[n - i - 1];for ( int i = 0;i < m;i ++ )B[i] = b[m - i - 1];polyinv ( n - m + 1, B, Binv );//根据推导的公式求出B在mod x^(n-m+1)意义下的逆len = 1;l = 0;while ( len < ( n << 1 ) - m ) {len <<= 1;l ++;}for ( int i = 0;i < len;i ++ )rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );NTT ( A, 1, len );NTT ( Binv, 1, len );for ( int i = 0;i < len;i ++ )//求反转的商D d[i] = A[i] * Binv[i] % mod;NTT ( d, -1, len ); for ( int i = 0, j = n - m;i < j;i ++, j -- )//把商反转成正常 swap ( d[i], d[j] );//接下来搞r memset ( B, 0, sizeof ( B ) );memset ( Binv, 0, sizeof ( Binv ) );memset ( rev, 0, sizeof ( rev ) );len = 1;l = 0;while ( len < n ) {l ++;len <<= 1;}for ( int i = 0;i < len;i ++ )rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) ); for ( int i = 0;i < m;i ++ )B[i] = b[i];for ( int i = 0;i < n - m + 1;i ++ )Binv[i] = d[i];NTT ( B, 1, len );NTT ( Binv, 1, len );for ( int i = 0;i < len;i ++ )r[i] = B[i] * Binv[i] % mod;NTT ( r, -1, len );for ( int i = 0;i < m;i ++ )r[i] = ( a[i] - r[i] + mod ) % mod;
}
例题:[luoguP4512]【模板】多项式除法
题目
给定一个 nnn 次多项式 F(x)F(x)F(x) 和一个 mmm 次多项式 G(x)G(x)G(x) ,请求出多项式 Q(x),R(x)Q(x), R(x)Q(x),R(x),满足以下条件:
Q(x)Q(x)Q(x) 次数为 n−mn−mn−m,R(x)R(x)R(x) 次数小于 mmm
F(x)=Q(x)∗G(x)+R(x)F(x) = Q(x) * G(x) + R(x)F(x)=Q(x)∗G(x)+R(x)
所有的运算在模 998244353998244353998244353 意义下进行。
输入格式
第一行两个整数nnn,mmm,意义如上。
第二行 n+1n+1n+1个整数,从低到高表示 F(x)F(x)F(x) 的各个系数。
第三行 mmm 个整数,从低到高表示 G(x)G(x)G(x) 的各个系数。
输出格式
第一行 n−m+1n-m+1n−m+1 个整数,从低到高表示 Q(x)Q(x)Q(x) 的各个系数。
第二行 mmm 个整数,从低到高表示 R(x)R(x)R(x) 的各个系数。
如果 R(x)R(x)R(x) 不足 m−1m-1m−1 次,多余的项系数补 000。
输入输出样例
输入
5 1
1 9 2 6 0 8
1 7
输出
237340659 335104102 649004347 448191342 855638018
760903695
说明/提示
对于所有数据,1≤m<n≤1051 \le m < n \le 10^51≤m<n≤105 ,给出的系数均属于 [0,998244353)∩Z[0, 998244353) \cap \mathbb{Z}[0,998244353)∩Z
code
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
#define mod 998244353
#define MAXN 300005
#define int long long
int l, pi, ni, len, inv;
int a[MAXN], b[MAXN], c[MAXN], r[MAXN], d[MAXN], rev[MAXN];
int A[MAXN], B[MAXN], Binv[MAXN];int qkpow ( int x, int y ) {int ans = 1;while ( y ) {if ( y & 1 )ans = ans * x % mod;x = x * x % mod;y >>= 1;}return ans;
}void NTT ( int *c, int f, int len ) {for ( int i = 0;i < len;i ++ )if ( i < rev[i] )swap ( c[i], c[rev[i]] );for ( int i = 1;i < len;i <<= 1 ) {int omega = qkpow ( f == 1 ? pi : ni, ( mod - 1 ) / ( i << 1 ) );for ( int j = 0;j < len;j += ( i << 1 ) ) {int w = 1;for ( int k = 0;k < i;k ++, w = w * omega % mod ) {int x = c[k + j], y = w * c[k + j + i] % mod;c[k + j] = ( x + y ) % mod;c[k + j + i] = ( x - y + mod ) % mod;}}}inv = qkpow ( len, mod - 2 ); if ( f == -1 )for ( int i = 0;i < len;i ++ )c[i] = c[i] * inv % mod;
}void polyinv ( int deg, int *a, int *b ) {if ( deg == 1 ) {b[0] = qkpow ( a[0], mod - 2 );return;}polyinv ( ( deg + 1 ) >> 1, a, b );len = 1;l = 0;while ( len < deg * 2 - 1 ) {len <<= 1;l ++;}inv = qkpow ( len, mod - 2 );for ( int i = 0;i < len;i ++ )rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );for ( int i = 0;i < deg;i ++ )c[i] = a[i];for ( int i = deg;i < len;i ++ )c[i] = 0;NTT ( c, 1, len );NTT ( b, 1, len );for ( int i = 0;i < len;i ++ )b[i] = ( 2 - b[i] * c[i] % mod + mod ) % mod * b[i] % mod; NTT ( b, -1, len );for ( int i = deg;i < len;i ++ )b[i] = 0;
}void polydiv ( int n, int m,int *d, int *a, int *b, int *r ) {for ( int i = 0;i < n;i ++ )//进行反转操作 A[i] = a[n - i - 1];for ( int i = 0;i < m;i ++ )B[i] = b[m - i - 1];polyinv ( n - m + 1, B, Binv );//根据推导的公式求出B在mod x^(n-m+1)意义下的逆len = 1;l = 0;while ( len < ( n << 1 ) - m ) {len <<= 1;l ++;}for ( int i = 0;i < len;i ++ )rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );NTT ( A, 1, len );NTT ( Binv, 1, len );for ( int i = 0;i < len;i ++ )//求反转的商D d[i] = A[i] * Binv[i] % mod;NTT ( d, -1, len ); for ( int i = 0, j = n - m;i < j;i ++, j -- )//把商反转成正常 swap ( d[i], d[j] );//接下来搞r memset ( B, 0, sizeof ( B ) );memset ( Binv, 0, sizeof ( Binv ) );memset ( rev, 0, sizeof ( rev ) );len = 1;l = 0;while ( len < n ) {l ++;len <<= 1;}for ( int i = 0;i < len;i ++ )rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) ); for ( int i = 0;i < m;i ++ )B[i] = b[i];for ( int i = 0;i < n - m + 1;i ++ )Binv[i] = d[i];NTT ( B, 1, len );NTT ( Binv, 1, len );for ( int i = 0;i < len;i ++ )r[i] = B[i] * Binv[i] % mod;NTT ( r, -1, len );for ( int i = 0;i < m;i ++ )r[i] = ( a[i] - r[i] + mod ) % mod;
}signed main() {pi = 3;ni = mod / pi + 1;int n, m;scanf ( "%lld %lld", &n, &m );n ++;m ++; for ( int i = 0;i < n;i ++ )scanf ( "%lld", &a[i] );for ( int i = 0;i < m;i ++ )scanf ( "%lld", &b[i] );polydiv ( n, m, d, a, b, r );for ( int i = 0;i < n - m + 1;i ++ )printf ( "%lld ", d[i] );printf ( "\n" );for ( int i = 0;i < m - 1;i ++ )printf ( "%lld ", r[i] );return 0;
}
多项式对数ln
理论推导
对于一个给定的多项式A(x)A(x)A(x),求
lnA(x)lnA(x)lnA(x)
直接运用求导和积分完成,需要保证A(x)A(x)A(x)的常数项为1
对于f(x),xf(x),xf(x),x点的导数:f′(x)=f(x+△)−f(x)△f'(x)=\frac{f(x+△)-f(x)}{△}f′(x)=△f(x+△)−f(x)
lnA(x)=∫ln(A(x))′=∫A(x)′A(x)lnA(x)=∫ln(A(x))'=∫\frac{A(x)'}{A(x)}lnA(x)=∫ln(A(x))′=∫A(x)A(x)′
反正对于小孩子而言,我是难以理解求导和积分的,可能这里会留个坑,等学到后面,再慢慢跟他们玩吧,如果有dalao能教教我(可私信),我一定感激死你
模板
void polyln ( int n, int *a, int *b ) {polyinv ( n, a, b );for ( int i = 1;i < n;i ++ )tmp[i - 1] = a[i] * i % mod;int len = 1, l = 0;while ( len < n * 2 - 1 ) {len <<= 1;l ++;}for ( int i = 0;i < len;i ++ )rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );NTT ( tmp, 1, len );NTT ( b, 1, len );for ( int i = 0;i < len;i ++ )b[i] = tmp[i] * b[i] % mod;NTT ( b, -1, len );for ( int i = n - 1;i; i -- )b[i] = b[i - 1] * qkpow ( i, mod - 2 ) % mod;b[0] = 0; for ( int i = n;i < len;i ++ )b[i] = 0;
}
例题
题目
luogu
给出 n−1n−1n−1 次多项式 A(x)A(x)A(x),求一个 modxn\bmod{\:x^n}modxn下的多项式B(x)B(x)B(x),满足 B(x)≡lnA(x)B(x) \equiv \ln A(x)B(x)≡lnA(x).
在 mod 998244353\text{mod } 998244353mod 998244353 下进行,且 ai∈[0,998244353]∩Za_i\in [0, 998244353] \cap \mathbb{Z}ai∈[0,998244353]∩Z
输入格式
第一行一个整数 nnn
下一行有 nnn 个整数,依次表示多项式的系数 a0,a1,⋯,an−1a_0, a_1, \cdots, a_{n-1}a0,a1,⋯,an−1
输出格式
输出 nnn 个整数,表示答案多项式中的系数 a0,a1,⋯,an−1a_0, a_1, \cdots, a_{n-1}a0,a1,⋯,an−1
输入输出样例
输入
6
1 927384623 878326372 3882 273455637 998233543
输出
0 927384623 817976920 427326948 149643566 610586717
说明/提示
对于 100%100\%100% 的数据,n≤105n \le 10^5n≤105
code
#include <cstdio>
#include <iostream>
using namespace std;
#define mod 998244353
#define MAXN 300005
#define int long long
int n, pi, ni, inv;
int a[MAXN], b[MAXN], c[MAXN], rev[MAXN], tmp[MAXN];int qkpow ( int x, int y ) {int ans = 1;while ( y ) {if ( y & 1 )ans = ans * x % mod;x = x * x % mod;y >>= 1;}return ans;
}void NTT ( int *c, int f, int len ) {for ( int i = 0;i < len;i ++ )if ( i < rev[i] )swap ( c[i], c[rev[i]] );for ( int i = 1;i < len;i <<= 1 ) {int omega = qkpow ( f == 1 ? pi : ni, ( mod - 1 ) / ( i << 1 ) );for ( int j = 0;j < len;j += ( i << 1 ) ) {int w = 1;for ( int k = 0;k < i;k ++, w = w * omega % mod ) {int x = c[k + j], y = w * c[k + j + i] % mod;c[k + j] = ( x + y ) % mod;c[k + j + i] = ( x - y + mod ) % mod;}}}inv = qkpow ( len, mod - 2 );if ( f == -1 )for ( int i = 0;i < len;i ++ )c[i] = c[i] * inv % mod;
}void polyinv ( int deg, int *a, int *b ) {if ( deg == 1 ) {b[0] = qkpow ( a[0], mod - 2 );return;}polyinv ( ( deg + 1 ) >> 1, a, b );int len = 1, l = 0;while ( len < deg * 2 - 1 ) {len <<= 1;l ++;}for ( int i = 0;i < len;i ++ )rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );for ( int i = 0;i < deg;i ++ )c[i] = a[i];for ( int i = deg;i < len;i ++ )c[i] = 0;NTT ( c, 1, len );NTT ( b, 1, len );for ( int i = 0;i < len;i ++ )b[i] = ( 2 - b[i] * c[i] % mod + mod ) % mod * b[i] % mod; NTT ( b, -1, len );for ( int i = deg;i < len;i ++ )b[i] = 0;
}void polyln ( int n, int *a, int *b ) {polyinv ( n, a, b );for ( int i = 1;i < n;i ++ )tmp[i - 1] = a[i] * i % mod;int len = 1, l = 0;while ( len < n * 2 - 1 ) {len <<= 1;l ++;}for ( int i = 0;i < len;i ++ )rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );NTT ( tmp, 1, len );NTT ( b, 1, len );for ( int i = 0;i < len;i ++ )b[i] = tmp[i] * b[i] % mod;NTT ( b, -1, len );for ( int i = n - 1;i; i -- )b[i] = b[i - 1] * qkpow ( i, mod - 2 ) % mod;b[0] = 0; for ( int i = n;i < len;i ++ )b[i] = 0;
} signed main() {pi = 3;ni = mod / pi + 1;scanf ( "%lld", &n );for ( int i = 0;i < n;i ++ )scanf ( "%lld", &a[i] );polyln ( n, a, b );for ( int i = 0;i < n;i ++ )printf ( "%lld ", b[i] );return 0;
}
多项式开根sqrt
理论推导
对于一个给定的多项式A(x)A(x)A(x),求B(x)B(x)B(x)满足
B2(x)≡A(x)(modxn)B^2(x)≡A(x) (mod\ x^n)B2(x)≡A(x)(mod xn)
如果n=1n=1n=1就直接是常数项开方
1.1.1.可以直接套牛顿迭代
B(x)=B0−B02(x)−A(x)2B0(x)=B0(x)+A(x)B0(x)2B(x)=B_0-\frac{B_0^2(x)-A(x)}{2B_0(x)}=\frac{B_0(x)+\frac{A(x)}{B_0(x)}}{2}B(x)=B0−2B0(x)B02(x)−A(x)=2B0(x)+B0(x)A(x)
2.2.2.也可以用逆元平方的思想
设B′(x)B'(x)B′(x)满足B′2(x)≡A(x)(modx⌈n2⌉)B'^2(x)\equiv A(x)(mod\ x^{\lceil \frac{n}{2}\rceil})B′2(x)≡A(x)(mod x⌈2n⌉),且我们已经知道它了
两式相减B2(x)−B′2(x)≡0(modx⌈n2⌉)B^2(x)-B'^2(x)\equiv 0(mod\ x^{\lceil \frac{n}{2}\rceil})B2(x)−B′2(x)≡0(mod x⌈2n⌉)
因式分解(B(x)−B′(x))(B(x)+B′(x)≡0(modx⌈n2⌉)(B(x)-B'(x))(B(x)+B'(x)\equiv 0(mod\ x^{\lceil \frac{n}{2}\rceil})(B(x)−B′(x))(B(x)+B′(x)≡0(mod x⌈2n⌉)
肯定只能是B(x)−B′(x)≡0(modx⌈n2⌉)B(x)-B'(x)\equiv 0(mod\ x^{\lceil \frac{n}{2}\rceil})B(x)−B′(x)≡0(mod x⌈2n⌉)不可能非负数系数开个根还能开出负数来,那你可以被飞签高等学府了
故伎重演,用求逆元的方式平方来一下,这里就直接打开了B2(x)+B′2(x)−2∗B(x)∗B′(x)≡0(modxn)B^2(x)+B'^2(x)-2*B(x)*B'(x)\equiv 0(mod\ x^n)B2(x)+B′2(x)−2∗B(x)∗B′(x)≡0(mod xn)
结合B2(x)B^2(x)B2(x)需要满足的条件,替换成A(x)A(x)A(x)
A(x)+B′2(x)−2∗B(x)∗B′(x)≡0(modxn)A(x)+B'^2(x)-2*B(x)*B'(x)\equiv 0(mod\ x^n)A(x)+B′2(x)−2∗B(x)∗B′(x)≡0(mod xn)
用B′(x),A(x)B'(x),A(x)B′(x),A(x)去表示B(x)B(x)B(x)
B(x)≡A(x)+B′(x)22∗B′(x)B(x)\equiv \frac{A(x)+B'(x)^2}{2*B'(x)}B(x)≡2∗B′(x)A(x)+B′(x)2
最后发现是孪生兄弟, 御用递归求解O(nlogn)O(nlogn)O(nlogn)
同时发现,只要常数项是二次项剩余且有逆元,这个多项式就可以开方
——摘自某dalao
模板
代码中常数项刚好是111,不是加强版,我太弱了
B(x)≡A(x)B′(x)+B′(x)2B(x)\equiv \frac{\frac{A(x)}{B'(x)}+B'(x)}{2}B(x)≡2B′(x)A(x)+B′(x)
void polysqrt ( int n, int *a, int *b ) {if ( n == 1 ) {b[0] = 1;return;}polysqrt ( ( n + 1 ) >> 1, a, b );static int tmp[MAXN], binv[MAXN];polyinv ( n, b, binv );//b就是公式的B'(x) tmp就是公式的A(x)int len = 1, l = 0;while ( len < n * 2 - 1 ) {len <<= 1;l ++;}for ( int i = 0;i < len;i ++ )rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );for ( int i = 0;i < n;i ++ )tmp[i] = a[i];for ( int i = n;i < len;i ++ )tmp[i] = 0;NTT ( tmp, 1, len );NTT ( b, 1, len );NTT ( binv, 1, len );for ( int i = 0;i < len;i ++ )//求解B(x) 用b表示b[i] = ( tmp[i] * binv[i] % mod + b[i] ) % mod * inv2 % mod;NTT ( b, -1, len );for ( int i = n;i < len;i ++ )b[i] = binv[i] = 0;for ( int i = 0;i < n;i ++ )binv[i] = 0;
}
例题
题目
luogu
给定一个n−1n-1n−1次多项式A(x)A(x)A(x),求一个在modxn\bmod\ x^nmod xn
意义下的多项式B(x)B(x)B(x),使得B2(x)≡A(x)(modxn)B^2(x) \equiv A(x) \ (\bmod\ x^n)B2(x)≡A(x) (mod xn)若有多解,请取零次项系数较小的作为答案。
多项式的系数在mod998244353\bmod\ 998244353mod 998244353的意义下进行运算。
输入格式
第一行一个正整数nnn
接下来nnn个整数,依次表示多项式的系数a0,a1,…,an−1a_0, a_1, \dots, a_{n-1}a0,a1,…,an−1
输出格式
输出nnn个整数,表示答案多项式的系数b0,b1,…,bn−1b_0, b_1, \dots, b_{n-1}b0,b1,…,bn−1
输入输出样例
输入
3
1 2 1
输出
1 1 0
输入
7
1 8596489 489489 4894 1564 489 35789489
输出
1 503420421 924499237 13354513 217017417 707895465 411020414
说明/提示
对于100%100\%100%的数据:n≤105ai∈[0,998244352]∩Zn \leq 10^5 \qquad a_i \in [0,998244352] \cap \mathbb{Z}n≤105ai∈[0,998244352]∩Z
code
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
#define mod 998244353
#define MAXN 300005
#define int long long
int pi, ni, inv, inv2;
int rev[MAXN];int qkpow ( int x, int y ) {int ans = 1;while ( y ) {if ( y & 1 )ans = ans * x % mod;x = x * x % mod;y >>= 1;}return ans;
}void NTT ( int *c, int f, int len ) {for ( int i = 0;i < len;i ++ )if ( i < rev[i] )swap ( c[i], c[rev[i]] );for ( int i = 1;i < len;i <<= 1 ) {int omega = qkpow ( f == 1 ? pi : ni, ( mod - 1 ) / ( i << 1 ) );for ( int j = 0;j < len;j += ( i << 1 ) ) {int w = 1;for ( int k = 0;k < i;k ++, w = w * omega % mod ) {int x = c[k + j], y = w * c[k + j + i] % mod;c[k + j] = ( x + y ) % mod;c[k + j + i] = ( x - y + mod ) % mod;}}}inv = qkpow ( len, mod - 2 );if ( f == -1 )for ( int i = 0;i < len;i ++ )c[i] = c[i] * inv % mod;
}void polyinv ( int deg, int *a, int *b ) {static int c[MAXN];if ( deg == 1 ) {b[0] = qkpow ( a[0], mod - 2 );return;}polyinv ( ( deg + 1 ) >> 1, a, b );int len = 1, l = 0;while ( len < deg * 2 - 1 ) {len <<= 1;l ++;}for ( int i = 0;i < len;i ++ )rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );for ( int i = 0;i < deg;i ++ )c[i] = a[i];for ( int i = deg;i < len;i ++ )c[i] = 0;NTT ( c, 1, len );NTT ( b, 1, len );for ( int i = 0;i < len;i ++ )b[i] = ( 2 - b[i] * c[i] % mod + mod ) % mod * b[i] % mod; NTT ( b, -1, len );for ( int i = deg;i < len;i ++ )b[i] = 0;
}void polysqrt ( int n, int *a, int *b ) {if ( n == 1 ) {b[0] = 1;return;}polysqrt ( ( n + 1 ) >> 1, a, b );static int tmp[MAXN], binv[MAXN];polyinv ( n, b, binv );int len = 1, l = 0;while ( len < n * 2 - 1 ) {len <<= 1;l ++;}for ( int i = 0;i < len;i ++ )rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );for ( int i = 0;i < n;i ++ )tmp[i] = a[i];for ( int i = n;i < len;i ++ )tmp[i] = 0;NTT ( tmp, 1, len );NTT ( b, 1, len );NTT ( binv, 1, len );for ( int i = 0;i < len;i ++ )b[i] = ( tmp[i] * binv[i] % mod + b[i] ) % mod * inv2 % mod;NTT ( b, -1, len );for ( int i = n;i < len;i ++ )b[i] = binv[i] = 0;for ( int i = 0;i < n;i ++ )binv[i] = 0;
}int a[MAXN], b[MAXN];signed main() {pi = 3;ni = mod / pi + 1;inv2 = qkpow ( 2, mod - 2 );int n;scanf ( "%lld", &n );for ( int i = 0;i < n;i ++ )scanf ( "%lld", &a[i] );polysqrt ( n, a, b );for ( int i = 0;i < n;i ++ )printf ( "%lld ", b[i] );return 0;
}
多项式指数exp
理论推导
对于一个给定的多项式A(x)A(x)A(x),求B(x)B(x)B(x)满足B(x)≡eA(x)(modxn)B(x) \equiv \text e^{A(x)}(mod\ x^n)B(x)≡eA(x)(mod xn)
取个对数,移个项ln(B(x))≡A(x)(modxn)——>ln(B(x))−A(x)≡0(modxn)ln(B(x))\equiv A(x)(mod\ x^n)——>ln(B(x))-A(x)\equiv 0(mod\ x^n)ln(B(x))≡A(x)(mod xn)——>ln(B(x))−A(x)≡0(mod xn)
对数求导公式:(lnx)′=1x(ln\ x)'=\frac{1}{x}(ln x)′=x1
别问我为什么,我也证不来,求教
接着套牛顿迭代B(x)≡B0(x)−g(B0(x))g′(B0(x))≡B0(x)−g(B0(x))1B0≡B0(x)−lnB0(x)−A(x)1B0(x)B(x)\equiv B_0(x)-\frac{g(B_0(x))}{g'(B_0(x))}\equiv B_0(x)-\frac{g(B_0(x))}{\frac{1}{B_0}}\equiv B_0(x)-\frac{lnB_0(x)-A(x)}{\frac{1}{B_0(x)}}B(x)≡B0(x)−g′(B0(x))g(B0(x))≡B0(x)−B01g(B0(x))≡B0(x)−B0(x)1lnB0(x)−A(x)≡B0(x)−B0(x)(lnB0(x)−A(x))≡B0(x)(1−lnB0(x)+A(x))(modx2n)\equiv B_0(x)-B_0(x)(lnB_0(x)-A(x))\equiv B_0(x)(1-lnB_0(x)+A(x))\ (mod\ x^{2n})≡B0(x)−B0(x)(lnB0(x)−A(x))≡B0(x)(1−lnB0(x)+A(x)) (mod x2n)
必须保证A(x)常数项为0,B(x)常数项为1,不知道为什么,求dalao教
模板
void polyexp ( int n, int *a, int *b ) {if ( n == 1 ) {b[0] = 1;return;}polyexp ( ( n + 1 ) >> 1, a, b );static int bln[MAXN];polyln ( n, b, bln );for ( int i = 0;i < n;i ++ )bln[i] = ( i == 0 ) + a[i] - bln[i] + mod;int len = 1, l = 0;while ( len < n * 2 - 1 ) {len <<= 1;l ++;}for ( int i = 0;i < len;i ++ )rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );NTT ( bln, 1, len );NTT ( b, 1, len );for ( int i = 0;i < len;i ++ )b[i] = b[i] * bln[i] % mod;NTT ( b, -1, len );for ( int i = n;i < len;i ++ )b[i] = bln[i] = 0;for ( int i = 0;i < n;i ++ )bln[i] = 0;
}
例题
luogu
给出 n−1n-1n−1 次多项式 A(x)A(x)A(x),求一个 modxn\bmod{\:x^n}modxn下的多项式 B(x)B(x)B(x),满足 B(x)≡eA(x)B(x) \equiv \text e^{A(x)}B(x)≡eA(x),系数对 998244353998244353998244353 取模
输入格式
第一行一个整数 nnn
下一行有 nnn 个整数,依次表示多项式的系数 a0,a1,⋯,an−1a_0, a_1, \cdots, a_{n-1}a0,a1,⋯,an−1
输出格式
输出 nnn 个整数,表示答案多项式中的系数 a0,a1,⋯,an−1a_0, a_1, \cdots, a_{n-1}a0,a1,⋯,an−1
输入输出样例
输入
6
0 927384623 817976920 427326948 149643566 610586717
输出
1 927384623 878326372 3882 273455637 998233543
说明/提示
对于 100%100\%100% 的数据,n≤105n \le 10^5n≤105
题目
code
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
#define mod 998244353
#define MAXN 300005
#define int long long
int pi, ni, inv;
int rev[MAXN];int qkpow ( int x, int y ) {int ans = 1;while ( y ) {if ( y & 1 )ans = ans * x % mod;x = x * x % mod;y >>= 1;}return ans;
}void NTT ( int *c, int f, int len ) {for ( int i = 0;i < len;i ++ )if ( i < rev[i] )swap ( c[i], c[rev[i]] );for ( int i = 1;i < len;i <<= 1 ) {int omega = qkpow ( f == 1 ? pi : ni, ( mod - 1 ) / ( i << 1 ) );for ( int j = 0;j < len;j += ( i << 1 ) ) {int w = 1;for ( int k = 0;k < i;k ++, w = w * omega % mod ) {int x = c[k + j], y = w * c[k + j + i] % mod;c[k + j] = ( x + y ) % mod;c[k + j + i] = ( x - y + mod ) % mod;}}}inv = qkpow ( len, mod - 2 );if ( f == -1 )for ( int i = 0;i < len;i ++ )c[i] = c[i] * inv % mod;
}void polyinv ( int deg, int *a, int *b ) {static int c[MAXN];if ( deg == 1 ) {b[0] = qkpow ( a[0], mod - 2 );return;}polyinv ( ( deg + 1 ) >> 1, a, b );int len = 1, l = 0;while ( len < deg * 2 - 1 ) {len <<= 1;l ++;}for ( int i = 0;i < len;i ++ )rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );for ( int i = 0;i < deg;i ++ )c[i] = a[i];for ( int i = deg;i < len;i ++ )c[i] = 0;NTT ( c, 1, len );NTT ( b, 1, len );for ( int i = 0;i < len;i ++ )b[i] = ( 2 - b[i] * c[i] % mod + mod ) % mod * b[i] % mod; NTT ( b, -1, len );for ( int i = deg;i < len;i ++ )b[i] = 0;
}void polyln ( int n, int *a, int *b ) {static int tmp[MAXN];polyinv ( n, a, b );for ( int i = 1;i < n;i ++ )tmp[i - 1] = a[i] * i % mod;int len = 1, l = 0;while ( len < n * 2 - 1 ) {len <<= 1;l ++;}for ( int i = 0;i < len;i ++ )rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );NTT ( tmp, 1, len );NTT ( b, 1, len );for ( int i = 0;i < len;i ++ )b[i] = tmp[i] * b[i] % mod;NTT ( b, -1, len );for ( int i = n - 1;i; i -- )b[i] = b[i - 1] * qkpow ( i, mod - 2 ) % mod;b[0] = 0;for ( int i = n;i < len;i ++ )b[i] = 0;
} void polyexp ( int n, int *a, int *b ) {if ( n == 1 ) {b[0] = 1;return;}polyexp ( ( n + 1 ) >> 1, a, b );static int bln[MAXN];polyln ( n, b, bln );for ( int i = 0;i < n;i ++ )bln[i] = ( i == 0 ) + a[i] - bln[i] + mod;int len = 1, l = 0;while ( len < n * 2 - 1 ) {len <<= 1;l ++;}for ( int i = 0;i < len;i ++ )rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );NTT ( bln, 1, len );NTT ( b, 1, len );for ( int i = 0;i < len;i ++ )b[i] = b[i] * bln[i] % mod;NTT ( b, -1, len );for ( int i = n;i < len;i ++ )b[i] = bln[i] = 0;for ( int i = 0;i < n;i ++ )bln[i] = 0;
}int a[MAXN], b[MAXN];signed main() {pi = 3;ni = mod / pi + 1;int n;scanf ( "%lld", &n );for ( int i = 0;i < n;i ++ )scanf ( "%lld", &a[i] );polyexp ( n, a, b );for ( int i = 0;i < n;i ++ )printf ( "%lld ", b[i] );return 0;
}
多项式快速幂
理论推导
通过上述多项式的基本操作后,我们直接暴力来一发拆公式
B(x)≡Ak(x)(modxn)B(x)≡A^k(x)\ (mod x^n)B(x)≡Ak(x) (modxn)
==>==>==> B(x)≡elnA(x)kB(x)\equiv e^{lnA(x)^k}B(x)≡elnA(x)k
==>==>==> B(x)≡ek∗lnA(x)B(x)\equiv e^{k*lnA(x)}B(x)≡ek∗lnA(x)
所以我们只用上述所讲的ln,expln,expln,exp即可完成快速幂
例题
题目
luogu
给定一个n−1n−1n−1次多项式A(x)A(x)A(x),求一个在modxn\bmod\ x^nmod xn
意义下的多项式B(x)B(x)B(x),使得B(x)≡Ak(x)(modxn)B(x) \equiv A^k(x) \ (\bmod\ x^n)B(x)≡Ak(x) (mod xn)
多项式的系数在mod998244353\bmod\ 998244353mod 998244353的意义下进行运算。
输入格式
第一行两个正整数n,kn,kn,k。
接下来nnn个整数,依次表示多项式的系数a0,a1,…,an−1a_0, a_1, \dots, a_{n-1}a0,a1,…,an−1
输出格式
输出nnn个整数,表示答案多项式的系数b0,b1,…,bn−1b_0, b_1, \dots, b_{n-1}b0,b1,…,bn−1
输入输出样例
输入
9 18948465
1 2 3 4 5 6 7 8 9
输出
1 37896930 597086012 720637306 161940419 360472177 560327751 446560856 524295016
输入
4 1
1 1 0 0
输出
1 1 0 0
输入
4 2
1 1 0 0
输出
1 2 1 0
输入
4 3
1 1 0 0
输出
1 3 3 1
说明/提示
对于100%100\%100%的数据:n≤1052≤k≤10105ai∈[0,998244352]∩Zn \leq 10^5 \qquad 2 \leq k \leq 10^{10^5}\qquad a_i \in [0,998244352] \cap \mathbb{Z}n≤1052≤k≤10105ai∈[0,998244352]∩Z
code
唯一的坑点🕳就是注意kkk很大,longlonglong\ longlong long装不下,必须快读取模
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
#define mod 998244353
#define MAXN 300005
#define int long long
int pi, ni, inv;
int rev[MAXN], result[MAXN];int qkpow ( int x, int y ) {int ans = 1;while ( y ) {if ( y & 1 )ans = ans * x % mod;x = x * x % mod;y >>= 1;}return ans;
}void NTT ( int *c, int f, int len ) {for ( int i = 0;i < len;i ++ )if ( i < rev[i] )swap ( c[i], c[rev[i]] );for ( int i = 1;i < len;i <<= 1 ) {int omega = qkpow ( f == 1 ? pi : ni, ( mod - 1 ) / ( i << 1 ) );for ( int j = 0;j < len;j += ( i << 1 ) ) {int w = 1;for ( int k = 0;k < i;k ++, w = w * omega % mod ) {int x = c[k + j], y = w * c[k + j + i] % mod;c[k + j] = ( x + y ) % mod;c[k + j + i] = ( x - y + mod ) % mod;}}}inv = qkpow ( len, mod - 2 );if ( f == -1 )for ( int i = 0;i < len;i ++ )c[i] = c[i] * inv % mod;
}void polyinv ( int deg, int *a, int *b ) {static int c[MAXN];if ( deg == 1 ) {b[0] = qkpow ( a[0], mod - 2 );return;}polyinv ( ( deg + 1 ) >> 1, a, b );int len = 1, l = 0;while ( len < deg * 2 - 1 ) {len <<= 1;l ++;}for ( int i = 0;i < len;i ++ )rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );for ( int i = 0;i < deg;i ++ )c[i] = a[i];for ( int i = deg;i < len;i ++ )c[i] = 0;NTT ( c, 1, len );NTT ( b, 1, len );for ( int i = 0;i < len;i ++ )b[i] = ( 2 - b[i] * c[i] % mod + mod ) % mod * b[i] % mod; NTT ( b, -1, len );for ( int i = deg;i < len;i ++ )b[i] = 0;
}void polyln ( int n, int *a, int *b ) {static int tmp[MAXN];polyinv ( n, a, b );for ( int i = 1;i < n;i ++ )tmp[i - 1] = a[i] * i % mod;int len = 1, l = 0;while ( len < n * 2 - 1 ) {len <<= 1;l ++;}for ( int i = 0;i < len;i ++ )rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );NTT ( tmp, 1, len );NTT ( b, 1, len );for ( int i = 0;i < len;i ++ )b[i] = tmp[i] * b[i] % mod;NTT ( b, -1, len );for ( int i = n - 1;i; i -- )b[i] = b[i - 1] * qkpow ( i, mod - 2 ) % mod;b[0] = 0;for ( int i = n;i < len;i ++ )b[i] = tmp[i] = 0;for ( int i = 0;i < n;i ++ )tmp[i] = 0;
} void polyexp ( int n, int *a, int *b ) {if ( n == 1 ) {b[0] = 1;return;}polyexp ( ( n + 1 ) >> 1, a, b );static int bln[MAXN];polyln ( n, b, bln );for ( int i = 0;i < n;i ++ )bln[i] = ( ( i == 0 ) + a[i] - bln[i] + mod ) % mod;int len = 1, l = 0;while ( len < n * 2 - 1 ) {len <<= 1;l ++;}for ( int i = 0;i < len;i ++ )rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );NTT ( bln, 1, len );NTT ( b, 1, len );for ( int i = 0;i < len;i ++ )b[i] = b[i] * bln[i] % mod;NTT ( b, -1, len );for ( int i = n;i < len;i ++ )b[i] = bln[i] = 0;for ( int i = 0;i < n;i ++ )bln[i] = 0;
}int a[MAXN], b[MAXN];void read ( int &x ) {x = 0;char s = getchar();while ( s < '0' || s > '9' )s = getchar();while ( s >= '0' && s <= '9' ) {x = ( ( x << 1 ) + ( x << 3 ) + ( s - '0' ) ) % mod;s = getchar();}
} signed main() {pi = 3;ni = mod / pi + 1;int n, k;scanf ( "%lld", &n );read ( k );//k太大了,必须要取模long long装不下 for ( int i = 0;i < n;i ++ )scanf ( "%lld", &a[i] );polyln ( n, a, b );for ( int i = 0;i < n;i ++ )b[i] = b[i] * k % mod;polyexp ( n, b, result );for ( int i = 0;i < n;i ++ )printf ( "%lld ", result[i] );return 0;
}
版权申明:
教会我多项式的blog
证明无代码
最后只想说一句,希望各个dalao能完善我的证明,甚至教教这个蒟蒻,部分实在证不来,百度百科长的又丑有问题欢迎评论,受教了