前言
FFT
定义啊证明啊先咕咕咕了
放个版溜了
typedef complex<double> cpd;
void FFT(cpd p[],int n,int dir){
for(int i=1,j=0;i<n-1;i++){
for(int s=n;j^=s>>=1,~j&s;) ;
if(i<j){
swap(p[i],p[j]);
}
}
for(int d=0;(1<<d)<n;d++){
int m=1<<d,m2=m*2;
double p0=pi/m*dir;
cpd wn(cos(p0),sin(p0));
for(int i=0;i<n;i+=m2){
cpd w(1,0);
for(int j=0;j<m;j++){
cpd t0=p[i+j],t1=p[i+j+m];
p[i+j]=t0+t1*w;
p[i+j+m]=t0-t1*w;
w=w*wn;
}
}
}
if(dir==-1)
for(int i=0;i<n;i++)
p[i].real()/=1.0*n;
}
这个版注意一下,图省事直接用的复数类,所以效率比较慢,而且real()根据编译器问题可能CE?
NTT
其实啊…跟FFT差不了多少
void NTT(int *p,int n,int dir){
for(int i=1,j=0;i<n-1;i++){
for(int s=n;j^=s>>=1,~j&s;);
if(i<j)
swap(p[i],p[j]);
}
for(int m=1;m<n;m<<=1){
int m2=m*2;
int wn=ksm(g,(mod-1)/m2);
if(dir==-1)
wn=ksm(wn,mod-2);
for(int i=0;i<n;i+=m2){
int w=1;
for(int j=0;j<m;j++){
int t0=p[i+j],t1=p[i+j+m];
p[i+j]=(t0+1ll*t1*w%mod)%mod;
p[i+j+m]=(t0-1ll*t1*w%mod+mod)%mod;
w=1ll*w*wn%mod;
}
}
}
if(dir==-1){
int inv=ksm(n,mod-2);
for(int i=0;i<n;i++)
p[i]=1ll*p[i]*inv%mod;
}
}
多项式的乘法
em…就是卷积嘛…
上FFT或者NTT即可
注意,一定注意清零的问题,下面的算法也是一样
//注意执行完这个函数后a、b会改变哦
void Mul(Poly &a,Poly &b,int n,Poly &res){
int len=1;
while(len<2*n)
len<<=1;
NTT(a,len,1);
NTT(b,len,1);
for(int i=0;i<len;i++)
res[i]=1ll*a[i]*b[i]%mod;
NTT(res,len,-1);
}
多项式求逆元
A
(
x
)
B
(
x
)
≡
1
(
m
o
d
x
n
)
A(x)B(x)\equiv1(mod\ x^n)
A(x)B(x)≡1(mod xn)
当n=1时,就剩常数项了,直接逆元
当n>1时,假设我们已经求出了关于
x
⌈
n
2
⌉
x^{\lceil \frac{n}{2}\rceil}
x⌈2n⌉ 的逆元,也就是说
A
(
x
)
B
0
(
x
)
≡
1
(
m
o
d
x
⌈
n
2
⌉
)
⋅
⋅
⋅
①
A(x)B_0(x)\equiv1(mod\ x^{\lceil \frac{n}{2}\rceil})···①
A(x)B0(x)≡1(mod x⌈2n⌉)⋅⋅⋅①
用原式-①,有
B
(
x
)
−
B
0
(
x
)
≡
0
(
m
o
d
x
⌈
n
2
⌉
)
B(x)-B_0(x)\equiv0(mod\ x^{\lceil \frac{n}{2}\rceil})
B(x)−B0(x)≡0(mod x⌈2n⌉)
两边平方,有
(
B
(
x
)
−
B
0
(
x
)
)
2
≡
0
(
m
o
d
x
⌈
n
2
⌉
)
(B(x)-B_0(x))^{2} \equiv 0(mod\ x^{\lceil \frac{n}{2}\rceil})
(B(x)−B0(x))2≡0(mod x⌈2n⌉)
由于,左边那个东西在0到(n-1)/2上都是0,所以平方之后,0到n-1上一定都是0
两边同时乘上
A
(
x
)
A(x)
A(x),移项,有
B
(
x
)
≡
2
B
0
(
x
)
−
A
(
x
)
B
0
2
(
x
)
(
m
o
d
x
n
)
B(x)\equiv2B_0(x)-A(x)B_0^2(x)(mod \ x^n)
B(x)≡2B0(x)−A(x)B02(x)(mod xn)
右边那一坨每次用NTT去搞就好了
void Inv(const Poly &a,Poly &inv,int n){
if(n==1){
inv[0]=ksm(a[0],mod-2);
return ;
}
Inv(a,inv,(n+1)/2);
int len=1;
while(len<2*n)
len<<=1;
copy(a,a+n,t);
fill(t+n,t+len,0);
fill(inv+n,inv+len,0);
NTT(t,len,1);
NTT(inv,len,1);
for(int i=0;i<len;i++)
inv[i]=1ll*inv[i]*((2-1ll*t[i]*inv[i]%mod+mod)%mod)%mod;
NTT(inv,len,-1);
fill(inv+n,inv+len,0);
}
多项式取模
A
(
x
)
=
P
(
x
)
Q
(
x
)
+
R
(
x
)
A(x)=P(x)Q(x)+R(x)
A(x)=P(x)Q(x)+R(x)
设
d
e
g
(
A
(
x
)
)
=
n
,
d
e
g
(
P
(
x
)
)
=
m
deg(A(x))=n,deg(P(x))=m
deg(A(x))=n,deg(P(x))=m
那么
d
e
g
(
Q
(
x
)
)
=
n
−
m
,
d
e
g
(
R
(
x
)
)
<
=
m
−
1
deg(Q(x))=n-m,deg(R(x))<=m-1
deg(Q(x))=n−m,deg(R(x))<=m−1
A
(
1
x
)
=
P
(
1
x
)
Q
(
1
x
)
+
R
(
1
x
)
A(\frac{1}{x})=P(\frac{1}{x})Q(\frac{1}{x})+R(\frac{1}{x})
A(x1)=P(x1)Q(x1)+R(x1)
x
n
A
(
1
x
)
=
x
n
P
(
1
x
)
Q
(
1
x
)
+
x
n
R
(
1
x
)
x^nA(\frac{1}{x})=x^nP(\frac{1}{x})Q(\frac{1}{x})+x^nR(\frac{1}{x})
xnA(x1)=xnP(x1)Q(x1)+xnR(x1)
x
n
A
(
1
x
)
=
x
m
P
(
1
x
)
x
n
−
m
Q
(
1
x
)
+
x
n
−
m
+
1
x
m
−
1
R
(
1
x
)
x^nA(\frac{1}{x})=x^mP(\frac{1}{x})x^{n-m}Q(\frac{1}{x})+x^{n-m+1}x^{m-1}R(\frac{1}{x})
xnA(x1)=xmP(x1)xn−mQ(x1)+xn−m+1xm−1R(x1)
观察可以发现
x
n
A
(
1
x
)
其
实
就
相
当
于
A
(
x
)
把
所
有
系
数
翻
转
了
,
我
们
记
翻
转
后
的
多
项
式
为
A
r
(
x
)
x^nA(\frac{1}{x})其实就相当于A(x)把所有系数翻转了,我们记翻转后的多项式为A^r(x)
xnA(x1)其实就相当于A(x)把所有系数翻转了,我们记翻转后的多项式为Ar(x)
A
r
(
x
)
=
P
r
(
x
)
Q
r
(
x
)
+
x
n
−
m
+
1
R
r
(
x
)
A^r(x)=P^r(x)Q^r(x)+x^{n-m+1}R^r(x)
Ar(x)=Pr(x)Qr(x)+xn−m+1Rr(x)
顺手取个模,
A
r
(
x
)
≡
P
r
(
x
)
Q
r
(
x
)
(
m
o
d
x
n
−
m
+
1
)
A^r(x)\equiv P^r(x)Q^r(x)(mod \ x^{n-m+1})
Ar(x)≡Pr(x)Qr(x)(mod xn−m+1)
Q
r
(
x
)
≡
A
r
(
x
)
P
r
(
x
)
(
m
o
d
x
n
−
m
+
1
)
Q^r(x)\equiv \frac{A^r(x)}{P^{r}(x)}(mod \ x^{n-m+1})
Qr(x)≡Pr(x)Ar(x)(mod xn−m+1)
右边就逆元搞定咯
然后就可以求出
Q
r
(
x
)
Q^r(x)
Qr(x)
然后再翻转,就得到了
Q
(
x
)
Q(x)
Q(x)
那么
R
(
x
)
=
A
(
x
)
−
P
(
x
)
Q
(
x
)
R(x)=A(x)-P(x)Q(x)
R(x)=A(x)−P(x)Q(x)搞定。
void Mod(Poly &a,const Poly &p,int n,int m){
copy(a,a+n,ra);
reverse(ra,ra+n);
copy(p,p+m,rp);
reverse(rp,rp+m);
Inv(rp,rpi,n-m+1);
int len=1;
while(len<(n-m+1)*2)
len<<=1;
fill(ra+n-m+1,ra+len,0);
fill(rpi+n-m+1,rpi+len,0);
Mul(ra,rpi,n-m+1,c);
len=1;
while(len<2*n)
len<<=1;
fill(c+n-m+1,c+len,0);
reverse(c,c+n-m+1);
copy(p,p+m,t);
fill(t+m,t+len,0);
Mul(c,t,n,c);
for(int i=0;i<n;i++)
a[i]=(a[i]-c[i]+mod)%mod;
}
多项式的多点求值
问题描述:给出一个多项式A(x),和n个横坐标xi,求出每个对应的A(xi)
我们令所有要求值的横坐标集合为X
然后定义
X
0
=
{
x
0
,
x
1
,
⋯
x
⌊
n
2
⌋
}
X_0=\{x_0,x_1,\cdots x_{\lfloor\frac{n}{2}\rfloor}\}
X0={x0,x1,⋯x⌊2n⌋}
X
1
=
{
x
⌊
n
2
⌋
+
1
,
x
⌊
n
2
⌋
+
2
,
⋯
x
n
−
1
}
X_1=\{x_{\lfloor\frac{n}{2}\rfloor+1},x_{\lfloor\frac{n}{2}\rfloor+2},\cdots x_{n-1}\}
X1={x⌊2n⌋+1,x⌊2n⌋+2,⋯xn−1}
记
X
0
和
X
1
插
值
得
到
的
多
项
式
分
别
为
A
0
,
A
1
记X_0和X_1插值得到的多项式分别为A_0,A_1
记X0和X1插值得到的多项式分别为A0,A1
然后我们定义两个多项式
P
0
(
x
)
=
∏
i
=
0
⌊
n
2
⌋
(
x
−
x
i
)
P_0(x)=\prod_{i=0}^{\lfloor\frac{n}{2}\rfloor}(x-x_i)
P0(x)=i=0∏⌊2n⌋(x−xi)
P
1
(
x
)
=
∏
i
=
⌊
n
2
⌋
+
1
n
−
1
(
x
−
x
i
)
P_1(x)=\prod_{i=\lfloor\frac{n}{2}\rfloor+1}^{n-1}(x-x_i)
P1(x)=i=⌊2n⌋+1∏n−1(x−xi)
那么可以有
A
(
x
)
=
Q
(
x
)
P
0
(
x
)
+
A
0
(
x
)
A(x)=Q(x)P_0(x)+A_0(x)
A(x)=Q(x)P0(x)+A0(x)
也即是说
A
0
(
x
)
≡
A
(
x
)
(
m
o
d
P
0
(
x
)
)
A_0(x)\equiv A(x)(mod \ P_0(x))
A0(x)≡A(x)(mod P0(x))
同理有
A
1
(
x
)
≡
A
(
x
)
(
m
o
d
P
1
(
x
)
)
A_1(x)\equiv A(x)(mod \ P_1(x))
A1(x)≡A(x)(mod P1(x))
然后往下带
然后到n=1的时候,就可以直接算出答案了
时间复杂度
O
(
n
l
o
g
2
n
)
O(nlog^2n)
O(nlog2n)
以上。
模板题
#include<cstdio>
#include<cstring>
#include<algorithm>
#define lch x<<1
#define rch x<<1|1
using namespace std;
typedef long long ll;
const int N=64005*8;
const int mod=998244353;
const int g=3;
int n,m;
void Debug(const int *p,char c,int n=8){
putchar(c);
for(int i=0;i<n;i++)
printf("%d ",p[i]);
puts("");
}
int ksm(int x,int y){
int res=1;
while(y){
if(y&1)
res=1ll*res*x%mod;
x=1ll*x*x%mod;
y>>=1;
}
return res;
}
void NTT(int *p,int n,int dir){
for(int i=1,j=0;i<n-1;i++){
for(int s=n;j^=s>>=1,~j&s;);
if(i<j)
swap(p[i],p[j]);
}
for(int m=1;m<n;m<<=1){
int m2=m*2;
int wn=ksm(g,(mod-1)/m2);
if(dir==-1)
wn=ksm(wn,mod-2);
for(int i=0;i<n;i+=m2){
int w=1;
for(int j=0;j<m;j++){
int t0=p[i+j],t1=p[i+j+m];
p[i+j]=(t0+1ll*t1*w%mod)%mod;
p[i+j+m]=(t0-1ll*t1*w%mod+mod)%mod;
w=1ll*w*wn%mod;
}
}
}
if(dir==-1){
int inv=ksm(n,mod-2);
for(int i=0;i<n;i++)
p[i]=1ll*p[i]*inv%mod;
}
}
void Mul(int *a,int *b,int n,int *res){
int len=1;
while(len<2*n)
len<<=1;
fill(a+n,a+len,0);
fill(b+n,b+len,0);
NTT(a,len,1);
NTT(b,len,1);
for(int i=0;i<len;i++)
res[i]=1ll*a[i]*b[i]%mod;
NTT(res,len,-1);
//fill(res+n,res+len,0);
}
void Inv(const int *a,int *inv,int n){
static int t[N];
if(n==1){
inv[0]=ksm(a[0],mod-2);
return ;
}
Inv(a,inv,(n+1)/2);
int len=1;
while(len<2*n)
len<<=1;
copy(a,a+n,t);
fill(t+n,t+len,0);
fill(inv+n,inv+len,0);
NTT(t,len,1);
NTT(inv,len,1);
for(int i=0;i<len;i++)
inv[i]=1ll*inv[i]*((2-1ll*inv[i]*t[i]%mod+mod)%mod)%mod;
NTT(inv,len,-1);
fill(inv+n,inv+len,0);
}
void Mod(int *a,const int *p,int n,int m){
static int ra[N],rp[N],rpi[N],c[N],d[N];
copy(a,a+n,ra);
reverse(ra,ra+n);
copy(p,p+m,rp);
reverse(rp,rp+m);
Inv(rp,rpi,n-m+1);
int len=1;
while(len<(n-m+1)*2)
len<<=1;
fill(ra+n-m+1,ra+len,0);
fill(rpi+n-m+1,rpi+len,0);
Mul(ra,rpi,n-m+1,c);
len=1;
while(len<2*n)
len<<=1;
fill(c+n-m+1,c+len,0);
reverse(c,c+n-m+1);
copy(p,p+m,d);
fill(d+m,d+len,0);
Mul(c,d,n,c);
for(int i=0;i<n;i++)
a[i]=(a[i]-c[i]+mod)%mod;
}
int *P[N*4],len[N*4];
int a[N],b[N];
void Init_P(int x,int l,int r){
static int t1[N],t2[N];
if(l==r){
len[x]=1;
P[x]=new int[2];
P[x][0]=(mod-b[l]%mod);
P[x][1]=1;
return ;
}
int mid=(l+r)>>1;
Init_P(lch,l,mid);
Init_P(rch,mid+1,r);
len[x]=(len[lch]+len[rch]);
//Debug(P[lch],'l',len[lch]+1);
//Debug(P[rch],'r',len[lch]+1);
P[x]=new int [len[x]+1];
int n=len[x]+1,Len=1;
while(Len<n*2)
Len<<=1;
copy(P[lch],P[lch]+len[lch]+1,t1);
fill(t1+len[lch]+1,t1+Len,0);
copy(P[rch],P[rch]+len[rch]+1,t2);
fill(t2+len[rch]+1,t2+Len,0);
//Debug(t1,'@');
//Debug(t2,'@');
Mul(t1,t2,n,t1);
//Debug(t1,'!');
for(int i=0;i<n;i++)
P[x][i]=t1[i];
}
int ans[N];
void calc(int x,int l,int r,const int *A){
if(l==r){
ans[l]=A[0];
return ;
}
int mid=(l+r)>>1;
int B[(len[x]+2)<<2];
int n=len[x]+1;
int Len=1;
while(Len<n*2)
Len<<=1;
copy(A,A+n,B);
fill(B+n,B+Len,0);
//Debug(B,'$');
//Debug(P[lch],'$',len[lch]+1);
Mod(B,P[lch],n,len[lch]+1);
//Debug(B,'#');
calc(lch,l,mid,B);
copy(A,A+n,B);
fill(B+n,B+Len,0);
Mod(B,P[rch],n,len[rch]+1);
calc(rch,mid+1,r,B);
}
void Read(int &x){
x=0;
char c=getchar();
while(c<'0'||c>'9')
c=getchar();
while(c>='0'&&c<='9')
x=x*10+c-'0',c=getchar();
}
int main()
{
//freopen("1.in","r",stdin);
//scanf("%d%d",&n,&m);
Read(n),Read(m);
for(int i=0;i<=n;i++){
//scanf("%d",&a[i]);
Read(a[i]);
}
for(int i=1;i<=m;i++){
//scanf("%d",&b[i]);
Read(b[i]);
}
Init_P(1,1,m);
/*if(n>=m)
Mod(a,P[1],n+1,m+1);*/
calc(1,1,m,a);
for(int i=1;i<=m;i++)
printf("%d\n",ans[i]);
}
拉格朗日插值法
问题描述:给出一个n-1次多项式上的n个点,求
f
(
x
)
f(x)
f(x)
然而似乎不能求出来
f
(
x
)
f(x)
f(x)的具体表达式,只能对于一个特定的
x
0
x_0
x0求出对应的
f
(
x
0
)
f(x_0)
f(x0)
直接上公式了,证明直接咕掉
f
(
x
)
=
∑
i
=
1
n
y
i
∏
j
=
1
,
j
≠
i
n
x
−
x
j
x
i
−
x
j
f(x)=\sum_{i=1}^ny_i\prod_{j=1,j\neq i}^{n}\frac{x-x_j}{x_i-x_j}
f(x)=i=1∑nyij=1,j̸=i∏nxi−xjx−xj
如果预处理逆元的话,就
O
(
n
2
)
O(n^2)
O(n2)咯
多项式的快速插值
法一
这个方法我没写哦,先咕了,就当是扩展思路?
还是按照上面类似的方法
把
X
分
成
X
0
和
Y
0
把X分成X_0和Y_0
把X分成X0和Y0
假如我们已经求出了
X
0
X_0
X0所对应的插值多项式
A
0
(
x
)
A_0(x)
A0(x)
设
P
(
x
)
=
∏
i
=
0
⌊
n
2
⌋
(
x
−
x
i
)
P(x)=\prod_{i=0}^{\lfloor\frac{n}{2}\rfloor}(x-x_i)
P(x)=i=0∏⌊2n⌋(x−xi)
那么再设
A
(
x
)
=
Q
(
x
)
P
(
x
)
+
A
0
(
x
)
A(x)=Q(x)P(x)+A_0(x)
A(x)=Q(x)P(x)+A0(x)
当
i
≤
⌊
n
2
⌋
i\leq\lfloor\frac{n}{2}\rfloor
i≤⌊2n⌋时,显然成立
那么当
i
>
⌊
n
2
⌋
i>\lfloor\frac{n}{2}\rfloor
i>⌊2n⌋时,有
y
i
=
Q
(
x
i
)
P
(
x
i
)
+
A
0
(
x
i
)
y_i=Q(x_i)P(x_i)+A_0(x_i)
yi=Q(xi)P(xi)+A0(xi)
也就是说
Q
(
x
i
)
=
y
i
−
A
0
(
x
i
)
P
(
x
i
)
Q(x_i)=\frac{y_i-A_0(x_i)}{P(x_i)}
Q(xi)=P(xi)yi−A0(xi)
后面那坨可以看做是一个新的点
(
x
,
y
i
−
A
0
(
x
i
)
P
(
x
i
)
)
(x,\frac{y_i-A_0(x_i)}{P(x_i)})
(x,P(xi)yi−A0(xi))在
Q
(
x
i
)
Q(x_i)
Q(xi)上,那么就也可以用快速插值解决
至于
P
(
x
i
)
和
A
0
(
x
i
)
P(x_i)和A_0(x_i)
P(xi)和A0(xi)的值,可以用多点求值解决。
总时间复杂度
O
(
n
l
o
g
3
n
)
O(nlog^3n)
O(nlog3n)
还有
O
(
n
l
o
g
2
n
)
O(nlog^2n)
O(nlog2n) 的算法呢
法二
首先,由拉格朗日插值法有
f
(
x
)
=
∑
i
=
1
n
y
i
∏
j
=
1
,
j
≠
i
n
x
−
x
j
x
i
−
x
j
f(x)=\sum_{i=1}^ny_i\prod_{j=1,j\neq i}^{n}\frac{x-x_j}{x_i-x_j}
f(x)=i=1∑nyij=1,j̸=i∏nxi−xjx−xj
变形,有
f
(
x
)
=
∑
i
=
1
n
y
i
∏
j
=
1
,
j
≠
i
n
(
x
i
−
x
j
)
∏
j
=
1
,
j
≠
=
i
n
(
x
−
x
j
)
f(x)=\sum_{i=1}^n\frac{y_i}{\prod_{j=1,j\neq i}^n(x_i-x_j)}\prod_{j=1,j\neq =i}^n(x-x_j)
f(x)=i=1∑n∏j=1,j̸=in(xi−xj)yij=1,j̸==i∏n(x−xj)
我们令
G
(
x
)
=
∏
j
=
1
n
(
x
−
x
j
)
G(x)=\prod_{j=1}^n(x-x_j)
G(x)=j=1∏n(x−xj)
那么,
∏
j
=
0
,
j
≠
i
n
(
x
−
x
j
)
=
G
(
x
)
x
−
x
i
\prod_{j=0,j\neq i}^n(x-x_j)=\frac{G(x)}{x-x_i}
j=0,j̸=i∏n(x−xj)=x−xiG(x)
记左式为H(x)
H
(
x
)
(
x
−
x
i
)
=
G
(
x
)
H(x)(x-x_i)=G(x)
H(x)(x−xi)=G(x)
求导
H
(
x
)
+
(
x
−
x
i
)
H
′
(
x
)
=
G
′
(
x
)
H(x)+(x-x_i)H'(x)=G'(x)
H(x)+(x−xi)H′(x)=G′(x)
令x=xi,有
H
(
x
i
)
=
G
′
(
x
i
)
H(x_i)=G'(x_i)
H(xi)=G′(xi)
那么,
F
(
x
)
=
∑
i
=
0
n
y
i
G
′
(
x
i
)
G
(
x
)
x
−
x
i
F(x)=\sum_{i=0}^n\frac{y_i}{G'(x_i)}\frac{G(x)}{x-x_i}
F(x)=i=0∑nG′(xi)yix−xiG(x)
然后呢,我们令
z
(
i
)
=
y
i
G
′
(
x
i
)
z(i)=\frac{y_i}{G'(x_i)}
z(i)=G′(xi)yi,
再令
L
(
x
)
=
∏
i
=
1
⌊
n
2
⌋
(
x
−
x
i
)
L(x)=\prod_{i=1}^{\lfloor\frac{n}{2}\rfloor}(x-x_i)
L(x)=i=1∏⌊2n⌋(x−xi)
R
(
x
)
=
∏
i
=
⌊
n
2
⌋
+
1
n
(
x
−
x
i
)
R(x)=\prod_{i=\lfloor\frac{n}{2}\rfloor+1}^{n}(x-x_i)
R(x)=i=⌊2n⌋+1∏n(x−xi)
P
0
(
x
)
=
∑
i
=
1
⌊
n
2
⌋
z
i
∏
j
≠
i
⌊
n
2
⌋
(
x
−
x
j
)
P_0(x)=\sum_{i=1}^{\lfloor\frac{n}{2}\rfloor}z_i\prod_{j\neq i}^{\lfloor\frac{n}{2}\rfloor}(x-x_j)
P0(x)=i=1∑⌊2n⌋zij̸=i∏⌊2n⌋(x−xj)
P
1
(
x
)
=
∑
i
=
⌊
n
2
⌋
+
1
n
z
i
∏
⌊
n
2
⌋
+
1
,
j
≠
i
n
(
x
−
x
j
)
P_1(x)=\sum_{i=\lfloor\frac{n}{2}\rfloor+1}^{n}z_i\prod_{\lfloor\frac{n}{2}\rfloor+1,j\neq i}^{n}(x-x_j)
P1(x)=i=⌊2n⌋+1∑nzi⌊2n⌋+1,j̸=i∏n(x−xj)
那么,
F
(
x
)
=
R
(
x
)
P
0
(
x
)
+
L
(
x
)
P
1
(
x
)
F(x)=R(x)P_0(x)+L(x)P_1(x)
F(x)=R(x)P0(x)+L(x)P1(x)
L
(
x
)
,
R
(
x
)
L(x),R(x)
L(x),R(x)可以预先NTT分治得到
z
(
i
)
z(i)
z(i)可以用多点求值得到
然后就再做一次分治即可。
总时间复杂度
O
(
n
l
o
g
2
n
)
O(nlog^2n)
O(nlog2n)
#include<cstdio>
#include<cstring>
#include<algorithm>
#define lch x<<1
#define rch x<<1|1
using namespace std;
typedef long long ll;
const int N=100005*4;
const int mod=998244353;
const int g=3;
int n,m;
void Debug(const int *p,char c,int n=8){
putchar(c);
for(int i=0;i<n;i++)
printf("%d ",p[i]);
puts("");
}
int ksm(int x,int y){
int res=1;
while(y){
if(y&1)
res=1ll*res*x%mod;
x=1ll*x*x%mod;
y>>=1;
}
return res;
}
int t_wn[N],t_wni[N];
void Prepare(){
for(int m=1;m<N;m<<=1){
int m2=m*2;
t_wn[m]=ksm(g,(mod-1)/m2);
t_wni[m]=ksm(t_wn[m],mod-2);
}
}
void NTT(int *p,int n,int dir){
for(int i=1,j=0;i<n-1;i++){
for(int s=n;j^=s>>=1,~j&s;);
if(i<j)
swap(p[i],p[j]);
}
for(int m=1;m<n;m<<=1){
int m2=m*2;
/*int wn=ksm(g,(mod-1)/m2);
if(dir==-1)
wn=ksm(wn,mod-2);*/
int wn;
if(dir==1){
wn=t_wn[m];
}
else{
wn=t_wni[m];
}
for(int i=0;i<n;i+=m2){
int w=1;
for(int j=0;j<m;j++){
int t0=p[i+j],t1=p[i+j+m];
p[i+j]=(t0+1ll*t1*w%mod)%mod;
p[i+j+m]=(t0-1ll*t1*w%mod+mod)%mod;
w=1ll*w*wn%mod;
}
}
}
if(dir==-1){
int inv=ksm(n,mod-2);
for(int i=0;i<n;i++)
p[i]=1ll*p[i]*inv%mod;
}
}
void Mul(int *a,int *b,int len,int *res){
NTT(a,len,1);
NTT(b,len,1);
for(int i=0;i<len;i++)
res[i]=1ll*a[i]*b[i]%mod;
NTT(res,len,-1);
//fill(res+n,res+len,0);
}
void Inv(const int *a,int *inv,int n){
static int t[N];
if(n==1){
inv[0]=ksm(a[0],mod-2);
return ;
}
Inv(a,inv,(n+1)/2);
int len=1;
while(len<n*2)
len<<=1;
copy(a,a+n,t);
fill(t+n,t+len,0);
fill(inv+n,inv+len,0);
NTT(t,len,1);
NTT(inv,len,1);
for(int i=0;i<len;i++)
inv[i]=1ll*inv[i]*((2-1ll*inv[i]*t[i]%mod+mod)%mod)%mod;
NTT(inv,len,-1);
fill(inv+n,inv+len,0);
}
void Mod(int *a,const int *p,int n,int m){
static int ra[N],rp[N],rpi[N],c[N],d[N];
copy(a,a+n,ra);
reverse(ra,ra+n);
copy(p,p+m,rp);
reverse(rp,rp+m);
Inv(rp,rpi,n-m+1);
int len=1;
while(len<(n-m+1)*2)
len<<=1;
fill(ra+n-m+1,ra+len,0);
fill(rpi+n-m+1,rpi+len,0);
Mul(ra,rpi,len,c);
len=1;
while(len<n)
len<<=1;
fill(c+n-m+1,c+len,0);
reverse(c,c+n-m+1);
copy(p,p+m,d);
fill(d+m,d+len,0);
Mul(c,d,len,c);
for(int i=0;i<n;i++)
a[i]=(a[i]-c[i]+mod)%mod;
}
int *P[N*4],len[N*4];
int z[N];
int X[N],Y[N];
void Init_P(int x,int l,int r){
static int t1[N],t2[N];
if(l==r){
len[x]=1;
P[x]=new int[2];
P[x][0]=(mod-X[l]%mod);
P[x][1]=1;
return ;
}
int mid=(l+r)>>1;
Init_P(lch,l,mid);
Init_P(rch,mid+1,r);
len[x]=(len[lch]+len[rch]);
//Debug(P[lch],'l',len[lch]+1);
//Debug(P[rch],'r',len[lch]+1);
P[x]=new int [len[x]+1];
int n=len[x]+1,Len=1;
while(Len<n)
Len<<=1;
int lsz=len[lch]+1,rsz=len[rch]+1;
copy(P[lch],P[lch]+lsz,t1);
fill(t1+lsz,t1+Len,0);
copy(P[rch],P[rch]+rsz,t2);
fill(t2+rsz,t2+Len,0);
//Debug(t1,'@');
//Debug(t2,'@');
Mul(t1,t2,Len,t1);
//Debug(t1,'!');
for(int i=0;i<n;i++)
P[x][i]=t1[i];
}
void calc(int x,int l,int r,const int *A){
if(l==r){
z[l]=(1ll*Y[l]*ksm(A[0],mod-2)%mod+mod)%mod;
return ;
}
int mid=(l+r)>>1;
int n=len[x]+1;
int Len=1;
while(Len<n)
Len<<=1;
int B[Len];
copy(A,A+n,B);
fill(B+n,B+Len,0);
//Debug(B,'$');
//Debug(P[lch],'$',len[lch]+1);
Mod(B,P[lch],n,len[lch]+1);
//Debug(B,'#');
calc(lch,l,mid,B);
copy(A,A+n,B);
fill(B+n,B+Len,0);
Mod(B,P[rch],n,len[rch]+1);
calc(rch,mid+1,r,B);
}
void qiudao(int *x,int n){
for(int i=0;i<n;i++){
if(i){
x[i-1]=1ll*i*x[i]%mod;
}
x[i]=0;
}
}
int G[N];
void Solve(int x,int l,int r,int *F){
static int p0[N],p1[N];
if(l==r){
F[0]=1ll*z[l]%mod;
return ;
}
int lsz=len[lch]+1;
int rsz=len[rch]+1;
int Len=1;
int n=len[x]+1;
while(Len<n)
Len<<=1;
int L[Len];
int R[Len];
fill(L,L+Len,0);
fill(R,R+Len,0);
int mid=(l+r)>>1;
Solve(lch,l,mid,L);
Solve(rch,mid+1,r,R);
//puts("!");
//Debug(L,'L',(len[lch]+1)<<1);
//Debug(R,'R',(len[rch]+1)<<1);
copy(P[lch],P[lch]+lsz,p0);
fill(p0+lsz,p0+Len,0);
copy(P[rch],P[rch]+rsz,p1);
fill(p1+rsz,p1+Len,0);
//Debug(p0,'(',Len);
//Debug(p1,')',Len);
fill(L+lsz,L+Len,0);
fill(R+rsz,R+Len,0);
Mul(p0,R,Len,p0);
Mul(p1,L,Len,p1);
//Debug(p0,'(',Len);
//Debug(p1,')',Len);
for(int i=0;i<n;i++)
F[i]=(p0[i]+p1[i])%mod;
}
int ans[N];
void Read(int &x){
x=0;
char c=getchar();
while(c<'0'||c>'9')
c=getchar();
while(c>='0'&&c<='9')
x=x*10+c-'0',c=getchar();
}
void Print(int x){
if(x>9)
Print(x/10);
putchar(x%10+'0');
}
int main()
{
//freopen("1.in","r",stdin);
//ReadOptmized();
Prepare();
Read(n);
for(int i=1;i<=n;i++){
Read(X[i]);
Read(Y[i]);
}
Init_P(1,1,n);
copy(P[1],P[1]+n+1,G);
//Debug(G,'!');
qiudao(G,n+1);
//Debug(G,'@');
calc(1,1,n,G);
Solve(1,1,n,ans);
//puts("?");
for(int i=0;i<n;i++){
//printf("%d ",ans[i]);
Print(ans[i]);
putchar(' ');
}
//return 0;
//printf("success\n");
}
多项式的开方
已知A(x),
B
2
(
x
)
≡
A
(
x
)
(
m
o
d
x
n
)
B^2(x)\equiv A(x) (mod\ x^n)
B2(x)≡A(x)(mod xn)求B(x)
假设有
B
0
2
(
x
)
≡
A
(
x
)
(
m
o
d
x
⌊
n
2
⌋
)
B_0^2(x)\equiv A(x) (mod\ x^{\lfloor \frac{n}{2}\rfloor})
B02(x)≡A(x)(mod x⌊2n⌋)
两式相减再平方,有
B
2
(
x
)
−
2
B
(
x
)
B
0
(
x
)
+
B
0
2
(
x
)
≡
0
(
m
o
d
x
n
)
B^2(x)-2B(x)B_0(x)+B^2_0(x)\equiv0(mod\ x^n)
B2(x)−2B(x)B0(x)+B02(x)≡0(mod xn)
此处模数平方的原因是跟求逆元时的解释一样的
那么最后就有
B
(
x
)
≡
A
(
x
)
+
B
0
2
(
x
)
2
B
0
(
x
)
(
m
o
d
x
n
)
B(x)\equiv\frac{A(x)+B_0^2(x)}{2B_0(x)}(mod\ x^n)
B(x)≡2B0(x)A(x)+B02(x)(mod xn)
void Sqrt(const int *a,int *b,int n){
static int bi[N],t[N];
if(n==1){
b[0]=sqrt(a[0]);
return ;
}
Sqrt(a,b,(n+1)>>1);
Inv(b,bi,n);
copy(a,a+n,t);
int len=1;
while(len<n*2)
len<<=1;
fill(t+n,t+len,0);
fill(b+n,b+len,0);
fill(bi+n,bi+len,0);
NTT(t,len,1);
NTT(b,len,1);
NTT(bi,len,1);
for(int i=0;i<len;i++)
b[i]=1ll*(1ll*b[i]*b[i]+t[i]%mod)%mod*bi[i]%mod*inv2%mod;
NTT(b,len,-1);
fill(b+n,b+len,0);
}
牛顿迭代
推荐先看一下这篇文章(虽然似乎跟我们即将讲的关联不大…)
已知G(x),
G
(
F
(
x
)
)
≡
0
(
m
o
d
x
n
)
G(F(x))\equiv 0 (mod\ x^n)
G(F(x))≡0(mod xn),求F(x)
假设已经求出F_0(x)满足
G
(
F
0
(
x
)
)
≡
0
(
m
o
d
x
⌊
n
2
⌋
)
G(F_0(x))\equiv 0 (mod\ x^{\lfloor \frac{n}{2}\rfloor})
G(F0(x))≡0(mod x⌊2n⌋)
那么就有
F
(
x
)
≡
F
0
(
x
)
−
G
(
F
0
(
x
)
)
G
′
(
F
0
(
x
)
)
(
m
o
d
x
n
)
F(x)\equiv F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))}(mod\ x^n)
F(x)≡F0(x)−G′(F0(x))G(F0(x))(mod xn)
证明先咕了
多项式的ln
B
(
x
)
=
l
n
(
A
(
x
)
)
B(x)=ln(A(x))
B(x)=ln(A(x))
同时求导
B
′
(
x
)
=
A
′
(
x
)
A
(
x
)
B'(x)=\frac{A'(x)}{A(x)}
B′(x)=A(x)A′(x)
所以就酱紫了
(常数项?不存在的)
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const int N=(1e5+5)*4;
const int g=3;
const int mod=998244353;
void Debug(const int *p,char c,int n=8){
putchar(c);
for(int i=0;i<n;i++)
printf("%d ",p[i]);
puts("");
}
int ksm(int x,int y){
int res=1;
while(y){
if(y&1){
res=1ll*res*x%mod;
}
x=1ll*x*x%mod;
y>>=1;
}
return res;
}
void NTT(int *a,int n,int dir){
for(int i=1,j=0;i<n-1;i++){
for(int s=n;j^=s>>=1,~j&s;);
if(i<j)
swap(a[i],a[j]);
}
for(int m=1;m<n;m<<=1){
int m2=m<<1;
int w0=ksm(g,(mod-1)/m2);
if(dir==-1)
w0=ksm(w0,mod-2);
for(int i=0;i<n;i+=m2){
int w=1;
for(int j=0;j<m;j++){
int t0=a[i+j],t1=a[i+j+m];
a[i+j]=(t0+1ll*t1*w%mod)%mod;
a[i+j+m]=(t0-1ll*t1*w%mod+mod)%mod;
w=1ll*w*w0%mod;
}
}
}
if(dir==-1){
int inv=ksm(n,mod-2);
for(int i=0;i<n;i++)
a[i]=1ll*a[i]*inv%mod;
}
}
void Inv(const int *a,int *inv,int n){
static int t[N];
if(n==1){
inv[0]=ksm(a[0],mod-2);
return ;
}
Inv(a,inv,(n+1)/2);
int len=1;
while(len<2*n)
len<<=1;
copy(a,a+n,t);
fill(t+n,t+len,0);
fill(inv+n,inv+len,0);
NTT(t,len,1);
NTT(inv,len,1);
for(int i=0;i<len;i++)
inv[i]=1ll*inv[i]*((2-1ll*t[i]*inv[i]%mod+mod)%mod)%mod;
NTT(inv,len,-1);
fill(inv+n,inv+len,0);
}
void qiudao(const int *a,int *b,int n){
fill(b,b+n,0);
for(int i=1;i<n;i++)
b[i-1]=1ll*i*a[i]%mod;
}
void jifen(const int *a,int *b,int n){
fill(b,b+n,0);
for(int i=1;i<n;i++){
int inv=ksm(i,mod-2);
b[i]=1ll*a[i-1]*inv%mod;
}
}
void Mul(int *a,int *b,int len){
NTT(a,len,1);
NTT(b,len,1);
for(int i=0;i<len;i++)
a[i]=1ll*a[i]*b[i]%mod;
NTT(a,len,-1);
}
int n;
int a[N];
int b[N];
void Log(const int *a,int *b,int n){
static int ai[N],a1[N];
Inv(a,ai,n);
qiudao(a,a1,n);
int len=1;
while(len<n*2)
len<<=1;
fill(ai+n,ai+len,0);
fill(a1+n,a1+len,0);
Mul(a1,ai,len);
jifen(a1,b,n);
}
int main()
{
scanf("%d",&n);
for(int i=0;i<n;i++)
scanf("%d",&a[i]);
Log(a,b,n);
for(int i=0;i<n;i++)
printf("%d ",b[i]);
}
多项式的exp
设
l
n
(
B
(
x
)
)
≡
A
(
x
)
(
m
o
d
x
n
)
ln(B(x))\equiv A(x)(mod\ x^n)
ln(B(x))≡A(x)(mod xn)
也就是说
l
n
(
B
(
x
)
)
−
A
(
x
)
≡
0
(
m
o
d
x
n
)
ln(B(x))-A(x)\equiv0(mod\ x^n)
ln(B(x))−A(x)≡0(mod xn)
由牛顿迭代法,有
B
(
x
)
≡
B
0
(
x
)
−
l
n
(
B
0
(
x
)
−
A
(
x
)
1
B
0
(
x
)
≡
B
0
(
x
)
(
1
−
l
n
(
B
0
(
x
)
)
+
A
(
x
)
)
(
m
o
d
x
n
)
B(x)\equiv B_0(x)-\frac{ln(B_0(x)-A(x)}{\frac{1}{B_0(x)}}\equiv B_0(x)(1-ln(B_0(x))+A(x))(mod\ x^n)
B(x)≡B0(x)−B0(x)1ln(B0(x)−A(x)≡B0(x)(1−ln(B0(x))+A(x))(mod xn)
void Exp(const int *a,int *b,int n){
static int t[N];
if(n==1){
b[0]=1;
return ;
}
Exp(a,b,(n+1)/2);
Ln(b,t,n);
for(int i=0;i<n;i++)
t[i]=a[i]-t[i];
t[0]++;
int len=1;
while(len<2*n)
len<<=1;
fill(b+n,b+len,0);
fill(t+n,t+len,0);
Mul(b,t,len);
}
多项式的快速幂
当然,更容易想到的方法是用多项式乘法和多项式取模,然后像正常整数快速幂那样去做
然而这里是另外一种方法
设
B
(
x
)
=
A
k
(
x
)
B(x)=A^k(x)
B(x)=Ak(x)
那么
l
n
(
B
(
x
)
)
=
k
l
n
(
A
(
x
)
)
ln(B(x))=kln(A(x))
ln(B(x))=kln(A(x))
于是我们求出ln(A(x)),乘上k,再exp即可
void Pow(const int *a,int *b,int n,int k){
static int t[N];
Ln(a,t,n);
for(int i=0;i<n;i++)
t[i]=1ll*t[i]*k%mod;
Exp(t,b,n);
}