众所周知,多项式全家桶指的是
写在前面
- 为了进一步优化,请区分 int 和 long long,尽量 int
NTT
虚单位根常熟大,考虑替代
模意义下原根可以替代
一般取mod=998244353,原根g=3
里面的一个操作(reverse)需要群论 证明,gun
Tips
- 由费马小定理 a p − 1 ≡ 1 ( m o d p ) a^{p-1}\equiv 1\pmod p ap−1≡1(modp)得任意 a a a在模意义下的逆元是 a p − 2 a^{p-2} ap−2
#include<bits/stdc++.h>
using namespace std;
#define in Read()
#define int long long
inline int in{
int i=0,f=1;char ch=0;
while(ch!='-'&&!isdigit(ch)) ch=getchar();
if(ch=='-') ch=getchar(),f=-1;
while(isdigit(ch)) i=(i<<1)+(i<<3)+ch-48,ch=getchar();
return i*f;
}
const int NNN=5e6+5;
const int MOD=998244353;
const int G=3;
int n,m,rev[NNN];
int A[NNN],B[NNN],C[NNN];
inline int qpow(int a,int x){
int res=1;
while(x){
if(x&1) res=res*a%MOD;
a=a*a%MOD;
x>>=1;
}
return res;
}
inline void NTT(int *f,int lim,int sgn){
for(int i=0;i<lim;++i)
if(i<rev[i]) swap(f[i],f[rev[i]]);
for(int len=1;len<lim;len<<=1){
int siz=len<<1;
int wn=qpow(G,(MOD-1)/siz);
for(int l=0;l<lim;l+=siz){
int w=1;
for(int i=l;i<l+len;++i){
int a=f[i],b=f[i+len]*w%MOD;
f[i]=(a+b)%MOD;
f[i+len]=(a-b+MOD)%MOD;
w=w*wn%MOD;
}
}
}
if(sgn) return;
int inv=qpow(lim,MOD-2);//inverse element of 'limit'
reverse(f+1,f+lim);//Just recite it babe
for(int i=0;i<lim;++i)
f[i]=f[i]*inv%MOD;
}
signed main(){
n=in,m=in;
for(int i=0;i<=n;++i) A[i]=in;
for(int i=0;i<=m;++i) B[i]=in;
int lim=1;
while(lim<=n+m) lim<<=1;
for(int i=0;i<lim;++i) rev[i]=(rev[i>>1]>>1)|((i&1)?lim>>1:0);
NTT(A,lim,1);
NTT(B,lim,1);
for(int i=0;i<lim;++i) C[i]=A[i]*B[i]%MOD;
// for(int i=1;i<=lim;++i)
// printf("%lld %lld\n",A[i],B[i]);
NTT(C,lim,0);
for(int i=0;i<=n+m;++i) printf("%lld ",C[i]);
return 0;
}
求逆
逆元的定义是模意义下,与本原积为1的东西,多项式这里当然就是多项式啦
A
(
x
)
B
(
x
)
≡
1
(
m
o
d
p
(
u
s
u
a
l
l
y
x
n
)
)
A(x)B(x)\equiv 1(\bmod p(\mathrm{usually}\ x^n))
A(x)B(x)≡1(modp(usually xn))
则称
B
(
x
)
B(x)
B(x)为
A
(
x
)
A(x)
A(x)的逆元,记为
A
−
1
(
x
)
A^{-1}(x)
A−1(x)
一顿暴搞得到
m
o
d
x
⌈
n
2
⌉
\bmod x^{\lceil\frac{n}{2}\rceil}
modx⌈2n⌉意义下的
f
0
−
1
f_0^{-1}
f0−1变为
m
o
d
x
n
\bmod x^{n}
modxn意义下的
f
−
1
f^{-1}
f−1的递归式是
f
−
1
≡
f
0
−
1
(
2
−
f
f
0
−
1
)
(
m
o
d
x
n
)
f^{-1}\equiv f_0^{-1}(2-ff_0^{-1}) \pmod {x^n}
f−1≡f0−1(2−ff0−1)(modxn)
这顿暴搞:
f
−
1
−
f
0
−
1
≡
0
(
m
o
d
x
⌈
n
2
⌉
)
平
方
f
−
2
−
2
f
−
1
f
0
−
1
+
f
0
−
2
≡
0
(
m
o
d
x
n
)
移
项
f
−
2
≡
2
f
−
1
f
0
−
1
−
f
0
−
2
(
m
o
d
x
n
)
同
乘
f
f
−
1
≡
2
f
0
−
1
−
f
f
0
−
2
f
−
1
≡
f
0
−
1
(
2
−
f
f
0
−
1
)
f^{-1}-f_0^{-1}\equiv0\pmod{x^{\lceil\frac{n}{2}\rceil}}\\ 平方f^{-2}-2f^{-1}f_0^{-1}+f_0^{-2}\equiv 0 \pmod{x^n}\\ 移项f^{-2}\equiv2f^{-1}f_0^{-1}-f_0^{-2}\pmod{x^n}\\ 同乘f\quad f^{-1}\equiv2f_0^{-1}-ff_0^{-2}\\ f^{-1}\equiv f_0^{-1}(2-ff_0^{-1})
f−1−f0−1≡0(modx⌈2n⌉)平方f−2−2f−1f0−1+f0−2≡0(modxn)移项f−2≡2f−1f0−1−f0−2(modxn)同乘ff−1≡2f0−1−ff0−2f−1≡f0−1(2−ff0−1)
操作就递归,时间复杂度就 O ( n log n ) O(n\log n) O(nlogn)
Tips
- 务必保护原多项式
- 在 x ⌈ λ ⌉ x^{\lceil \lambda\rceil} x⌈λ⌉以外的全没有,务必赋零
- 为了防止爆 int 之类的东西,请在乘的时候多乘个 1ll
- 求逆的时候有个数组拿来临时存装原多项式的信息,如果在函数里面定义它就炸了,定义在外面就没问题(由于递归的性质,定义在外面是不影响正确性的)
#include<bits/stdc++.h>
using namespace std;
#define in Read()
typedef long long ll;
inline int in{
int i=0,f=1;char ch=0;
while(ch!='-'&&!isdigit(ch)) ch=getchar();
if(ch=='-') ch=getchar(),f=-1;
while(isdigit(ch)) i=(i<<1)+(i<<3)+ch-48,ch=getchar();
return i*f;
}
const int NNN=5e5+5;
const ll MOD=998244353;
int n,rev[NNN];
ll A[NNN],I[NNN],h[NNN];
inline ll qpow(ll a,int x){
ll res=1;
while(x){
if(x&1) res=res*a%MOD;
a=a*a%MOD;
x>>=1;
}
return res;
}
inline void NTT(ll *f,int lim,int sgn){
for(int i=0;i<lim;++i)
if(i<rev[i]) swap(f[i],f[rev[i]]);
for(int len=1;len<lim;len<<=1){
int siz=len<<1;
ll wn=qpow(3,(MOD-1)/siz);
for(int l=0;l<lim;l+=siz){
ll w=1;
for(int i=l;i<l+len;++i){
ll a=f[i],b=f[i+len]*w%MOD;
f[i]=(a+b)%MOD;
f[i+len]=(a-b+MOD)%MOD;
w=w*wn%MOD;
}
}
}
if(sgn) return;
int inv=qpow(lim,MOD-2);
reverse(f+1,f+lim);
for(int i=0;i<lim;++i)
f[i]=f[i]*inv%MOD;
}
inline void Inv(int len,ll *f,ll *g){
if(len==1){
g[0]=qpow(f[0],MOD-2);
return;
}
Inv((len+1)>>1,f,g);//ceil len
int lim=1;
while(lim<len) lim<<=1;lim<<=1;
for(int i=0;i<lim;++i) rev[i]=(rev[i>>1]>>1)|((i&1)?lim>>1:0);
for(int i=0;i<len;++i) h[i]=f[i];
for(int i=len;i<lim;++i) h[i]=0;
NTT(h,lim,1);
NTT(g,lim,1);
for(int i=0;i<lim;++i) g[i]=1ll*(1ll*2-1ll*g[i]*h[i]%MOD+MOD)%MOD*g[i]%MOD;
NTT(g,lim,0);
for(int i=len;i<lim;++i) g[i]=0;
return;
}
int main(){
n=in;
for(int i=0;i<n;++i) A[i]=in;
Inv(n,A,I);
for(int i=0;i<n;++i) printf("%d ",I[i]);
return 0;
}
ln
求个导再积个分
U
n
d
e
r
(
m
o
d
x
n
)
B
=
ln
A
⇒
B
′
=
(
ln
A
)
′
×
A
′
=
A
′
A
⇒
B
=
∫
B
′
d
x
=
∫
A
′
A
d
x
\begin{aligned} &Under\pmod{x^{n}}\\ &B=\ln A\\ \Rightarrow&B'=(\ln A)'\times A'=\frac{A'}{A}\\ \Rightarrow&B=\int B'\mathrm{d}x=\int \frac{A'}{A}\mathrm{d}x \end{aligned}
⇒⇒Under(modxn)B=lnAB′=(lnA)′×A′=AA′B=∫B′dx=∫AA′dx
求导,积分都是 O ( n ) O(n) O(n),求逆(求 1 A \frac{1}{A} A1)乘法都是 O ( n log n ) O(n\log n) O(nlogn)
没有除法!
直接
i
n
v
≡
i
p
−
2
inv\equiv i^{p-2}
inv≡ip−2复杂度过高
那么需要逆元,介绍一种方法:
对
于
{
a
n
}
(
通
常
是
a
n
=
n
的
数
列
)
,
记
p
r
i
=
∏
j
=
1
i
a
j
,
第
一
遍
,
i
n
v
n
=
p
r
n
p
−
2
(
即
p
r
n
的
逆
元
)
≡
1
∏
i
=
1
n
a
i
i
n
v
i
=
i
n
v
i
+
1
×
a
i
+
1
,
递
推
下
去
第
二
遍
,
i
n
v
i
=
i
n
v
i
×
p
r
i
−
1
时
间
复
杂
度
O
(
N
)
\begin{aligned} &对于\{a_n\}(通常是a_n=n的数列),记pr_i=\prod_{j=1}^ia_j,\\ &第一遍,inv_n=pr_n^{p-2}(即pr_n的逆元)\equiv\frac{1}{\prod_{i=1}^na_i}\\ &inv_i=inv_{i+1}\times a_{i+1},递推下去\\ &第二遍,inv_i=inv_i\times pr_{i-1} \\\\&时间复杂度O(N) \end{aligned}
对于{an}(通常是an=n的数列),记pri=j=1∏iaj,第一遍,invn=prnp−2(即prn的逆元)≡∏i=1nai1invi=invi+1×ai+1,递推下去第二遍,invi=invi×pri−1时间复杂度O(N)
不过我似乎不会这种方法,那么再来这种好写的
模
p
意
义
下
,
求
i
的
逆
元
p
=
k
i
+
r
(
k
=
⌊
p
i
⌋
,
r
=
p
m
o
d
i
)
k
i
+
r
≡
0
(
m
o
d
p
)
i
−
1
r
−
1
(
k
i
+
r
)
≡
0
(
m
o
d
p
)
i
−
1
+
k
r
−
1
≡
0
(
m
o
d
p
)
i
−
1
≡
−
k
r
−
1
(
m
o
d
p
)
由
于
r
=
p
m
o
d
i
,
则
r
<
i
,
求
i
时
r
已
求
出
,
可
以
递
归
求
解
边
界
1
−
1
=
1
\begin{aligned} &模p意义下,求i的逆元\\ &p=ki+r\quad(k=\lfloor\frac{p}{i}\rfloor,r=p\bmod i)\\ &ki+r\equiv 0\pmod p\\ &i^{-1}r^{-1}(ki+r)\equiv 0\pmod p\\ &i^{-1}+kr^{-1}\equiv 0\pmod p\\ &i^{-1}\equiv -kr^{-1}\pmod p\\ &由于r=p\bmod i,则r<i,求i时r已求出,可以递归求解\\ &边界1^{-1}=1 \end{aligned}
模p意义下,求i的逆元p=ki+r(k=⌊ip⌋,r=pmodi)ki+r≡0(modp)i−1r−1(ki+r)≡0(modp)i−1+kr−1≡0(modp)i−1≡−kr−1(modp)由于r=pmodi,则r<i,求i时r已求出,可以递归求解边界1−1=1
#include<bits/stdc++.h>
using namespace std;
#define in Read()
typedef long long ll;
inline int in{
int i=0,f=1;char ch=0;
while(ch!='-'&&!isdigit(ch)) ch=getchar();
if(ch=='-') ch=getchar(),f=-1;
while(isdigit(ch)) i=(i<<1)+(i<<3)+ch-48,ch=getchar();
return i*f;
}
const int NNN=5e5+5;
const ll MOD=998244353;
int n,rev[NNN];
ll int_inv[NNN],prod[NNN];
ll A[NNN],B[NNN],C[NNN],D[NNN],E[NNN];
//B=lnA target polynomial
//C=A' differential coefficient
//D is inversion of A
//B is integral of A
//C<-A, D<-A, E=C*D, B<-E
inline ll qpow(ll a,int x){
ll res=1;
while(x){
if(x&1) res=res*a%MOD;
a=a*a%MOD;
x>>=1;
}
return res;
}
ll h[NNN];
inline void NTT(ll *f,int lim,int sgn){
for(int i=0;i<lim;++i)
if(i<rev[i]) swap(f[i],f[rev[i]]);
for(int len=1;len<lim;len<<=1){
int siz=len<<1;
ll wn=qpow(3,(MOD-1)/siz);
for(int l=0;l<lim;l+=siz){
ll w=1;
for(int i=l;i<l+len;++i){
ll a=f[i],b=f[i+len]*w%MOD;
f[i]=(a+b)%MOD;
f[i+len]=(a-b+MOD)%MOD;
w=w*wn%MOD;
}
}
}
if(sgn) return;
ll inv=qpow(lim,MOD-2);
reverse(f+1,f+lim);
for(int i=0;i<lim;++i)
f[i]=f[i]*inv%MOD;
return;
}
inline void Inv(int len,ll *f,ll *g){
if(len==1){
g[0]=qpow(f[0],MOD-2);
return;
}
Inv((len+1)>>1,f,g);
int lim=1;
while(lim<len) lim<<=1;lim<<=1;
for(int i=0;i<lim;++i) rev[i]=(rev[i>>1]>>1)|((i&1)?lim>>1:0);
for(int i=0;i<len;++i) h[i]=f[i];
for(int i=len;i<lim;++i) h[i]=0;
NTT(h,lim,1);
NTT(g,lim,1);
for(int i=0;i<lim;++i) g[i]=1ll*(2ll-1ll*g[i]*h[i]%MOD+MOD)%MOD*g[i]%MOD;
NTT(g,lim,0);
for(int i=len;i<lim;++i) g[i]=0;
return;
}
int main(){
n=in;
for(int i=0;i<n;++i) A[i]=in;
for(int i=0;i<n;++i) C[i]=A[i+1]*(i+1)%MOD;
Inv(n,A,D);
int lim=1;
while(lim<(n<<1)-1) lim<<=1;
for(int i=0;i<lim;++i) rev[i]=(rev[i>>1]>>1)|((i&1)?lim>>1:0);
NTT(C,lim,1);
NTT(D,lim,1);
for(int i=0;i<lim;++i) E[i]=C[i]*D[i];
NTT(E,lim,0);
int_inv[1]=1;
for(int i=2;i<lim;++i) int_inv[i]=(MOD-MOD/i)*int_inv[MOD%i]%MOD;
for(int i=0;i<lim;++i) B[i+1]=(E[i]%MOD+MOD)%MOD*int_inv[i+1]%MOD;
for(int i=0;i<n;++i) printf("%d ",B[i]);
}
牛顿迭代
心态爆炸,想学一个,发现前置知识是另一个,然后就递归了qwq
话说这不是拓扑排序吗qwqwq
已知
g
(
x
)
g(x)
g(x),求
f
(
x
)
f(x)
f(x),其中
g
(
f
(
x
)
)
≡
0
(
m
o
d
x
n
)
g(f(x))\equiv0\pmod {x^n}
g(f(x))≡0(modxn)
前置知识Talor’s Expansion(写这篇博客的时候本蒟蒻还很菜,当然现在也还很菜,大家将就着看)
牛顿迭代,考虑倍增
已知
g
(
f
0
)
≡
0
(
m
o
d
x
⌈
n
2
⌉
)
g(f_0)\equiv 0\pmod {x^{\lceil\frac{n}{2}\rceil}}
g(f0)≡0(modx⌈2n⌉)
要求
g
(
f
)
≡
0
(
m
o
d
x
n
)
g(f)\equiv 0\pmod{x^n}
g(f)≡0(modxn)
对
g
(
f
)
g(f)
g(f)在
f
0
f_0
f0处泰勒展开
∑
i
=
1
+
∞
g
(
i
)
(
f
0
)
i
!
(
f
−
f
0
)
i
≡
0
(
m
o
d
x
n
)
\sum_{i=1}^{+\infin}\frac{g^{(i)}(f_0)}{i!}(f-f_0)^i\equiv0\pmod{x^n}\\
i=1∑+∞i!g(i)(f0)(f−f0)i≡0(modxn)
容易知道
f
−
f
0
f-f_0
f−f0是
x
⌈
n
2
⌉
→
x
n
x^{\lceil\frac{n}{2}\rceil}\to x^n
x⌈2n⌉→xn次的,那么当
i
≥
2
i\geq 2
i≥2时,
(
f
−
f
0
)
i
≡
0
(
m
o
d
x
n
)
(f-f_0)^i\equiv0\pmod{x^n}
(f−f0)i≡0(modxn)
所以
g
(
f
0
)
+
g
′
(
f
0
)
⋅
(
f
−
f
0
)
≡
0
(
m
o
d
x
n
)
f
⋅
g
′
(
f
0
)
≡
g
′
(
f
0
)
⋅
f
0
−
g
(
f
0
)
(
m
o
d
x
n
)
f
≡
f
0
−
g
(
f
0
)
g
′
(
f
0
)
(
m
o
d
x
n
)
\begin{aligned} g(f_0)+g'(f_0)\cdot (f-f_0)&\equiv0\pmod{x^n}\\ f\cdot g'(f_0)&\equiv g'(f_0)\cdot f_0-g(f_0)\pmod{x^n}\\ f&\equiv f_0-\frac{g(f_0)}{g'(f_0)}\pmod{x^n} \end{aligned}
g(f0)+g′(f0)⋅(f−f0)f⋅g′(f0)f≡0(modxn)≡g′(f0)⋅f0−g(f0)(modxn)≡f0−g′(f0)g(f0)(modxn)
最下面的那个柿子是结论,建议背住
此外,边界情况依题而定
exp
已知
A
A
A,求
B
=
e
A
B=e^A
B=eA(保证
A
(
0
)
=
0
A(0)=0
A(0)=0,即
A
A
A的常数项为
0
0
0,因为如果不是这样你求出来的东西就不是整数了)
B
=
e
A
⇒
ln
B
=
A
ln
B
−
A
=
0
,
兴
奋
地
把
它
看
做
一
个
函
数
,
B
为
变
量
,
A
为
常
数
记
为
g
(
B
)
=
ln
B
−
A
=
0
(
即
求
零
点
,
于
是
考
虑
牛
顿
迭
代
)
B
=
B
0
−
g
(
B
0
)
g
′
(
B
0
)
g
′
(
B
)
=
d
g
d
B
(
特
别
写
出
来
,
注
意
变
量
)
=
1
B
=
B
−
1
B
=
B
0
−
ln
B
0
−
A
B
0
−
1
=
B
0
−
(
ln
B
0
−
A
)
B
0
=
B
0
(
1
+
ln
B
0
−
A
)
\begin{aligned} &B=e^A\Rightarrow\ln B=A\\ &\ln B-A=0,兴奋地把它看做一个函数,B为变量,A为常数\\ &记为g(B)=\ln B-A=0(即求零点,于是考虑牛顿迭代)\\ &B=B_0-\frac{g(B_0)}{g'(B_0)}\\ &g'(B)=\frac{\mathrm{d}g}{\mathrm{d}B}(特别写出来,注意变量)=\frac{1}{B}=B^{-1}\\ B&=B_0-\frac{\ln B_0-A}{B_0^{-1}}\\ &=B_0-(\ln B_0-A)B_0\\ &=B_0(1+\ln B_0-A) \end{aligned}
BB=eA⇒lnB=AlnB−A=0,兴奋地把它看做一个函数,B为变量,A为常数记为g(B)=lnB−A=0(即求零点,于是考虑牛顿迭代)B=B0−g′(B0)g(B0)g′(B)=dBdg(特别写出来,注意变量)=B1=B−1=B0−B0−1lnB0−A=B0−(lnB0−A)B0=B0(1+lnB0−A)
倍增做
边界:
n
=
0
时
,
B
=
e
0
=
1
n=0时, B=e^0=1
n=0时,B=e0=1
注意数组清零
#include<bits/stdc++.h>
using namespace std;
#define in Read()
typedef long long ll;
inline int in{
int i=0,f=1;char ch=0;
while(ch!='-'&&!isdigit(ch)) ch=getchar();
if(ch=='-') ch=getchar(),f=-1;
while(isdigit(ch)) i=(i<<1)+(i<<3)+ch-48,ch=getchar();
return i*f;
}
const int NNN=5e5+5;
const ll mod=998244353;
int rev[NNN],lim;
inline ll qpow(ll a,int x){
ll res=1;
while(x){
if(x&1) res=res*a%mod;
a=a*a%mod;
x>>=1;
}
return res;
}
inline void Init(int len){
lim=1;
while(lim<len) lim<<=1;
for(int i=0;i<lim;++i)
rev[i]=(rev[i>>1]>>1)|((i&1)?lim>>1:0);
return;
}
inline void NTT(ll *f,int sgn){
for(int i=0;i<lim;++i)
if(i<rev[i]) swap(f[i],f[rev[i]]);
for(int len=1;len<lim;len<<=1){
int siz=len<<1;
ll wn=qpow(3,(mod-1)/siz);
for(int l=0;l<lim;l+=siz){
ll w=1;
for(int i=l;i<l+len;++i){
ll a=f[i],b=f[i+len]*w%mod;
f[i]=(a+b)%mod;
f[i+len]=(a-b+mod)%mod;
w=w*wn%mod;
}
}
}
if(sgn) return;
int inv=qpow(lim,mod-2);
reverse(f+1,f+lim);
for(int i=0;i<lim;++i)
f[i]=f[i]*inv%mod;
return;
}
//After successfully debuged, come back and try to make it `inline ll *Con(...)`
inline void Con(ll *f,ll *g,int lf,int lg){//Convolve h=f*g a.k.a.卷积
Init(lf+lg);
NTT(f,1);
NTT(g,1);
for(int i=0;i<lim;++i) f[i]=f[i]*g[i]%mod;
NTT(f,0);
return;
}
ll tmp_inv[NNN];
inline void Inv(ll *f,ll *g,int len){//inversion: fg=1, f=>g
#define h tmp_inv
if(len==1){
g[0]=qpow(f[0],mod-2);
return;
}
Inv(f,g,(len+1)>>1);
Init(len<<1);
for(int i=0;i<len;++i) h[i]=f[i];
for(int i=len;i<lim;++i) h[i]=0;
NTT(g,1);
NTT(h,1);
for(int i=0;i<lim;++i) g[i]=1ll*(2ll-1ll*h[i]*g[i]%mod+mod)%mod*g[i]%mod;
NTT(g,0);
for(int i=len;i<lim;++i) g[i]=0;
return;
#undef h
}
inline void Der(ll *f,ll *g,int len){//derivation: g=f' a.k.a.求导
for(int i=0;i<len-1;++i) g[i]=f[i+1]*(i+1)%mod;
g[len-1]=0;
return;
}
ll tmp1_int[NNN],tmp2_int[NNN];
inline void Int(ll *f,ll *g,int len){//integral: g=\int f a.k.a.积分
#define inv tmp1_int
#define pr tmp2_int
pr[1]=1;
for(int i=2;i<=len;++i) pr[i]=pr[i-1]*i%mod;
inv[len]=qpow(pr[len],mod-2);
for(int i=len-1;i;--i) inv[i]=inv[i+1]*(i+1)%mod;
inv[1]=1;
for(int i=2;i<=len;++i) inv[i]=inv[i]*pr[i-1]%mod;
for(int i=1;i<len;++i) g[i]=f[i-1]*inv[i]%mod;
g[0]=0;
for(int i=0;i<=len;++i) inv[i]=pr[i]=0;
return;
#undef inv
#undef pr
}
ll tmp1_ln[NNN],tmp2_ln[NNN];
inline void Ln(ll *f,ll *g,int len){//napierian: g=ln f a.k.a.ln
#define p1 tmp1_ln//p1 for derivation
#define p2 tmp2_ln//p2 for inversion
Der(f,p1,len);
Inv(f,p2,len);
Con(p1,p2,len-1,len);
Int(p1,g,len);
for(int i=len;i<lim;++i) g[i]=0;
for(int i=0;i<lim;++i) p1[i]=p2[i]=0;
return;
#undef p1
#undef p2
}
ll tmp1_exp[NNN];
inline void Exp(ll *f,ll *g,int len){//exponential: g=e^f a.k.a.exp
#define h tmp1_exp
if(len==1){
g[0]=1;
return;
}
Exp(f,g,(len+1)>>1);
Ln(g,h,len);
for(int i=0;i<len;++i) h[i]=(f[i]-h[i]+mod)%mod;
for(int i=len;i<lim;++i) h[i]=0;
h[0]=(h[0]+1)%mod;
Con(g,h,len,len);
for(int i=0;i<lim;++i) h[i]=0;
for(int i=len;i<lim;++i) g[i]=0;
return;
#undef h
}
ll A[NNN],B[NNN];
int n;
int main(){
n=in;
for(int i=0;i<n;++i) A[i]=in;
Exp(A,B,n);
for(int i=0;i<n;++i) printf("%lld ",B[i]);
return 0;
}
多项式快速幂
弱化版
显然不能按普通数快速幂的方法算,这
k
k
k太大了
考虑推导
B
=
A
k
B
=
e
k
ln
A
B=A^k\\ B=e^{k\ln A}
B=AkB=eklnA
那就先ln,然后exp
显然ln的定义域不能有0,也就是弱化版和加强版的区别
由于k太大,利用某个模的性质 x k ≡ x k m o d p ( m o d p ) x^k\equiv x^{k\bmod p}\pmod p xk≡xkmodp(modp)
#include<bits/stdc++.h>
using namespace std;
#define in Read()
typedef long long ll;
inline ll in{
ll i=0,f=1;char ch=0;
while(ch!='-'&&!isdigit(ch)) ch=getchar();
if(ch=='-') ch=getchar(),f=-1;
while(isdigit(ch)) i=(i<<1)+(i<<3)+ch-48,ch=getchar();
return i*f;
}
const int NNN=5e5+5;
const ll mod=998244353;
int lim,rev[NNN];
inline ll get_k(){
ll i=0;char ch=0;
while(!isdigit(ch)) ch=getchar();
while(isdigit(ch)) i=(i*10%mod+ch-48)%mod,ch=getchar();
return i;
}
inline ll qpow(ll a,int x){
ll res=1;
while(x){
if(x&1) res=res*a%mod;
a=a*a%mod;
x>>=1;
}
return res;
}
inline void Init(int len){
lim=1;
while(lim<len) lim<<=1;
for(int i=0;i<lim;++i)
rev[i]=(rev[i>>1]>>1)|((i&1)?lim>>1:0);
return;
}
inline void NTT(ll *f,int sgn){
for(int i=0;i<lim;++i)
if(i<rev[i]) swap(f[i],f[rev[i]]);
for(int len=1;len<lim;len<<=1){
int siz=len<<1;
ll wn=qpow(3,(mod-1)/siz);
for(int l=0;l<lim;l+=siz){
ll w=1;
for(int i=l;i<l+len;++i){
ll a=f[i],b=f[i+len]*w%mod;
f[i]=(a+b)%mod;
f[i+len]=(a-b+mod)%mod;
w=w*wn%mod;
}
}
}
if(sgn) return;
ll inv=qpow(lim,mod-2);
reverse(f+1,f+lim);
for(int i=0;i<lim;++i)
f[i]=f[i]*inv%mod;
return;
}
inline void Con(ll *f,ll *g,int lf,int lg){
Init(lf+lg);
NTT(f,1);
NTT(g,1);
for(int i=0;i<lim;++i) f[i]=f[i]*g[i]%mod;
NTT(f,0);
return;
}
inline void Der(ll *f,ll *g,int len){
for(int i=0;i<len-1;++i) g[i]=f[i+1]*(i+1)%mod;
g[len-1]=0;
return;
}
ll tmp_inv[NNN];
inline void Inv(ll *f,ll *g,int len){
#define h tmp_inv
if(len==1){
g[0]=qpow(f[0],mod-2);
return;
}
Inv(f,g,(len+1)>>1);
Init(len<<1);
for(int i=0;i<len;++i) h[i]=f[i];
for(int i=len;i<lim;++i) h[i]=0;
NTT(g,1);
NTT(h,1);
for(int i=0;i<lim;++i) g[i]=1ll*(2ll-1ll*g[i]*h[i]%mod+mod)%mod*g[i]%mod;
NTT(g,0);
for(int i=len;i<lim;++i) g[i]=0;
return;
#undef h
}
ll tmp_int[NNN],pro_int[NNN];
inline void Int(ll *f,ll *g,int len){
#define inv tmp_int
#define pro pro_int
pro[1]=1;
for(int i=2;i<=len;++i) pro[i]=pro[i-1]*i%mod;
inv[len]=qpow(pro[len],mod-2);
for(int i=len-1;i;--i) inv[i]=inv[i+1]*(i+1)%mod;
inv[1]=1;
for(int i=2;i<=len;++i) inv[i]=inv[i]*pro[i-1]%mod;
for(int i=1;i<len;++i) g[i]=f[i-1]*inv[i]%mod;
g[0]=0;
for(int i=0;i<=len;++i) inv[i]=pro[i]=0;
return;
#undef inv
#undef pro
}
ll tmp1_ln[NNN],tmp2_ln[NNN];
inline void Ln(ll *f,ll *g,int len){
#define p1 tmp1_ln
#define p2 tmp2_ln
Der(f,p1,len);
Inv(f,p2,len);
Con(p1,p2,len-1,len);
Int(p1,g,len);
for(int i=len;i<lim;++i) g[i]=0;
for(int i=0;i<lim;++i) p1[i]=p2[i]=0;
return;
#undef p1
#undef p2
}
ll tmp_exp[NNN];
inline void Exp(ll *f,ll *g,int len){
#define h tmp_exp
if(len==1){
g[0]=1;
return;
}
Exp(f,g,(len+1)>>1);
Ln(g,h,len);
for(int i=0;i<len;++i) h[i]=(f[i]-h[i]+mod)%mod;
for(int i=len;i<lim;++i) h[i]=0;
h[0]=(h[0]+1)%mod;
Con(g,h,len,len);
for(int i=0;i<lim;++i) h[i]=0;
for(int i=len;i<lim;++i) g[i]=0;
return;
#undef h
}
ll tmp_pow[NNN];
inline void polypow(ll *f,ll *g,int len,int k){
#define tmp tmp_pow
Ln(f,tmp,len);
for(int i=0;i<len;++i) tmp[i]=tmp[i]*k%mod;
Exp(tmp,g,len);
return;
#undef tmp
}
int n;
ll k,A[NNN],B[NNN];
int main(){
n=in,k=get_k();
for(int i=0;i<n;++i) A[i]=in;
polypow(A,B,n,k);
for(int i=0;i<n;++i) printf("%lld ",B[i]);
return 0;
}
加强版
- 想办法避开
0
0
0,有
0
0
0的位置跳过,一直跳到没零为止
跳 0 0 0前最低位若为 x n x^n xn,那么快速幂之后最低位应为 x n k x^{nk} xnk,跳回去即可 - a 0 ≠ 1 a_0\neq1 a0=1时,你也不知道递归到 1 1 1时的次幂,那就除掉它
注意模了的
k
k
k不是原来的
k
k
k
综上所述,我们需要记录的内容有没有模过的
k
k
k和乘除
a
0
a_0
a0的幂次(注意费马小定理)
多项式开平方根
弱化版
推柿子(就懒得写同余了,反正都能理解)
B
=
A
B
2
−
A
=
0
,
记
B
2
−
A
=
g
(
B
)
显
然
牛
顿
迭
代
B
=
B
0
−
g
(
B
0
)
g
′
B
0
=
B
0
−
B
0
2
−
A
2
B
0
=
1
2
(
B
0
+
A
B
0
)
\begin{aligned} &B=\sqrt A\\ &B^2-A=0,记B^2-A=g(B)\\ &显然牛顿迭代\\ B&=B_0-\frac{g(B_0)}{g'B_0}\\ &=B_0-\frac{B_0^2-A}{2B_0}\\ &=\frac{1}{2}(B_0+\frac{A}{B_0}) \end{aligned}
BB=AB2−A=0,记B2−A=g(B)显然牛顿迭代=B0−g′B0g(B0)=B0−2B0B02−A=21(B0+B0A)
后面两个先咕了,去学习一下sszx码风/sszx秘传码风