多项式求逆
定义
对于一个多项式A(x),如果存在另一个多项式B(x),有deg \ B(x) \leq deg \ A(x),且A(x)B(x) \equiv 1 (\mod x^{n}),那么称B(x)为A(x)在\mod x^{n}下的逆元,记为A^{-1}(x)
求法
当n=1时,A(x) \equiv c (\mod x),此时A^{-1}(x) \equiv c^{-1}
不妨设A(x)B^{'}(x) \equiv 1 (\mod x^{\lceil \frac{n}{2} \rceil}),A(x)B(x) \equiv 1 (\mod x^{n})
显然,也有A(x)B(x) \equiv 1 (\mod x^{\lceil \frac{n}{2} \rceil})
和第一个式子相减,可得B(x) \equiv B^{'}(x) (\mod x^{\lceil \frac{n}{2} \rceil})
移项,两边平方,有B^{2}(x)-2B(x)B^{'}(x)+B^{'2}(x) \equiv 0 (\mod x^{n})
两边同乘A(x),化简可得B(x) \equiv 2B^{'}(x)-A(x)B^{'2}(x) (\mod x^{n})
于是就可以根据最后一个式子递归计算了.在计算的时候需要拿NTT来实现多项式乘法的过程.总的时间复杂度为O(n\log n).
代码模板
IL void PolyInv(LL x[],LL y[],int len){
int i=0;
if (len==1) {
y[0]=Mi(x[0],MOD-2);
return;
}
PolyInv(x,y,len>>1);
for (i=0;i<len;i++) X[i]=x[i],Y[i]=y[i];
NTT(X,len<<1,1); NTT(Y,len<<1,1);
for (i=0;i<(len<<1);i++)
X[i]=((X[i]*Y[i])%MOD*Y[i])%MOD;
NTT(X,len<<1,-1);
for (i=0;i<len;i++) y[i]=((y[i]<<1)%MOD+MOD-X[i])%MOD;
}
例题
#include<cstring>
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<vector>
#include<cmath>
#include<map>
#define l(x) (x<<1)
#define r(x) ((x<<1)+1)
#define IL inline
#define reg register
#define LL long long
#define N 400010
#define MOD 998244353
#define INF 0x3f3f3f3f
using namespace std;
int n,i,len;
LL a[N],b[N],X[N],Y[N];
IL int Abs(int x){return (x<0)?-x:x;}
IL void Swap(LL &a,LL &b){a^=b^=a^=b;}
IL int Min(int a,int b){return (a<b)?a:b;}
IL int Max(int a,int b){return (a>b)?a:b;}
IL LL read(){
LL p=0,f=1; char c=getchar();
while (c<48||c>57) {if (c=='-') f=-1; c=getchar();}
while (c>=48&&c<=57) p=(p<<1)+(p<<3)+c-48,c=getchar();
return p*f;
}
IL LL Mi(LL x,LL y){
LL p=x,t=1,Res=1;
for (;t<=y;(t&y)?Res=(Res*p)%MOD:0,p=(p*p)%MOD,t<<=1);
return Res;
}
IL void Rader(LL y[],int len){
int i=0,j=0,k=0;
for (i=1,j=len>>1;i<len-1;i++){
if (i<j) Swap(y[i],y[j]);
k=len>>1;
while (j>=k) j-=k,k>>=1;
if (j<k) j+=k;
}
}
IL void NTT(LL a[],int len,int f){
int i=0,j=0,k=0;
LL wn=0,w=0,u=0,t=0;
Rader(a,len);
for (k=2;k<=len;k<<=1){
wn=Mi(3,(MOD-1)/k);
if (f==-1) wn=Mi(wn,MOD-2);
for (i=0;i<len;i+=k){
w=1;
for (j=i;j<i+k/2;j++){
u=a[j]; t=(w*a[j+k/2])%MOD;
a[j]=(u+t)%MOD; a[j+k/2]=(u-t+MOD)%MOD;
w=(w*wn)%MOD;
}
}
}
if (f==-1)
for (i=0;i<len;i++) a[i]=(a[i]*Mi(len,MOD-2))%MOD;
}
IL void PolyInv(LL x[],LL y[],int len){
int i=0;
if (len==1) {
y[0]=Mi(x[0],MOD-2);
return;
}
PolyInv(x,y,len>>1);
for (i=0;i<len;i++) X[i]=x[i],Y[i]=y[i];
NTT(X,len<<1,1); NTT(Y,len<<1,1);
for (i=0;i<(len<<1);i++)
X[i]=((X[i]*Y[i])%MOD*Y[i])%MOD;
NTT(X,len<<1,-1);
for (i=0;i<len;i++) y[i]=((y[i]<<1)%MOD+MOD-X[i])%MOD;
}
int main(){
#ifdef __Marvolo
freopen("zht.in","r",stdin);
freopen("zht.out","w",stdout);
#endif
n=read();
for (i=0;i<n;i++) a[i]=(read()+MOD)%MOD;
for (len=1;len<n;len<<=1);
PolyInv(a,b,len);
for (i=0;i<n;i++) printf("%lld ",b[i]);
return 0;
}
多项式除法&取模
给出两个多项式F(x),G(x),求D(x),R(x),使得F(x)=D(x)G(x)+R(x).其中,F(x)为n次多项式,G(x)为m次多项式,m \leq n.要求求出的D(x)为n-m次多项式.
首先定义翻转操作:对于一个n次多项式A(x),它的翻转多项式为A^{r}(x)=x^{n}A(\frac{1}{x}).假如说A(x)=x^{2}+2x+3,那么A^{r}(x)=3x^{2}+2x+1,也就是把系数翻转了一下.
定义了翻转操作后,对上面这个多项式除法式进行一下变形,将\frac{1}{x}替代x,两边同乘x^{n}.
得到x^{n}F(\frac{1}{x})=x^{m}G(\frac{1}{x})x^{n-m}D(\frac{1}{x})+x^{n-m+1}x^{m-1}R(\frac{1}{x})
化简,有F^{r}(x)=G^{r}(x)D^{r}(x)+x^{n-m+1}R^{r}(x)
两边同模x^{n-m+1},则R^{r}(x)这一项显然会被消掉,只剩下
F^{r}(x) \equiv G^{r}(x)D^{r}(x) (\mod x^{n-m+1})
已知D(x)的次数是n-m次,也就是说在上面的模意义下,D(x)的所有项都会保留下来.进一步变形,就有
D^{r}(x) \equiv F^{r}(x)G^{-1r}(x) (\mod x^{n-m+1})
然后利用上面的求逆元过程,算出D^{r}(x),翻转就可以得到D(x)了.然后带回到原式,就可以算出R(x).总的时间复杂度也是O(n \log n).
代码模板
IL void PolyMul(LL a[],LL b[],LL c[],int l1,int l2){
reg int i=0,L=Max(l1,l2),len=1;
for (;len<=L;len<<=1);
len<<=1;
for (i=0;i<=l1;i++) X[i]=a[i];
for (i=0;i<=l2;i++) Y[i]=b[i];
NTT(X,len,1); NTT(Y,len,1);
for (i=0;i<=len;i++) c[i]=(X[i]*Y[i])%MOD,X[i]=Y[i]=0;
NTT(c,len,-1);
}
IL void PolyDiv(LL x[],LL y[],LL a[],LL b[]){
reg int i=0,len=1;
reverse(x,x+1+n); reverse(y,y+1+m);
for (;len<=(n-m);len<<=1);
PolyInv(y,s,len);
memset(X,0,sizeof(X)); memset(Y,0,sizeof(Y));
PolyMul(x,s,a,n-m,n-m);
reverse(a,a+n-m+1); reverse(x,x+1+n); reverse(y,y+1+m);
PolyMul(a,y,b,n-m,m);
for (i=0;i<m;i++) b[i]=(f[i]-b[i]+MOD)%MOD;
}
例题
#include<cstring>
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<vector>
#include<cmath>
#include<map>
#define l(x) (x<<1)
#define r(x) ((x<<1)+1)
#define IL inline
#define reg register
#define LL long long
#define N 400010
#define MOD 998244353
#define INF 0x3f3f3f3f
using namespace std;
int n,m,i;
LL f[N],g[N],d[N],r[N],s[N],X[N],Y[N],t[N];
IL int Abs(int x){return (x<0)?-x:x;}
IL void Swap(LL &a,LL &b){a^=b^=a^=b;}
IL int Min(int a,int b){return (a<b)?a:b;}
IL int Max(int a,int b){return (a>b)?a:b;}
IL int read(){
int p=0,f=1; char c=getchar();
while (c<48||c>57) {if (c=='-') f=-1; c=getchar();}
while (c>=48&&c<=57) p=(p<<1)+(p<<3)+c-48,c=getchar();
return p*f;
}
IL LL Mi(LL x,LL y){
LL p=x,t=1,Res=1;
for (;t<=y;(t&y)?Res=(Res*p)%MOD:0,p=(p*p)%MOD,t<<=1);
return Res;
}
IL void Rader(LL y[],int len){
int i=0,j=0,k=0;
for (i=1,j=len>>1;i<len-1;i++){
if (i<j) Swap(y[i],y[j]);
k=len>>1;
while (j>=k) j-=k,k>>=1;
if (j<k) j+=k;
}
}
IL void NTT(LL a[],int len,int f){
reg int i=0,j=0,k=0;
LL wn=0,w=0,u=0,t=0;
Rader(a,len);
for (k=2;k<=len;k<<=1){
wn=Mi(3,(MOD-1)/k);
if (f==-1) wn=Mi(wn,MOD-2);
for (i=0;i<len;i+=k){
w=1;
for (j=i;j<i+k/2;j++){
u=a[j]; t=(w*a[j+k/2])%MOD;
a[j]=(u+t)%MOD; a[j+k/2]=(u-t+MOD)%MOD;
w=(w*wn)%MOD;
}
}
}
if (f==-1)
for (i=0;i<len;i++) a[i]=(a[i]*Mi(len,MOD-2))%MOD;
}
IL void PolyInv(LL x[],LL y[],int len){
reg int i=0;
if (len==1) {
y[0]=Mi(x[0],MOD-2);
return;
}
PolyInv(x,y,len>>1);
for (i=0;i<len;i++) X[i]=x[i],Y[i]=y[i];
NTT(X,len<<1,1); NTT(Y,len<<1,1);
for (i=0;i<(len<<1);i++)
X[i]=((X[i]*Y[i])%MOD*Y[i])%MOD;
NTT(X,len<<1,-1);
for (i=0;i<len;i++) y[i]=((y[i]<<1)%MOD+MOD-X[i])%MOD;
}
IL void PolyMul(LL a[],LL b[],LL c[],int l1,int l2){
reg int i=0,L=Max(l1,l2),len=1;
for (;len<=L;len<<=1);
len<<=1;
for (i=0;i<=l1;i++) X[i]=a[i];
for (i=0;i<=l2;i++) Y[i]=b[i];
NTT(X,len,1); NTT(Y,len,1);
for (i=0;i<=len;i++) c[i]=(X[i]*Y[i])%MOD,X[i]=Y[i]=0;
NTT(c,len,-1);
}
IL void PolyDiv(LL x[],LL y[],LL a[],LL b[]){
reg int i=0,len=1;
reverse(x,x+1+n); reverse(y,y+1+m);
for (;len<=(n-m);len<<=1);
PolyInv(y,s,len);
memset(X,0,sizeof(X)); memset(Y,0,sizeof(Y));
PolyMul(x,s,a,n-m,n-m);
reverse(a,a+n-m+1); reverse(x,x+1+n); reverse(y,y+1+m);
PolyMul(a,y,b,n-m,m);
for (i=0;i<m;i++) b[i]=(f[i]-b[i]+MOD)%MOD;
}
int main(){
#ifdef __Marvolo
freopen("zht.in","r",stdin);
freopen("zht.out","w",stdout);
#endif
n=read(); m=read();
for (i=0;i<=n;i++) f[i]=read();
for (i=0;i<=m;i++) g[i]=read();
PolyDiv(f,g,d,r);
for (i=0;i<=n-m;i++) printf("%lld ",d[i]);
putchar(10);
for (i=0;i<m;i++) printf("%lld ",r[i]);
return 0;
}
多项式开根
假如说求B(x),使得B(x)^{2} \equiv A(x) (\mod x^{n}).不妨设B'(x)^{2} \equiv A(x) (\mod x^{\lceil \frac{n}{2} \rceil}).同时,已知B(x)^{2} \equiv A(x) (\mod x^{\lceil \frac{n}{2} \rceil}).
两个等式相减,再平方,可得B'(x)^{4}-2B'(x)^{2}B(x)^{2}+B(x)^{4} \equiv 0 (\mod x^{n})
做一下变形,有B'(x)^{4}+2B'(x)^{2}B(x)^{2}+B(x)^{4} \equiv 4B'(x)^{2}B(x)^{2} (\mod x^{n})
故有B'(x)^{2}+B(x)^{2}\equiv 2B(x)B'(x) (\mod x^{n})
将已知条件代入,有B(x) \equiv \frac{A(x)+B'(x)^{2}}{2B'(x)}
和多项式求逆那样递归计算即可.
有一个问题是,当递归到n=1的时候,要求常数项在模意义下开根.洛谷上的例题规定了常数项一定为1,所以保证有解.如果不规定常数项的话,还需要通过二次剩余来判断解的存在性.
代码模板
注意数组大小要开到8倍以上.
IL void PolySqrt(LL a[],LL b[],int len){
reg int i=0;
if (len==1){
b[0]=1; return;
}
PolySqrt(a,b,(len+1)>>1);
memset(e,0,sizeof(e));
for (i=0;i<len;i++) c[i]=(b[i]<<1)%MOD;
PolyInv(c,e,len);
PolyMul(b,b,d,len,len);
for (i=0;i<len;i++) b[i]=(d[i]+a[i])%MOD;
PolyMul(b,e,b,len,len);
}
例题
#include<cstring>
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<vector>
#include<cmath>
#include<map>
#define l(x) (x<<1)
#define r(x) ((x<<1)+1)
#define IL inline
#define reg register
#define LL long long
#define N 800010
#define MOD 998244353
#define INF 0x3f3f3f3f
using namespace std;
int n,i,len;
LL f[N],g[N],X[N],Y[N],c[N],d[N],e[N],p[N];
IL int Abs(int x){return (x<0)?-x:x;}
IL void Swap(LL &a,LL &b){a^=b^=a^=b;}
IL int Min(int a,int b){return (a<b)?a:b;}
IL int Max(int a,int b){return (a>b)?a:b;}
IL LL read(){
LL p=0,f=1; char c=getchar();
while (c<48||c>57) {if (c=='-') f=-1; c=getchar();}
while (c>=48&&c<=57) p=(p<<1)+(p<<3)+c-48,c=getchar();
return p*f;
}
IL LL Mi(LL x,LL y){
LL p=x,t=1,Res=1;
for (;t<=y;(t&y)?Res=(Res*p)%MOD:0,p=(p*p)%MOD,t<<=1);
return Res;
}
IL void Rader(LL y[],int len){
int i=0,j=0,k=0;
for (i=1,j=len>>1;i<len-1;i++){
if (i<j) Swap(y[i],y[j]);
k=len>>1;
while (j>=k) j-=k,k>>=1;
if (j<k) j+=k;
}
}
IL void NTT(LL a[],int len,int f){
int i=0,j=0,k=0;
LL wn=0,w=0,u=0,t=0;
Rader(a,len);
for (k=2;k<=len;k<<=1){
wn=Mi(3,(MOD-1)/k);
if (f==-1) wn=Mi(wn,MOD-2);
for (i=0;i<len;i+=k){
w=1;
for (j=i;j<i+k/2;j++){
u=a[j]; t=(w*a[j+k/2])%MOD;
a[j]=(u+t)%MOD; a[j+k/2]=(u-t+MOD)%MOD;
w=(w*wn)%MOD;
}
}
}
t=Mi(len,MOD-2);
if (f==-1)
for (i=0;i<len;i++) a[i]=(a[i]*t)%MOD;
}
IL void PolyInv(LL x[],LL y[],int len){
int i=0;
if (len==1) {
y[0]=Mi(x[0],MOD-2);
return;
}
PolyInv(x,y,len>>1);
for (i=0;i<len;i++) X[i]=x[i],Y[i]=y[i];
NTT(X,len<<1,1); NTT(Y,len<<1,1);
for (i=0;i<(len<<1);i++)
X[i]=((X[i]*Y[i])%MOD*Y[i])%MOD;
NTT(X,len<<1,-1);
for (i=0;i<len;i++) y[i]=((y[i]<<1)%MOD+MOD-X[i])%MOD;
}
IL void PolyMul(LL a[],LL b[],LL c[],int l1,int l2){
reg int i=0,L=Max(l1,l2),len=1;
for (;len<=L;len<<=1);
len<<=1;
for (i=0;i<=l1;i++) X[i]=a[i];
for (i=0;i<=l2;i++) Y[i]=b[i];
NTT(X,len,1); NTT(Y,len,1);
for (i=0;i<=len;i++) c[i]=(X[i]*Y[i])%MOD,X[i]=Y[i]=0;
NTT(c,len,-1);
}
IL void PolySqrt(LL a[],LL b[],int len){
reg int i=0;
if (len==1){
b[0]=1; return;
}
PolySqrt(a,b,(len+1)>>1);
memset(e,0,sizeof(e));
for (i=0;i<len;i++) c[i]=(b[i]<<1)%MOD;
PolyInv(c,e,len);
PolyMul(b,b,d,len,len);
for (i=0;i<len;i++) b[i]=(d[i]+a[i])%MOD;
PolyMul(b,e,b,len,len);
}
int main(){
#ifdef __Marvolo
freopen("zht.in","r",stdin);
freopen("zht.out","w",stdout);
#endif
n=read();
for (i=0;i<n;i++) f[i]=read();
for (len=1;len<=n;len<<=1);
PolySqrt(f,g,len);
for (i=0;i<n;i++) printf("%lld ",g[i]);
return 0;
}
多项式求导&积分
原理应该不用多讲了...直接算就行了,积分的时候除上逆元.
代码模板
IL void PolyDx(LL a[],LL b[],LL len){
reg LL i=0;
for (i=0;i<=len;i++) b[i]=(a[i+1]*(i+1))%MOD;
}
IL void PolyInte(LL a[],LL b[],LL len){
reg LL i=0;
b[0]=0;
for (i=1;i<=len;i++)
b[i]=(a[i-1]*Mi(i,MOD-2))%MOD;
}
多项式Ln
给出A(x),求B(x),使得B(x) \equiv ln(A(x)) (\mod x^{n})
对原式两边求导,有B'(x) \equiv \frac{A'(x)}{A(x)} (\mod x^{n})
通过多项式的求导和求逆算出B'(x)后,再积分即可得到B(x).
代码模板
IL void PolyLn(LL a[],LL b[],LL len){
PolyDx(a,d,len);
PolyInv(a,c,len);
PolyMul(d,c,d,len,len);
PolyInte(d,b,len);
}
例题
#include<cstring>
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<vector>
#include<cmath>
#include<map>
#define l(x) (x<<1)
#define r(x) ((x<<1)+1)
#define IL inline
#define reg register
#define LL long long
#define N 800010
#define MOD 998244353
#define INF 0x3f3f3f3f
using namespace std;
int n,len,i;
LL f[N],g[N],c[N],d[N],X[N],Y[N];
IL int Abs(int x){return (x<0)?-x:x;}
IL void Swap(LL &a,LL &b){a^=b^=a^=b;}
IL int Min(int a,int b){return (a<b)?a:b;}
IL int Max(int a,int b){return (a>b)?a:b;}
IL LL read(){
LL p=0,f=1; char c=getchar();
while (c<48||c>57) {if (c=='-') f=-1; c=getchar();}
while (c>=48&&c<=57) p=(p<<1)+(p<<3)+c-48,c=getchar();
return p*f;
}
IL LL Mi(LL x,LL y){
LL p=x,t=1,Res=1;
for (;t<=y;(t&y)?Res=(Res*p)%MOD:0,p=(p*p)%MOD,t<<=1);
return Res;
}
IL void Rader(LL y[],int len){
int i=0,j=0,k=0;
for (i=1,j=len>>1;i<len-1;i++){
if (i<j) Swap(y[i],y[j]);
k=len>>1;
while (j>=k) j-=k,k>>=1;
if (j<k) j+=k;
}
}
IL void NTT(LL a[],int len,int f){
int i=0,j=0,k=0;
LL wn=0,w=0,u=0,t=0;
Rader(a,len);
for (k=2;k<=len;k<<=1){
wn=Mi(3,(MOD-1)/k);
if (f==-1) wn=Mi(wn,MOD-2);
for (i=0;i<len;i+=k){
w=1;
for (j=i;j<i+k/2;j++){
u=a[j]; t=(w*a[j+k/2])%MOD;
a[j]=(u+t)%MOD; a[j+k/2]=(u-t+MOD)%MOD;
w=(w*wn)%MOD;
}
}
}
t=Mi(len,MOD-2);
if (f==-1)
for (i=0;i<len;i++) a[i]=(a[i]*t)%MOD;
}
IL void PolyInv(LL x[],LL y[],int len){
int i=0;
if (len==1) {
y[0]=Mi(x[0],MOD-2);
return;
}
PolyInv(x,y,len>>1);
for (i=0;i<len;i++) X[i]=x[i],Y[i]=y[i];
NTT(X,len<<1,1); NTT(Y,len<<1,1);
for (i=0;i<(len<<1);i++)
X[i]=((X[i]*Y[i])%MOD*Y[i])%MOD;
NTT(X,len<<1,-1);
for (i=0;i<len;i++) y[i]=((y[i]<<1)%MOD+MOD-X[i])%MOD;
}
IL void PolyMul(LL a[],LL b[],LL c[],int l1,int l2){
reg int i=0,L=Max(l1,l2),len=1;
for (;len<=L;len<<=1);
len<<=1;
for (i=0;i<=l1;i++) X[i]=a[i];
for (i=0;i<=l2;i++) Y[i]=b[i];
NTT(X,len,1); NTT(Y,len,1);
for (i=0;i<=len;i++) c[i]=(X[i]*Y[i])%MOD,X[i]=Y[i]=0;
NTT(c,len,-1);
}
IL void PolyDx(LL a[],LL b[],LL len){
reg LL i=0;
for (i=0;i<=len;i++) b[i]=(a[i+1]*(i+1))%MOD;
}
IL void PolyInte(LL a[],LL b[],LL len){
reg LL i=0;
b[0]=0;
for (i=1;i<=len;i++)
b[i]=(a[i-1]*Mi(i,MOD-2))%MOD;
}
IL void PolyLn(LL a[],LL b[],LL len){
PolyDx(a,d,len);
PolyInv(a,c,len);
PolyMul(d,c,d,len,len);
PolyInte(d,b,len);
}
int main(){
#ifdef __Marvolo
freopen("zht.in","r",stdin);
freopen("zht.out","w",stdout);
#endif
n=read();
for (i=0;i<n;i++) f[i]=read();
for (len=1;len<=n;len<<=1);
PolyLn(f,g,len);
for (i=0;i<n;i++) printf("%lld ",g[i]);
return 0;
}
多项式Exp
需要一些多项式牛顿迭代的基础.
求一个多项式G(x),使得F(G(x)) \equiv 0 (\mod x^{n}).如果说已经求出来了F(G_{0}(x)) \equiv 0 (\mod x^{\lceil \frac{n}{2} \rceil}),将F(G(x))在G_{0}(x)处展开,因为G(x)-G_{0}(x)的最低次数已经是\lceil \frac{n}{2} \rceil了,所以这一项平方之后,在模x^{n}意义下为0.故泰勒展开只需要保留前两项,也就是:
F(G(x)) = F(G_{0}(x)) +(G(x)-G_{0}(x))F'(G_{0}(x)) \equiv 0 (\mod x^{n})
其中的F'(x)表示一阶导.化简,有G(x) \equiv G_{0}(x)-\frac{F(G_{0}(x))}{F'(G_{0}(x))} (\mod x^{n}).
然后回到这个问题上.现在是要求B(x) \equiv e^{A(x)} (\mod x^{n}).两边求对数,化简,有ln(B(x))-A(x) \equiv 0 (\mod x^{n}).
这个时候,构造F(G(x))=lnG(x)-A(x).利用牛顿迭代不断求解就行了.
式子的话,(F(G(x)))'=\frac{G'(x)}{G(x)}.把它带到迭代 的式子中,化简一下就有G(x)=G_{0}(x)-\frac{G_{0}(x)(1-lnG_{0}(x)+A(x))}{G'_{0}(x)}
最后根据上面这个式子就可以递归计算了.
因为要多次算Ln,所以中间用到的两个数组需要清空一下
不知道哪个环节有问题,这个板子的常数极大,而且是别人平均时间的两倍左右...有待优化.
IL void PolyLn(LL a[],LL b[],LL len){
PolyDx(a,d,len);
PolyInv(a,c,len);
PolyMul(d,c,d,len,len);
PolyInte(d,b,len);
for (LL i=0;i<=len;i++) c[i]=d[i]=0;
}
IL void PolyExp(LL a[],LL b[],LL len){
reg LL i=0;
if (len==1){
b[0]=1; return;
}
PolyExp(a,b,len>>1);
for (i=0;i<=len;i++) e[i]=0;
PolyLn(b,e,len);
for (i=0;i<len;i++)
e[i]=(((i==0)?1:0)-e[i]+a[i]+MOD)%MOD;
PolyMul(b,e,b,len,len);
}
例题
#include<cstring>
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<vector>
#include<cmath>
#include<map>
#define l(x) (x<<1)
#define r(x) ((x<<1)+1)
#define IL inline
#define reg register
#define LL long long
#define N 800010
#define MOD 998244353
#define INF 0x3f3f3f3f
using namespace std;
int n,i,len;
LL f[N],g[N],c[N],d[N],e[N],X[N],Y[N];
IL int Abs(int x){return (x<0)?-x:x;}
IL void Swap(LL &a,LL &b){a^=b^=a^=b;}
IL int Min(int a,int b){return (a<b)?a:b;}
IL int Max(int a,int b){return (a>b)?a:b;}
IL LL read(){
LL p=0,f=1; char c=getchar();
while (c<48||c>57) {if (c=='-') f=-1; c=getchar();}
while (c>=48&&c<=57) p=(p<<1)+(p<<3)+c-48,c=getchar();
return p*f;
}
IL LL Mi(LL x,LL y){
LL p=x,t=1,Res=1;
for (;t<=y;(t&y)?Res=(Res*p)%MOD:0,p=(p*p)%MOD,t<<=1);
return Res;
}
IL void Rader(LL y[],int len){
int i=0,j=0,k=0;
for (i=1,j=len>>1;i<len-1;i++){
if (i<j) Swap(y[i],y[j]);
k=len>>1;
while (j>=k) j-=k,k>>=1;
if (j<k) j+=k;
}
}
IL void NTT(LL a[],int len,int f){
int i=0,j=0,k=0;
LL wn=0,w=0,u=0,t=0;
Rader(a,len);
for (k=2;k<=len;k<<=1){
wn=Mi(3,(MOD-1)/k);
if (f==-1) wn=Mi(wn,MOD-2);
for (i=0;i<len;i+=k){
w=1;
for (j=i;j<i+k/2;j++){
u=a[j]; t=(w*a[j+k/2])%MOD;
a[j]=(u+t)%MOD; a[j+k/2]=(u-t+MOD)%MOD;
w=(w*wn)%MOD;
}
}
}
t=Mi(len,MOD-2);
if (f==-1)
for (i=0;i<len;i++) a[i]=(a[i]*t)%MOD;
}
IL void PolyInv(LL x[],LL y[],int len){
int i=0;
if (len==1) {
y[0]=Mi(x[0],MOD-2);
return;
}
PolyInv(x,y,len>>1);
for (i=0;i<len;i++) X[i]=x[i],Y[i]=y[i];
NTT(X,len<<1,1); NTT(Y,len<<1,1);
for (i=0;i<(len<<1);i++)
X[i]=((X[i]*Y[i])%MOD*Y[i])%MOD;
NTT(X,len<<1,-1);
for (i=0;i<len;i++) y[i]=((y[i]<<1)%MOD+MOD-X[i])%MOD;
}
IL void PolyMul(LL a[],LL b[],LL c[],int l1,int l2){
reg int i=0,L=Max(l1,l2),len=1;
for (;len<=L;len<<=1);
len<<=1;
for (i=0;i<=l1;i++) X[i]=a[i];
for (i=0;i<=l2;i++) Y[i]=b[i];
NTT(X,len,1); NTT(Y,len,1);
for (i=0;i<=len;i++) c[i]=(X[i]*Y[i])%MOD,X[i]=Y[i]=0;
NTT(c,len,-1);
}
IL void PolyDx(LL a[],LL b[],LL len){
reg LL i=0;
for (i=0;i<=len;i++) b[i]=(a[i+1]*(i+1))%MOD;
}
IL void PolyInte(LL a[],LL b[],LL len){
reg LL i=0;
b[0]=0;
for (i=1;i<=len;i++)
b[i]=(a[i-1]*Mi(i,MOD-2))%MOD;
}
IL void PolyLn(LL a[],LL b[],LL len){
PolyDx(a,d,len);
PolyInv(a,c,len);
PolyMul(d,c,d,len,len);
PolyInte(d,b,len);
for (LL i=0;i<=len;i++) c[i]=d[i]=0;
}
IL void PolyExp(LL a[],LL b[],LL len){
reg LL i=0;
if (len==1){
b[0]=1; return;
}
PolyExp(a,b,len>>1);
memset(e,0,sizeof(e));
PolyLn(b,e,len);
for (i=0;i<len;i++)
e[i]=(((i==0)?1:0)-e[i]+a[i]+MOD)%MOD;
PolyMul(b,e,b,len,len);
}
int main(){
#ifdef __Marvolo
freopen("zht.in","r",stdin);
freopen("zht.out","w",stdout);
#endif
n=read();
for (i=0;i<n;i++) f[i]=read();
for (len=1;len<=n;len<<=1);
PolyExp(f,g,len);
for (i=0;i<n;i++) printf("%lld ",g[i]);
return 0;
}
多项式快速幂
会了求导和求对数之后就很简单了.
B(x) \equiv A(x)^{k},求导,有lnB(x) \equiv klnA(x).
求导,乘以k,再exp回去就行了.
代码模板
因为同时使用了对数和指数的原因,常数极大.
IL void PolyMi(LL a[],LL b[],LL len,LL k){
reg int i=0;
PolyLn(a,p,len);
for (i=0;i<=len;i++) p[i]=(p[i]*k)%MOD;
PolyExp(p,b,len);
}
例题
#include<cstring>
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<vector>
#include<cmath>
#include<map>
#define l(x) (x<<1)
#define r(x) ((x<<1)+1)
#define IL inline
#define reg register
#define LL long long
#define N 800010
#define MOD 998244353
#define INF 0x3f3f3f3f
using namespace std;
LL n,i,k,len;
LL f[N],g[N],c[N],d[N],e[N],X[N],Y[N],p[N];
IL int Abs(int x){return (x<0)?-x:x;}
IL void Swap(LL &a,LL &b){a^=b^=a^=b;}
IL int Min(int a,int b){return (a<b)?a:b;}
IL int Max(int a,int b){return (a>b)?a:b;}
IL LL read(){
LL p=0,f=1; char c=getchar();
while (c<48||c>57) {if (c=='-') f=-1; c=getchar();}
while (c>=48&&c<=57) p=((p<<1)+(p<<3)+c-48)%MOD,c=getchar();
return p*f;
}
IL LL Mi(LL x,LL y){
LL p=x,t=1,Res=1;
for (;t<=y;(t&y)?Res=(Res*p)%MOD:0,p=(p*p)%MOD,t<<=1);
return Res;
}
IL void Rader(LL y[],int len){
int i=0,j=0,k=0;
for (i=1,j=len>>1;i<len-1;i++){
if (i<j) Swap(y[i],y[j]);
k=len>>1;
while (j>=k) j-=k,k>>=1;
if (j<k) j+=k;
}
}
IL void NTT(LL a[],int len,int f){
int i=0,j=0,k=0;
LL wn=0,w=0,u=0,t=0;
Rader(a,len);
for (k=2;k<=len;k<<=1){
wn=Mi(3,(MOD-1)/k);
if (f==-1) wn=Mi(wn,MOD-2);
for (i=0;i<len;i+=k){
w=1;
for (j=i;j<i+k/2;j++){
u=a[j]; t=(w*a[j+k/2])%MOD;
a[j]=(u+t)%MOD; a[j+k/2]=(u-t+MOD)%MOD;
w=(w*wn)%MOD;
}
}
}
t=Mi(len,MOD-2);
if (f==-1)
for (i=0;i<len;i++) a[i]=(a[i]*t)%MOD;
}
IL void PolyInv(LL x[],LL y[],int len){
int i=0;
if (len==1) {
y[0]=Mi(x[0],MOD-2);
return;
}
PolyInv(x,y,len>>1);
for (i=0;i<len;i++) X[i]=x[i],Y[i]=y[i];
NTT(X,len<<1,1); NTT(Y,len<<1,1);
for (i=0;i<(len<<1);i++)
X[i]=((X[i]*Y[i])%MOD*Y[i])%MOD;
NTT(X,len<<1,-1);
for (i=0;i<len;i++) y[i]=((y[i]<<1)%MOD+MOD-X[i])%MOD;
}
IL void PolyMul(LL a[],LL b[],LL c[],int l1,int l2){
reg int i=0,L=Max(l1,l2),len=1;
for (;len<=L;len<<=1);
len<<=1;
for (i=0;i<=l1;i++) X[i]=a[i];
for (i=0;i<=l2;i++) Y[i]=b[i];
NTT(X,len,1); NTT(Y,len,1);
for (i=0;i<=len;i++) c[i]=(X[i]*Y[i])%MOD,X[i]=Y[i]=0;
NTT(c,len,-1);
}
IL void PolyDx(LL a[],LL b[],LL len){
reg LL i=0;
for (i=0;i<=len;i++) b[i]=(a[i+1]*(i+1))%MOD;
}
IL void PolyInte(LL a[],LL b[],LL len){
reg LL i=0;
b[0]=0;
for (i=1;i<=len;i++)
b[i]=(a[i-1]*Mi(i,MOD-2))%MOD;
}
IL void PolyLn(LL a[],LL b[],LL len){
PolyDx(a,d,len);
PolyInv(a,c,len);
PolyMul(d,c,d,len,len);
PolyInte(d,b,len);
for (LL i=0;i<=len;i++) c[i]=d[i]=0;
}
IL void PolyExp(LL a[],LL b[],LL len){
reg LL i=0;
if (len==1){
b[0]=1; return;
}
PolyExp(a,b,len>>1);
for (i=0;i<=len;i++) e[i]=0;
PolyLn(b,e,len);
for (i=0;i<len;i++)
e[i]=(((i==0)?1:0)-e[i]+a[i]+MOD)%MOD;
PolyMul(b,e,b,len,len);
}
IL void PolyMi(LL a[],LL b[],LL len,LL k){
reg int i=0;
PolyLn(a,p,len);
for (i=0;i<=len;i++) p[i]=(p[i]*k)%MOD;
PolyExp(p,b,len);
}
int main(){
#ifdef __Marvolo
freopen("zht.in","r",stdin);
freopen("zht.out","w",stdout);
#endif
n=read(); k=read();
for (i=0;i<n;i++) f[i]=read();
for (len=1;len<=n;len<<=1);
PolyMi(f,g,len,k);
for (i=0;i<n;i++) printf("%lld ",g[i]);
return 0;
}
多项式三角函数
根据著名的欧拉公式:e^{ix}=\cos x+i\sin x.代入正负x,然后解方程,有:
\sin x=\frac{e^{ix}-e^{-ix}}{2i}
\cos x=\frac{e^{ix}+e^{-ix}}{2}
然后就是一波玄学操作:i^{2}=-1,i^{2} \equiv -1 (\mod 998244353),i \equiv 86583718 (\mod 998244353),i^{-1} \equiv 954952494 (\mod 998244353)
分子分母上的i就搞定了.
没学过复变,但是感觉复变老师看到这一段可能想打人
代码模板
这里参考大佬的博客给NTT加了一个预处理的优化,但是效果不是特别理想.
IL void Ready(){
reg LL i=0;
for (i=2;i<=800000;i<<=1)
Wn[i]=Mi(3,(MOD-1)/i),WN[i]=Mi(Wn[i],MOD-2);
}
所以还是依靠O2,O3优化靠谱一些.慎重使用这个板子吧.
IL void PolySin(LL a[],LL b[],LL len){
reg LL i=0,u=86583718,invu=Mi(u,MOD-2),inv2=Mi(2,MOD-2),inv=0;
for (i=0;i<=len;i++) a[i]=(a[i]*u)%MOD;
PolyExp(a,p,len); PolyInv(p,b,len);
inv=(inv2*invu)%MOD;
for (i=0;i<=len;i++) b[i]=(p[i]-b[i]+MOD)%MOD;
for (i=0;i<=len;i++) b[i]=(b[i]*inv)%MOD;
}
IL void PolyCos(LL a[],LL b[],LL len){
reg LL i=0,u=86583718,invu=Mi(u,MOD-2),inv2=Mi(2,MOD-2);
for (i=0;i<=len;i++) a[i]=(a[i]*u)%MOD;
PolyExp(a,p,len); PolyInv(p,b,len);
for (i=0;i<=len;i++) b[i]=(p[i]+b[i])%MOD;
for (i=0;i<=len;i++) b[i]=(b[i]*inv2)%MOD;
}
例题
#include<cstring>
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<vector>
#include<cmath>
#include<map>
#define l(x) (x<<1)
#define r(x) ((x<<1)+1)
#define IL inline
#define reg register
#define LL long long
#define N 800010
#define MOD 998244353
#define INF 0x3f3f3f3f
using namespace std;
LL n,k,i,len;
LL Wn[N],WN[N];
LL f[N],g[N],c[N],d[N],e[N],X[N],Y[N],p[N],q[N];
IL int Abs(int x){return (x<0)?-x:x;}
IL void Swap(LL &a,LL &b){a^=b^=a^=b;}
IL int Min(int a,int b){return (a<b)?a:b;}
IL int Max(int a,int b){return (a>b)?a:b;}
IL LL read(){
LL p=0,f=1; char c=getchar();
while (c<48||c>57) {if (c=='-') f=-1; c=getchar();}
while (c>=48&&c<=57) p=(p<<1)+(p<<3)+c-48,c=getchar();
return p*f;
}
IL LL Mi(LL x,LL y){
LL p=x,t=1,Res=1;
for (;t<=y;(t&y)?Res=(Res*p)%MOD:0,p=(p*p)%MOD,t<<=1);
return Res;
}
IL void Ready(){
reg LL i=0;
for (i=2;i<=800000;i<<=1)
Wn[i]=Mi(3,(MOD-1)/i),WN[i]=Mi(Wn[i],MOD-2);
}
IL void Rader(LL y[],int len){
int i=0,j=0,k=0;
for (i=1,j=len>>1;i<len-1;i++){
if (i<j) Swap(y[i],y[j]);
k=len>>1;
while (j>=k) j-=k,k>>=1;
if (j<k) j+=k;
}
}
IL void NTT(LL a[],int len,int f){
int i=0,j=0,k=0;
LL wn=0,w=0,u=0,t=0;
Rader(a,len);
for (k=2;k<=len;k<<=1){
wn=Wn[k];
if (f==-1) wn=WN[k];
for (i=0;i<len;i+=k){
w=1;
for (j=i;j<i+k/2;j++){
u=a[j]; t=(w*a[j+k/2])%MOD;
a[j]=(u+t)%MOD; a[j+k/2]=(u-t+MOD)%MOD;
w=(w*wn)%MOD;
}
}
}
t=Mi(len,MOD-2);
if (f==-1)
for (i=0;i<len;i++) a[i]=(a[i]*t)%MOD;
}
IL void PolyInv(LL x[],LL y[],int len){
int i=0;
if (len==1) {
y[0]=Mi(x[0],MOD-2);
return;
}
PolyInv(x,y,len>>1);
for (i=0;i<len;i++) X[i]=x[i],Y[i]=y[i];
NTT(X,len<<1,1); NTT(Y,len<<1,1);
for (i=0;i<(len<<1);i++)
X[i]=((X[i]*Y[i])%MOD*Y[i])%MOD;
NTT(X,len<<1,-1);
for (i=0;i<len;i++) y[i]=((y[i]<<1)%MOD+MOD-X[i])%MOD;
}
IL void PolyMul(LL a[],LL b[],LL c[],int l1,int l2){
reg int i=0,L=Max(l1,l2),len=1;
for (;len<=L;len<<=1);
len<<=1;
for (i=0;i<=l1;i++) X[i]=a[i];
for (i=0;i<=l2;i++) Y[i]=b[i];
NTT(X,len,1); NTT(Y,len,1);
for (i=0;i<=len;i++) c[i]=(X[i]*Y[i])%MOD,X[i]=Y[i]=0;
NTT(c,len,-1);
}
IL void PolyDx(LL a[],LL b[],LL len){
reg LL i=0;
for (i=0;i<=len;i++) b[i]=(a[i+1]*(i+1))%MOD;
}
IL void PolyInte(LL a[],LL b[],LL len){
reg LL i=0;
b[0]=0;
for (i=1;i<=len;i++)
b[i]=(a[i-1]*Mi(i,MOD-2))%MOD;
}
IL void PolyLn(LL a[],LL b[],LL len){
PolyDx(a,d,len);
PolyInv(a,c,len);
PolyMul(d,c,d,len,len);
PolyInte(d,b,len);
for (LL i=0;i<=len;i++) c[i]=d[i]=0;
}
IL void PolyExp(LL a[],LL b[],LL len){
reg LL i=0;
if (len==1){
b[0]=1; return;
}
PolyExp(a,b,len>>1);
for (i=0;i<=len;i++) e[i]=0;
PolyLn(b,e,len);
for (i=0;i<len;i++)
e[i]=(((i==0)?1:0)-e[i]+a[i]+MOD)%MOD;
PolyMul(b,e,b,len,len);
}
IL void PolySin(LL a[],LL b[],LL len){
reg LL i=0,u=86583718,invu=Mi(u,MOD-2),inv2=Mi(2,MOD-2),inv=0;
for (i=0;i<=len;i++) a[i]=(a[i]*u)%MOD;
PolyExp(a,p,len); PolyInv(p,b,len);
inv=(inv2*invu)%MOD;
for (i=0;i<=len;i++) b[i]=(p[i]-b[i]+MOD)%MOD;
for (i=0;i<=len;i++) b[i]=(b[i]*inv)%MOD;
}
IL void PolyCos(LL a[],LL b[],LL len){
reg LL i=0,u=86583718,invu=Mi(u,MOD-2),inv2=Mi(2,MOD-2);
for (i=0;i<=len;i++) a[i]=(a[i]*u)%MOD;
PolyExp(a,p,len); PolyInv(p,b,len);
for (i=0;i<=len;i++) b[i]=(p[i]+b[i])%MOD;
for (i=0;i<=len;i++) b[i]=(b[i]*inv2)%MOD;
}
int main(){
#ifdef __Marvolo
freopen("zht.in","r",stdin);
freopen("zht.out","w",stdout);
#endif
Ready();
n=read(); k=read();
for (i=0;i<n;i++) f[i]=read();
for (len=1;len<=n;len<<=1);
(k)?PolyCos(f,g,len):PolySin(f,g,len);
for (i=0;i<n;i++) printf("%lld ",g[i]);
return 0;
}
多项式反三角函数
首先,有
\frac{d}{dx}\arcsin x=\frac{1}{\sqrt{1-x^{2}}}
\frac{d}{dx}\arctan x=\frac{1}{1+x^{2}}
根据这个式子,将多项式代入,再积分,就有
B(x)=\int \frac{A'(x)}{\sqrt{1-A'(x)^{2}}}
B(x)=\int \frac{A'(x)}{1+A'(x)^{2}}
所以只需要多项式求导,求逆,开根,积分就能完成反三角函数的求解.
代码模板
IL void PolyArcsin(LL a[],LL b[],LL len){
reg int i=0;
PolyDx(a,p,len);
PolyMul(a,a,q,len,len);
for (i=0;i<=len;i++) q[i]=(MOD-q[i])%MOD;
q[0]=Upd(q[0]+1-MOD);
PolySqrt(q,w,len);
memset(q,0,sizeof(q));
PolyInv(w,q,len); PolyMul(p,q,p,len,len);
PolyInte(p,b,len);
}
IL void PolyArctan(LL a[],LL b[],LL len){
PolyDx(a,p,len);
PolyMul(a,a,q,len,len);
q[0]=(q[0]+1)%MOD;
PolyInv(q,w,len); PolyMul(p,w,p,len,len);
PolyInte(p,b,len);
}
例题
存在明显的卡常情况,最后的解决办法是优化加法取模的过程.
#pragma GCC optimize(3)
#include<cstring>
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<vector>
#include<cmath>
#include<map>
#define l(x) (x<<1)
#define r(x) ((x<<1)+1)
#define IL inline
#define reg register
#define LL long long
#define N 800010
#define INF 0x3f3f3f3f
using namespace std;
LL n,k,i,len;
const LL MOD=998244353;
LL Wn[N],WN[N];
LL f[N],g[N],c[N],d[N],e[N],X[N],Y[N],p[N],q[N],w[N];
IL int Abs(int x){return (x<0)?-x:x;}
IL void Swap(LL &a,LL &b){a^=b^=a^=b;}
IL int Min(int a,int b){return (a<b)?a:b;}
IL int Max(int a,int b){return (a>b)?a:b;}
char buf[1<<15],*fs,*ft;
inline char gc(){
return (fs==ft&&(ft=(fs=buf)+fread(buf,1,1<<15,stdin),fs==ft))?0:*fs++;}
inline int gint(){
int x=0,ch=gc();
while(ch<'0'||ch>'9') ch=gc();
while(ch>='0'&&ch<='9') x=(x<<1)+(x<<3)+ch-'0',ch=gc();
return x;
} //Fastest Read
IL LL Upd(int x) {
return x+(x>>31&MOD);
}
IL LL Mi(LL x,LL y){
reg LL p=x,t=1,Res=1;
for (;t<=y;(t&y)?Res=(Res*p)%MOD:0,p=(p*p)%MOD,t<<=1);
return Res;
}
IL void Ready(){
reg LL i=0;
for (i=2;i<=800000;i<<=1)
Wn[i]=Mi(3,(MOD-1)/i),WN[i]=Mi(Wn[i],MOD-2);
}
IL void NTT(LL a[],int len,int f){
reg int i=0,j=0,k=0;
reg LL wn=0,w=0,u=0,t=0;
for (i=1,j=len>>1;i<len-1;i++){
if (i<j) Swap(a[i],a[j]);
k=len>>1;
while (j>=k) j-=k,k>>=1;
(j<k)?j+=k:0;
}
for (k=2;k<=len;k<<=1){
wn=(f==-1)?WN[k]:Wn[k];
for (i=0;i<len;i+=k){
w=1;
for (j=i;j<i+k/2;j++){
u=a[j]; t=(w*a[j+k/2])%MOD;
a[j]=(u+t)%MOD; a[j+k/2]=(u-t+MOD)%MOD;
w=(w*wn)%MOD;
}
}
}
t=Mi(len,MOD-2);
if (f==-1)
for (i=0;i<len;i++) a[i]=(a[i]*t)%MOD;
}
IL void PolyInv(LL x[],LL y[],int len){
reg int i=0;
if (len==1) {
y[0]=Mi(x[0],MOD-2);
return;
}
PolyInv(x,y,len>>1);
for (i=0;i<len;i++) X[i]=x[i],Y[i]=y[i];
NTT(X,len<<1,1); NTT(Y,len<<1,1);
for (i=0;i<(len<<1);i++)
X[i]=((X[i]*Y[i])%MOD*Y[i])%MOD;
NTT(X,len<<1,-1);
for (i=0;i<len;i++) y[i]=((y[i]<<1)%MOD+MOD-X[i])%MOD;
}
IL void PolyMul(LL a[],LL b[],LL c[],int l1,int l2){
reg int i=0,L=Max(l1,l2),len=1;
for (;len<=L;len<<=1); len<<=1;
for (i=0;i<=L;i++) X[i]=a[i],Y[i]=b[i];
NTT(X,len,1); NTT(Y,len,1);
for (i=0;i<=len;i++) c[i]=(X[i]*Y[i])%MOD,X[i]=Y[i]=0;
NTT(c,len,-1);
}
IL void PolyDx(LL a[],LL b[],LL len){
reg LL i=0;
for (i=0;i<=len;i++) b[i]=(a[i+1]*(i+1))%MOD;
}
IL void PolyInte(LL a[],LL b[],LL len){
reg LL i=0;
b[0]=0;
for (i=1;i<=len;i++)
b[i]=(a[i-1]*Mi(i,MOD-2))%MOD;
}
IL void PolySqrt(LL a[],LL b[],int len){
reg int i=0;
if (len==1){
b[0]=1; return;
}
PolySqrt(a,b,(len+1)>>1);
for (i=0;i<=len;i++) e[i]=0;
for (i=0;i<len;i++) c[i]=Upd((b[i]<<1)-MOD);
PolyInv(c,e,len);
PolyMul(b,b,d,len,len);
for (i=0;i<len;i++) b[i]=Upd(d[i]+a[i]-MOD);
PolyMul(b,e,b,len,len);
}
IL void PolyArcsin(LL a[],LL b[],LL len){
reg int i=0;
PolyDx(a,p,len);
PolyMul(a,a,q,len,len);
for (i=0;i<=len;i++) q[i]=(MOD-q[i])%MOD;
q[0]=Upd(q[0]+1-MOD);
PolySqrt(q,w,len);
memset(q,0,sizeof(q));
PolyInv(w,q,len); PolyMul(p,q,p,len,len);
PolyInte(p,b,len);
}
IL void PolyArctan(LL a[],LL b[],LL len){
PolyDx(a,p,len);
PolyMul(a,a,q,len,len);
q[0]=(q[0]+1)%MOD;
PolyInv(q,w,len); PolyMul(p,w,p,len,len);
PolyInte(p,b,len);
}
int main(){
#ifdef __Marvolo
freopen("zht.in","r",stdin);
freopen("zht.out","w",stdout);
#endif
Ready();
n=gint(); k=gint();
for (i=0;i<n;i++) f[i]=gint();
for (len=1;len<=n;len<<=1);
(k)?PolyArctan(f,g,len):PolyArcsin(f,g,len);
for (i=0;i<n;i++) printf("%lld ",g[i]);
return 0;
}
例题
题意:求
G(x) \equiv ((1+ln(2+F(x)-F(0)-exp(\int _{0}^{x}\frac{1}{\sqrt{F(x)}})))^{k})'
题解:
多项式大杂烩,套一下上面的板子就行了.如果不优化的话常数特别大,会T.强开O3后总时间可以卡到30s.用之前的NTT预处理优化+Upd优化后,不开O3时间可以压到28s.
代码:
#pragma GCC optimize(3)
#include<cstring>
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<vector>
#include<cmath>
#include<map>
#define l(x) (x<<1)
#define r(x) ((x<<1)+1)
#define IL inline
#define reg register
#define LL long long
#define N 800010
#define MOD 998244353
#define INF 0x3f3f3f3f
using namespace std;
LL n,i,len,re,k;
LL Wn[N],WN[N];
LL f[N],g[N],X[N],Y[N],c[N],d[N],e[N],p[N],q[N];
struct F{
LL a,b,w;
F(){};
F(const LL &x,const LL &y,const LL &c){a=x; b=y; w=c;}
};
IL int Abs(int x){return (x<0)?-x:x;}
IL void Swap(LL &a,LL &b){a^=b^=a^=b;}
IL LL Min(LL a,LL b){return (a<b)?a:b;}
IL int Max(int a,int b){return (a>b)?a:b;}
IL int read(){
int p=0,f=1; char c=getchar();
while (c<48||c>57) {if (c=='-') f=-1; c=getchar();}
while (c>=48&&c<=57) p=(p<<1)+(p<<3)+c-48,c=getchar();
return p*f;
}
IL LL Upd(int x) {
return x+(x>>31&MOD);
}
IL LL mi(LL x,LL y,LL mod){
LL p=x,t=1,Res=1;
for (;t<=y;(t&y)?Res=(Res*p)%mod:0,p=(p*p)%mod,t<<=1);
return Res;
}
IL LL Mi(LL x,LL y){
LL p=x,t=1,Res=1;
for (;t<=y;(t&y)?Res=(Res*p)%MOD:0,p=(p*p)%MOD,t<<=1);
return Res;
}
IL F Mul(F a,F b,LL mod){
F c=F(0,0,a.w);
c.a=((a.a*b.a)%mod+(((a.b*b.b)%mod)*a.w)%mod)%mod;
c.b=(b.a*a.b%mod+a.a*b.b%mod)%mod;
return c;
}
IL F Pow(F x,LL y,LL mod){
LL t=1;
F p=x,Res=F(1,0,x.w);
for (;t<=y;Res=(t&y)?Mul(Res,p,mod):Res,p=Mul(p,p,mod),t<<=1);
return Res;
}
IL LL Solve(int x,int p){
LL a=0,w=0;
if (!x) return 0;
if (mi(x,(p-1)>>1,p)==p-1) return -1;
while (1){
a=rand(); w=(a*a%p-x+p)%p;
if (mi(w,(p-1)>>1,p)==p-1){
F s=Pow(F(a,1,w),(p+1)>>1,p);
return s.a;
}
}
}
IL void Ready(){
reg LL i=0;
for (i=2;i<=800000;i<<=1)
Wn[i]=Mi(3,(MOD-1)/i),WN[i]=Mi(Wn[i],MOD-2);
}
IL void Rader(LL y[],int len){
int i=0,j=0,k=0;
for (i=1,j=len>>1;i<len-1;i++){
if (i<j) Swap(y[i],y[j]);
k=len>>1;
while (j>=k) j-=k,k>>=1;
if (j<k) j+=k;
}
}
IL void NTT(LL a[],int len,int f){
int i=0,j=0,k=0;
LL wn=0,w=0,u=0,t=0;
Rader(a,len);
for (k=2;k<=len;k<<=1){
wn=Wn[k];
if (f==-1) wn=WN[k];
for (i=0;i<len;i+=k){
w=1;
for (j=i;j<i+k/2;j++){
u=a[j]; t=(w*a[j+k/2])%MOD;
a[j]=Upd(u+t-MOD); a[j+k/2]=Upd(u-t);
w=(w*wn)%MOD;
}
}
}
t=Mi(len,MOD-2);
if (f==-1)
for (i=0;i<len;i++) a[i]=(a[i]*t)%MOD;
}
IL void PolyInv(LL x[],LL y[],int len){
int i=0;
if (len==1) {
y[0]=Mi(x[0],MOD-2);
return;
}
PolyInv(x,y,len>>1);
for (i=0;i<len;i++) X[i]=x[i],Y[i]=y[i];
NTT(X,len<<1,1); NTT(Y,len<<1,1);
for (i=0;i<(len<<1);i++)
X[i]=((X[i]*Y[i])%MOD*Y[i])%MOD;
NTT(X,len<<1,-1);
for (i=0;i<len;i++) y[i]=Upd(Upd((y[i]<<1)-MOD)-X[i]);
}
IL void PolyMul(LL a[],LL b[],LL c[],int l1,int l2){
reg int i=0,L=Max(l1,l2),len=1;
for (;len<=L;len<<=1);
len<<=1;
for (i=0;i<=l1;i++) X[i]=a[i];
for (i=0;i<=l2;i++) Y[i]=b[i];
NTT(X,len,1); NTT(Y,len,1);
for (i=0;i<=len;i++) c[i]=(X[i]*Y[i])%MOD,X[i]=Y[i]=0;
NTT(c,len,-1);
}
IL void PolySqrt(LL a[],LL b[],int len){
reg int i=0;
if (len==1){
b[0]=re; return;
}
PolySqrt(a,b,(len+1)>>1);
memset(e,0,sizeof(e));
for (i=0;i<len;i++) c[i]=Upd((b[i]<<1)-MOD);
PolyInv(c,e,len);
PolyMul(b,b,d,len,len);
for (i=0;i<len;i++) b[i]=Upd(d[i]+a[i]-MOD);
PolyMul(b,e,b,len,len);
}
IL void PolyDx(LL a[],LL b[],LL len){
reg LL i=0;
for (i=0;i<=len;i++) b[i]=(a[i+1]*(i+1))%MOD;
}
IL void PolyInte(LL a[],LL b[],LL len){
reg LL i=0;
b[0]=0;
for (i=1;i<=len;i++)
b[i]=(a[i-1]*Mi(i,MOD-2))%MOD;
}
IL void PolyLn(LL a[],LL b[],LL len){
PolyDx(a,d,len);
PolyInv(a,c,len);
PolyMul(d,c,d,len,len);
PolyInte(d,b,len);
for (LL i=0;i<=len;i++) c[i]=d[i]=0;
}
IL void PolyExp(LL a[],LL b[],LL len){
reg LL i=0;
if (len==1){
b[0]=1; return;
}
PolyExp(a,b,len>>1);
for (i=0;i<=len;i++) e[i]=0;
PolyLn(b,e,len);
for (i=0;i<len;i++)
e[i]=Upd(((i==0)?1:0)-e[i]+a[i]);
PolyMul(b,e,b,len,len);
}
IL void PolyMi(LL a[],LL b[],LL len,LL k){
reg int i=0;
PolyLn(a,p,len);
for (i=0;i<=len;i++) p[i]=(p[i]*k)%MOD;
PolyExp(p,b,len);
}
IL void Work(){
reg LL i=0;
PolySqrt(f,g,len);
memset(c,0,sizeof(c));
memset(d,0,sizeof(d));
PolyInv(g,p,len);
PolyInte(p,g,len);
memset(p,0,sizeof(p));
memset(X,0,sizeof(X));
memset(Y,0,sizeof(Y));
PolyExp(g,p,len);
memset(g,0,sizeof(g));
for (i=0;i<=len;i++) g[i]=(f[i]-p[i]+MOD)%MOD;
g[0]=(g[0]+2-f[0]+MOD)%MOD;
memset(p,0,sizeof(p));
PolyLn(g,p,len); p[0]=(p[0]+1)%MOD;
memset(g,0,sizeof(g));
PolyMi(p,g,len,k);
PolyDx(g,p,len);
}
int main(){
#ifdef __Marvolo
freopen("zht.in","r",stdin);
freopen("zht.out","w",stdout);
#endif
Ready();
n=read(); k=read(); ++n;
for (i=0;i<n;i++) f[i]=read();
for (len=1;len<=n;len<<=1);
re=Solve(f[0],MOD);
re=Min(re,MOD-re);
Work();
for (i=0;i<n-1;i++) printf("%lld ",p[i]);
return 0;
}
参考资料
[1] https://www.cnblogs.com/bztMinamoto/p/9743310.html
[2] http://blog.miskcoo.com/2015/05/polynomial-inverse
[3] https://www.cnblogs.com/bztMinamoto/p/9744682.html
[4] https://blog.csdn.net/ezoixx118/article/details/86763688
声明:
本文采用
BY-NC-SA
协议进行授权，如无注明均为原创，转载请注明转自
Marvolo
本文地址: 多项式全家桶学习笔记
本文地址: 多项式全家桶学习笔记
