一些约定:下面 f i ( x ) f^i(x) fi(x) 表示 i i i 阶导数。 f ( x ) i f(x)^i f(x)i 表示幂次。若不说明绝大部分除一个多项式时都是代表乘上它的逆。
多项式加减,求导积分
过于简单不讲。
(
x
a
)
′
=
a
x
a
−
1
,
∫
x
a
d
x
=
x
a
+
1
a
+
1
+
C
(x^a)'=ax^{a-1},\displaystyle\int x^a {\rm d}x=\dfrac{x^{a+1}}{a+1}+C
(xa)′=axa−1,∫xadx=a+1xa+1+C。实际操作时
C
=
0
C=0
C=0,正确性是好证的。(如果你实在闲可以每次证一下)
多项式乘法,任意模数多项式乘法
目前常数很小而且代码不长的 NTT:参考这里。
任意模数请自行学习,基本 useless。
普通 $\texttt{NTT code}$#include<bits/stdc++.h>
#define fr(x) freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);
using namespace std;
const int mod=998244353,N=4e6+5;
int n,m,a[N],b[N],w[N],mmax;
inline int rd()
{
int x=0,zf=1;
char ch=getchar();
while(ch<'0'||ch>'9') (ch=='-')and(zf=-1),ch=getchar();
while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
return x*zf;
}
inline void wr(int x)
{
if(x==0) return putchar('0'),putchar(' '),void();
int num[35],len=0;
while(x) num[++len]=x%10,x/=10;
for(int i=len;i>=1;i--) putchar(num[i]+'0');
putchar(' ');
}
inline int bger(int x){return x|=x>>1,x|=x>>2,x|=x>>4,x|=x>>8,x|=x>>16,x+1;}
inline int md(int x){return x>=mod?x-mod:x;}
inline int ksm(int x,int p){int s=1;for(;p;(p&1)&&(s=1ll*s*x%mod),x=1ll*x*x%mod,p>>=1);return s;}
inline void init(int mmax)
{
for(int i=1,j,k;i<mmax;i<<=1)
for(w[j=i]=1,k=ksm(3,(mod-1)/(i<<1)),j++;j<(i<<1);j++)
w[j]=1ll*w[j-1]*k%mod;
}
inline void DNT(int *a,int mmax)
{
for(int i,j,k=mmax>>1,L,*W,*x,*y,z;k;k>>=1)
for(L=k<<1,i=0;i<mmax;i+=L)
for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
*y=1ll*(*x+mod-(z=*y))* *W%mod,*x=md(*x+z);
}
inline void IDNT(int *a,int mmax)
{
for(int i,j,k=1,L,*W,*x,*y,z;k<mmax;k<<=1)
for(L=k<<1,i=0;i<mmax;i+=L)
for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
z=1ll* *W* *y%mod,*y=md(*x+mod),*x=md(*x);
reverse(a+1,a+mmax);
for(int inv=ksm(mmax,mod-2),i=0;i<mmax;i++) a[i]=1ll*a[i]*inv%mod;
}
inline void NTT(int *a,int *b,int n,int m)
{
mmax=bger(n+m);init(mmax);DNT(a,mmax);DNT(b,mmax);
for(int i=0;i<mmax;i++) a[i]=1ll*a[i]*b[i]%mod;IDNT(a,mmax);
}
int main()
{
n=rd();m=rd();
for(int i=0;i<=n;i++) a[i]=rd();
for(int i=0;i<=m;i++) b[i]=rd();
NTT(a,b,n,m);
for(int i=0;i<=n+m;i++) wr(a[i]);
return 0;
}
多项式牛顿迭代
给定多项式 g ( x ) g(x) g(x),已知 g ( f ( x ) ) ≡ 0 ( m o d x n ) g(f(x))\equiv 0\pmod{x^n} g(f(x))≡0(modxn),要解出 f ( x ) f(x) f(x)。其中 g ( x ) g(x) g(x) 的形式较为简单(具体看下面求逆)。(注:这没有模板,只是一个前置知识,会应用到各个模板)
考虑倍增。
若已经求出 f 0 ( x ) f_0(x) f0(x) 满足: g ( f 0 ( x ) ) ≡ 0 ( m o d x ⌈ n / 2 ⌉ ) g(f_0(x))\equiv 0\pmod{x^{\lceil n/2\rceil}} g(f0(x))≡0(modx⌈n/2⌉),下面考虑求出 f ( x ) f(x) f(x)。
把 g ( f ( x ) ) g(f(x)) g(f(x)) 在 f 0 ( x ) f_0(x) f0(x) 出泰勒展开:
0 ≡ g ( f ( x ) ) = ∑ i ≥ 0 1 i ! g i ( f 0 ( x ) ) ( f ( x ) − f 0 ( x ) ) i ( m o d x n ) 0\equiv g(f(x))=\sum\limits_{i\ge 0} \dfrac{1}{i!}g^i(f_0(x))(f(x)-f_0(x))^i\pmod {x^n} 0≡g(f(x))=i≥0∑i!1gi(f0(x))(f(x)−f0(x))i(modxn)。
注意到 f ( x ) − f 0 ( x ) f(x)-f_0(x) f(x)−f0(x) 的最低次项次数 ≥ ⌈ n / 2 ⌉ \ge \lceil n/2\rceil ≥⌈n/2⌉,于是 ∀ i ≥ 2 , ( f ( x ) − f 0 ( x ) ) i ≡ 0 ( m o d x n ) \forall i\ge 2,(f(x)-f_0(x))^i\equiv 0\pmod {x^n} ∀i≥2,(f(x)−f0(x))i≡0(modxn)。
于是 g ′ ( f 0 ( x ) ) ( f ( x ) − f 0 ( x ) ) + g ( f 0 ( x ) ) ≡ 0 ( m o d x n ) ⇒ f ( x ) ≡ f 0 ( x ) − g ( f 0 ( x ) ) g ′ ( f 0 ( x ) ) ( m o d x n ) g'(f_0(x))(f(x)-f_0(x))+g(f_0(x))\equiv 0\pmod {x^n}\Rightarrow f(x)\equiv f_0(x)-\dfrac{g(f_0(x))}{g'(f_0(x))}\pmod {x^n} g′(f0(x))(f(x)−f0(x))+g(f0(x))≡0(modxn)⇒f(x)≡f0(x)−g′(f0(x))g(f0(x))(modxn)。
你可能觉得这不是还要求逆吗?你为啥不先讲求逆?你先别急。
多项式求逆
给定一个多项式 f ( x ) f(x) f(x),请求出一个多项式 g ( x ) g(x) g(x),满足 f ( x ) × g ( x ) ≡ 1 ( m o d x n ) f(x)\times g(x)\equiv1\pmod{x^n} f(x)×g(x)≡1(modxn)。系数对 998244353 998244353 998244353 取模。
$\texttt{code}$考虑牛顿迭代。
G ( X ) = f ( x ) − 1 X G(X)=f(x)-\dfrac{1}{X} G(X)=f(x)−X1,则 G ( g ( x ) ) ≡ 0 ( m o d x n ) G(g(x))\equiv0 \pmod {x^n} G(g(x))≡0(modxn)。
假设求出了 g 0 ( x ) g_0(x) g0(x) 使得 G ( g 0 ( x ) ) ≡ 0 ( m o d x ⌈ n / 2 ⌉ ) G(g_0(x))\equiv 0\pmod{x^{\lceil n/2\rceil}} G(g0(x))≡0(modx⌈n/2⌉),则:
g ( x ) ≡ g 0 ( x ) − G ( g 0 ( x ) ) G ′ ( g 0 ( x ) ) = g 0 ( x ) − f ( x ) − 1 g 0 ( x ) 1 g 0 ( x ) 2 ≡ 2 g 0 ( x ) − f ( x ) g 0 ( x ) 2 ( m o d x n ) g(x)\equiv g_0(x)-\dfrac{G(g_0(x))}{G'(g_0(x))}=g_0(x)-\dfrac{f(x)-\dfrac{1}{g_0(x)}}{\dfrac{1}{g_0(x)^2}}\equiv 2g_0(x)-f(x)g_0(x)^2\pmod {x^n} g(x)≡g0(x)−G′(g0(x))G(g0(x))=g0(x)−g0(x)21f(x)−g0(x)1≡2g0(x)−f(x)g0(x)2(modxn)。
用两次多项式乘法即可,但是可以通过点值进行特殊变换做到一次多项式乘法。
即一般的多项式乘法设两个对应的点值为 x , y x,y x,y,则变换 f ( x , y ) = x y f(x,y)=xy f(x,y)=xy,这里 f ( x , y ) = 2 y − x y 2 f(x,y)=2y-xy^2 f(x,y)=2y−xy2 即可。
复杂度 T ( n ) = T ( n / 2 ) + O ( n log n ) = O ( n log n ) T(n)=T(n/2)+O(n\log n)=O(n\log n) T(n)=T(n/2)+O(nlogn)=O(nlogn)。
#include<bits/stdc++.h>
#define LL long long
using namespace std;
const int mod=998244353,N=4e5+5;
int n,m,a[N],b[N],w[N],mmax;
inline int rd()
{
int x=0,zf=1;
char ch=getchar();
while(ch<'0'||ch>'9') (ch=='-')and(zf=-1),ch=getchar();
while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
return x*zf;
}
inline void wr(int x)
{
if(x==0) return putchar('0'),putchar(' '),void();
int num[35],len=0;
while(x) num[++len]=x%10,x/=10;
for(int i=len;i>=1;i--) putchar(num[i]+'0');
putchar(' ');
}
inline int bger(int x){return x|=x>>1,x|=x>>2,x|=x>>4,x|=x>>8,x|=x>>16,x+1;}
inline int md(int x){return x>=mod?x-mod:x;}
inline int ksm(int x,int p){int s=1;for(;p;(p&1)&&(s=1ll*s*x%mod),x=1ll*x*x%mod,p>>=1);return s;}
inline void init(int mmax)
{
for(int i=1,j,k;i<mmax;i<<=1)
for(w[j=i]=1,k=ksm(3,(mod-1)/(i<<1)),j++;j<(i<<1);j++)
w[j]=1ll*w[j-1]*k%mod;
}
inline void DNT(int *a,int mmax)
{
for(int i,j,k=mmax>>1,L,*W,*x,*y,z;k;k>>=1)
for(L=k<<1,i=0;i<mmax;i+=L)
for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
*y=1ll*(*x+mod-(z=*y))* *W%mod,*x=md(*x+z);
}
inline void IDNT(int *a,int mmax)
{
for(int i,j,k=1,L,*W,*x,*y,z;k<mmax;k<<=1)
for(L=k<<1,i=0;i<mmax;i+=L)
for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
z=1ll* *W* *y%mod,*y=md(*x+mod-z),*x=md(*x+z);
reverse(a+1,a+mmax);
for(int inv=ksm(mmax,mod-2),i=0;i<mmax;i++) a[i]=1ll*a[i]*inv%mod;
}
void INV(int num,int *a,int *b)
{
if(num==1) return b[0]=ksm(a[0],mod-2),void();
INV((num+1)>>1,a,b);
int mmax=bger(num<<1);init(mmax);
static int c[N];
for(int i=0;i<num;i++) c[i]=a[i];for(int i=num;i<mmax;i++) c[i]=0;
DNT(c,mmax);DNT(b,mmax);
for(int i=0;i<mmax;i++) b[i]=1ll*(2-1ll*c[i]*b[i]%mod+mod)%mod*b[i]%mod;
IDNT(b,mmax);
for(int i=num;i<mmax;i++) b[i]=0;
}
signed main()
{
n=rd();
for(int i=0;i<n;i++) a[i]=rd();
INV(n,a,b);
for(int i=0;i<n;i++) wr(b[i]);
return 0;
}
分治 FFT
给定序列 g 1 , 2 … n − 1 g_{1,2\dots n - 1} g1,2…n−1,求序列 f 0 , 1 … n − 1 f_{0,1\dots n - 1} f0,1…n−1。
其中 f i = ∑ j = 1 i f i − j g j f_i=\sum\limits_{j=1}^if_{i-j}g_j fi=j=1∑ifi−jgj,边界为 f 0 = 1 f_0=1 f0=1。答案对 998244353 998244353 998244353 取模。
$\texttt{code}$设 F ( x ) = ∑ i ≥ 0 f i x i , G ( x ) = ∑ i ≥ 0 g i x i F(x)=\sum\limits_{i\ge 0}f_ix^i,G(x)=\sum\limits_{i\ge 0}g_ix^i F(x)=i≥0∑fixi,G(x)=i≥0∑gixi,其中令 g 0 = 0 g_0=0 g0=0。
则 F ( x ) × G ( x ) = ∑ i ≥ 0 x i ∑ j = 0 i f i − j g j = ∑ i ≥ 0 x i ∑ j = 0 i f i − j g j = ∑ i ≥ 1 x i f i = F ( x ) − f 0 = F ( x ) − 1 F(x)\times G(x)=\sum\limits_{i\ge 0}x^i\sum\limits_{j=0}^i f_{i-j}g_j=\sum\limits_{i\ge 0}x^i\sum\limits_{j=0}^i f_{i-j}g_j=\sum\limits_{i\ge 1}x^if_i=F(x)-f_0=F(x)-1 F(x)×G(x)=i≥0∑xij=0∑ifi−jgj=i≥0∑xij=0∑ifi−jgj=i≥1∑xifi=F(x)−f0=F(x)−1。
于是 G ( x ) = 1 1 − F ( x ) G(x)=\dfrac{1}{1-F(x)} G(x)=1−F(x)1。多项式求逆即可,复杂度 O ( n log n ) O(n\log n) O(nlogn)。
#include<bits/stdc++.h>
#define fr(x) freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);
using namespace std;
const int mod=998244353,N=4e5+5;
int n,m,a[N],b[N],w[N],mmax;
inline int rd()
{
int x=0,zf=1;
char ch=getchar();
while(ch<'0'||ch>'9') (ch=='-')and(zf=-1),ch=getchar();
while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
return x*zf;
}
inline void wr(int x)
{
if(x==0) return putchar('0'),putchar(' '),void();
int num[35],len=0;
while(x) num[++len]=x%10,x/=10;
for(int i=len;i>=1;i--) putchar(num[i]+'0');
putchar(' ');
}
inline int bger(int x){return x|=x>>1,x|=x>>2,x|=x>>4,x|=x>>8,x|=x>>16,x+1;}
inline int md(int x){return x>=mod?x-mod:x;}
inline int ksm(int x,int p){int s=1;for(;p;(p&1)&&(s=1ll*s*x%mod),x=1ll*x*x%mod,p>>=1);return s;}
inline void init(int mmax)
{
for(int i=1,j,k;i<mmax;i<<=1)
for(w[j=i]=1,k=ksm(3,(mod-1)/(i<<1)),j++;j<(i<<1);j++)
w[j]=1ll*w[j-1]*k%mod;
}
inline void DNT(int *a,int mmax)
{
for(int i,j,k=mmax>>1,L,*W,*x,*y,z;k;k>>=1)
for(L=k<<1,i=0;i<mmax;i+=L)
for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
*y=1ll*(*x+mod-(z=*y))* *W%mod,*x=md(*x+z);
}
inline void IDNT(int *a,int mmax)
{
for(int i,j,k=1,L,*W,*x,*y,z;k<mmax;k<<=1)
for(L=k<<1,i=0;i<mmax;i+=L)
for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
z=1ll* *W* *y%mod,*y=md(*x+mod-z),*x=md(*x+z);
reverse(a+1,a+mmax);
for(int inv=ksm(mmax,mod-2),i=0;i<mmax;i++) a[i]=1ll*a[i]*inv%mod;
}
void dfs(int num,int *a,int *b)
{
if(num==1) return b[0]=ksm(a[0],mod-2),void();
dfs((num+1)>>1,a,b);
int mmax=bger(num<<1);init(mmax);
static int c[N];
for(int i=0;i<num;i++) c[i]=a[i];for(int i=num;i<mmax;i++) c[i]=0;
DNT(c,mmax);DNT(b,mmax);
for(int i=0;i<mmax;i++) b[i]=1ll*(2-1ll*c[i]*b[i]%mod+mod)%mod*b[i]%mod;
IDNT(b,mmax);
for(int i=num;i<mmax;i++) b[i]=0;
}
signed main()
{
n=rd();
for(int i=1;i<n;i++) a[i]=(mod-rd())%mod;a[0]=1;
dfs(n,a,b);
for(int i=0;i<n;i++) wr(b[i]);
return 0;
}
多项式除法
给定一个 n n n 次多项式 F ( x ) F(x) F(x) 和一个 m m m 次多项式 G ( x ) G(x) G(x) ,请求出多项式 Q ( x ) Q(x) Q(x), R ( x ) R(x) R(x),满足以下条件:
- Q ( x ) Q(x) Q(x) 次数为 n − m n-m n−m, R ( x ) R(x) R(x) 次数小于 m m m
- F ( x ) = Q ( x ) × G ( x ) + R ( x ) F(x) = Q(x) \times G(x) + R(x) F(x)=Q(x)×G(x)+R(x)
所有的运算在模 998244353 998244353 998244353 意义下进行。保证 m < n m<n m<n。
$\texttt{code}$这种 trick 基本不会用第二次了。注意到我们只需求出 Q ( x ) Q(x) Q(x) 即可一次多项式乘法解决。
设 F R ( x ) F_R(x) FR(x) 表示 F ( x ) F(x) F(x) 翻转系数形成的多项式,即 F R ( x ) = x deg f F ( 1 x ) F_R(x)=x^{\deg f}F(\frac{1}{x}) FR(x)=xdegfF(x1)。而 deg F R ( x ) = F ( x ) \deg F_R(x)=F(x) degFR(x)=F(x)。
F ( x ) = Q ( x ) × G ( x ) + R ( x ) ⇒ x n F ( 1 x ) = ( x m G ( 1 x ) ( x n − m Q ( 1 x ) ) + R ( x ) ⇒ F R ( x ) = G R ( x ) Q R ( x ) + x n − m + 1 R R ( x ) F(x) = Q(x) \times G(x) + R(x)\Rightarrow x^nF(\frac{1}{x})=(x^mG(\frac{1}{x})(x^{n-m}Q(\frac{1}{x}))+R(x)\Rightarrow F_R(x)=G_R(x)Q_R(x)+x^{n-m+1}R_R(x) F(x)=Q(x)×G(x)+R(x)⇒xnF(x1)=(xmG(x1)(xn−mQ(x1))+R(x)⇒FR(x)=GR(x)QR(x)+xn−m+1RR(x)。
我们消去 R R ( x ) R_R(x) RR(x) 的影响: F R ( x ) ≡ G R ( x ) Q R ( x ) ( m o d x n − m + 1 ) F_R(x)\equiv G_R(x)Q_R(x)\pmod {x^{n-m+1}} FR(x)≡GR(x)QR(x)(modxn−m+1)。
于是多项式求逆即可,注意细节。复杂度 O ( n log n ) O(n\log n) O(nlogn)。
#include<bits/stdc++.h>
#define fr(x) freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);
#define LL long long
using namespace std;
const int mod=998244353,N=4e5+5;
int n,m,a[N],b[N],c[N],d[N],e[N],w[N],mmax;
inline int rd()
{
int x=0,zf=1;
char ch=getchar();
while(ch<'0'||ch>'9') (ch=='-')and(zf=-1),ch=getchar();
while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
return x*zf;
}
inline void wr(int x)
{
if(x==0) return putchar('0'),putchar(' '),void();
int num[35],len=0;
while(x) num[++len]=x%10,x/=10;
for(int i=len;i>=1;i--) putchar(num[i]+'0');
putchar(' ');
}
inline int bger(int x){return x|=x>>1,x|=x>>2,x|=x>>4,x|=x>>8,x|=x>>16,x+1;}
inline int md(int x){return x>=mod?x-mod:x;}
inline int ksm(int x,int p){int s=1;for(;p;(p&1)&&(s=1ll*s*x%mod),x=1ll*x*x%mod,p>>=1);return s;}
inline void init(int mmax)
{
for(int i=1,j,k;i<mmax;i<<=1)
for(w[j=i]=1,k=ksm(3,(mod-1)/(i<<1)),j++;j<(i<<1);j++)
w[j]=1ll*w[j-1]*k%mod;
}
inline void DNT(int *a,int mmax)
{
for(int i,j,k=mmax>>1,L,*W,*x,*y,z;k;k>>=1)
for(L=k<<1,i=0;i<mmax;i+=L)
for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
*y=1ll*(*x+mod-(z=*y))* *W%mod,*x=md(*x+z);
}
inline void IDNT(int *a,int mmax)
{
for(int i,j,k=1,L,*W,*x,*y,z;k<mmax;k<<=1)
for(L=k<<1,i=0;i<mmax;i+=L)
for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
z=1ll* *W* *y%mod,*y=md(*x+mod-z),*x=md(*x+z);
reverse(a+1,a+mmax);
for(int inv=ksm(mmax,mod-2),i=0;i<mmax;i++) a[i]=1ll*a[i]*inv%mod;
}
inline void NTT(int *a,int *b,int n,int m)
{
mmax=bger(n+m);init(mmax);
DNT(a,mmax);DNT(b,mmax);
for(int i=0;i<mmax;i++) a[i]=1ll*a[i]*b[i]%mod;
IDNT(a,mmax);
}
void INV(int num,int *a,int *b)
{
if(num==1) return b[0]=ksm(a[0],mod-2),void();
INV((num+1)>>1,a,b);
int mmax=bger(num<<1);init(mmax);
static int c[N];
for(int i=0;i<num;i++) c[i]=a[i];for(int i=num;i<mmax;i++) c[i]=0;
DNT(c,mmax);DNT(b,mmax);
for(int i=0;i<mmax;i++) b[i]=1ll*(2-1ll*c[i]*b[i]%mod+mod)%mod*b[i]%mod;
IDNT(b,mmax);
for(int i=num;i<mmax;i++) b[i]=0;
}
inline void div(int *a,int *b,int n,int m,int *d)//余式是a
{
reverse(a,a+1+n);reverse(b,b+1+m);INV(n-m+1,b,c);
int mmax=bger(2*n-m+1);static int e[N];init(mmax);
for(int i=0;i<=n;i++) d[i]=a[i];
DNT(d,mmax);DNT(c,mmax);
for(int i=0;i<mmax;i++) d[i]=1ll*d[i]*c[i]%mod;
IDNT(d,mmax);reverse(d,d+1+n-m);for(int i=n-m+1;i<=n;i++) d[i]=0;
for(int i=0;i<n-m+1;i++) e[i]=d[i];
mmax=bger(n<<1);init(mmax);
reverse(a,a+1+n);reverse(b,b+1+m);
DNT(b,mmax);DNT(e,mmax);
for(int i=0;i<mmax;i++) b[i]=1ll*b[i]*e[i]%mod;
IDNT(b,mmax);
for(int i=0;i<m;i++) a[i]=(mod+a[i]-b[i])%mod;
}
signed main()
{
n=rd();m=rd();
for(int i=0;i<=n;i++) a[i]=rd();
for(int i=0;i<=m;i++) b[i]=rd();
div(a,b,n,m,d);
for(int i=0;i<n-m+1;i++) wr(d[i]);puts("");
for(int i=0;i<m;i++) wr(a[i]);
return 0;
}
多项式 ln/exp
多项式 ln
给出 n − 1 n-1 n−1 次多项式 A ( x ) A(x) A(x),求一个 m o d x n \bmod{\:x^n} modxn 下的多项式 B ( x ) B(x) B(x),满足 B ( x ) ≡ ln A ( x ) B(x) \equiv \ln A(x) B(x)≡lnA(x)
在 mod 998244353 \text{mod } 998244353 mod 998244353 下进行,且 a i ∈ [ 0 , 998244353 ) ∩ Z a_i \in [0, 998244353) \cap \mathbb{Z} ai∈[0,998244353)∩Z
$\texttt{code}$注意到 ( ln f ( x ) ) ′ = f ′ ( x ) f ( x ) (\ln f(x))'=\dfrac{f'(x)}{f(x)} (lnf(x))′=f(x)f′(x),于是 ln ( A ( x ) ) = ∫ A ′ ( x ) A ( x ) d x \ln(A(x))=\displaystyle\int \dfrac{A'(x)}{A(x)} {\rm d}x ln(A(x))=∫A(x)A′(x)dx。
只需一次求导,一次求逆,一次积分即可。复杂度 O ( n log n ) O(n\log n) O(nlogn)。
//洛谷 P4725
//https://www.luogu.com.cn/problem/P4725
#include<bits/stdc++.h>
#define fr(x) freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);
using namespace std;
const int mod=998244353,N=4e5+5;
int n,a[N],w[N],mmax;
inline int rd()
{
int x=0,zf=1;
char ch=getchar();
while(ch<'0'||ch>'9') (ch=='-')and(zf=-1),ch=getchar();
while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
return x*zf;
}
inline void wr(int x)
{
if(x==0) return putchar('0'),putchar(' '),void();
int num[35],len=0;
while(x) num[++len]=x%10,x/=10;
for(int i=len;i>=1;i--) putchar(num[i]+'0');
putchar(' ');
}
inline int bger(int x){return x|=x>>1,x|=x>>2,x|=x>>4,x|=x>>8,x|=x>>16,x+1;}
inline int md(int x){return x>=mod?x-mod:x;}
inline int ksm(int x,int p){int s=1;for(;p;(p&1)&&(s=1ll*s*x%mod),x=1ll*x*x%mod,p>>=1);return s;}
inline void dao(int *a,int n){for(int i=1;i<n;i++) a[i-1]=1ll*i*a[i]%mod;a[n-1]=0;}
inline void ji(int *a,int n){for(int i=n-1;i>=1;i--) a[i]=1ll*ksm(i,mod-2)*a[i-1]%mod;a[0]=0;}
inline void init(int mmax)
{
for(int i=1,j,k;i<mmax;i<<=1)
for(w[j=i]=1,k=ksm(3,(mod-1)/(i<<1)),j++;j<(i<<1);j++)
w[j]=1ll*w[j-1]*k%mod;
}
inline void DNT(int *a,int mmax)
{
for(int i,j,k=mmax>>1,L,*W,*x,*y,z;k;k>>=1)
for(L=k<<1,i=0;i<mmax;i+=L)
for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
*y=1ll*(*x+mod-(z=*y))* *W%mod,*x=md(*x+z);
}
inline void IDNT(int *a,int mmax)
{
for(int i,j,k=1,L,*W,*x,*y,z;k<mmax;k<<=1)
for(L=k<<1,i=0;i<mmax;i+=L)
for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
z=1ll* *W* *y%mod,*y=md(*x+mod-z),*x=md(*x+z);
reverse(a+1,a+mmax);
for(int inv=ksm(mmax,mod-2),i=0;i<mmax;i++) a[i]=1ll*a[i]*inv%mod;
}
inline void NTT(int *a,int *b,int n,int m)
{
mmax=bger(n+m);init(mmax);
DNT(a,mmax);DNT(b,mmax);
for(int i=0;i<mmax;i++) a[i]=1ll*a[i]*b[i]%mod;
IDNT(a,mmax);
}
void INV(int num,int *a,int *b)
{
if(num==1) return b[0]=ksm(a[0],mod-2),void();
INV((num+1)>>1,a,b);
int mmax=bger(num<<1);init(mmax);
static int c[N];
for(int i=0;i<num;i++) c[i]=a[i];for(int i=num;i<mmax;i++) c[i]=0;
DNT(c,mmax);DNT(b,mmax);
for(int i=0;i<mmax;i++) b[i]=1ll*(2-1ll*c[i]*b[i]%mod+mod)%mod*b[i]%mod;
IDNT(b,mmax);
for(int i=num;i<mmax;i++) b[i]=0;
}
inline void Ln(int *a,int n){static int b[N];for(int i=0;i<bger(n<<1);i++) b[i]=0;INV(n,a,b);dao(a,n);NTT(a,b,n,n);ji(a,n);for(int i=n;i<bger(n<<1);i++) a[i]=0;}
int main()
{
n=rd();for(int i=0;i<n;i++) a[i]=rd();
Ln(a,n);for(int i=0;i<n;i++) wr(a[i]);
return 0;
}
多项式 exp
给出 n − 1 n-1 n−1 次多项式 A ( x ) A(x) A(x),求一个 m o d x n \bmod{\:x^n} modxn 下的多项式 B ( x ) B(x) B(x),满足 B ( x ) ≡ e A ( x ) B(x) \equiv \text e^{A(x)} B(x)≡eA(x)。系数对 998244353 998244353 998244353 取模。
$\texttt{code}$考虑牛顿迭代。
构造 g ( X ) = ln ( X ) − A ( x ) g(X)=\ln(X)-A(x) g(X)=ln(X)−A(x),则 g ( B ( x ) ) ≡ 0 ( m o d x n ) g(B(x))\equiv 0\pmod {x^n} g(B(x))≡0(modxn)。
设 B 0 ( x ) B_0(x) B0(x) 满足 g ( B 0 ( x ) ) ≡ 0 ( m o d x ⌈ x / 2 ⌉ ) g(B_0(x))\equiv 0\pmod {x^{\lceil x/2\rceil}} g(B0(x))≡0(modx⌈x/2⌉)
则 B ( x ) ≡ B 0 ( x ) − g ( B 0 ( x ) ) g ′ ( B 0 ( x ) ) = B 0 ( x ) − ln ( B 0 ( x ) ) − A ( x ) 1 B 0 ( x ) = B 0 ( x ) ( 1 − ln ( B 0 ( x ) ) + A ( x ) ) ( m o d x n ) B(x)\equiv B_0(x)-\dfrac{g(B_0(x))}{g'(B_0(x))}=B_0(x)-\dfrac{\ln(B_0(x))-A(x)}{\dfrac{1}{B_0(x)}}=B_0(x)(1-\ln(B_0(x))+ A(x))\pmod {x^n} B(x)≡B0(x)−g′(B0(x))g(B0(x))=B0(x)−B0(x)1ln(B0(x))−A(x)=B0(x)(1−ln(B0(x))+A(x))(modxn)
这其中需要一次多项式乘法,一次多项式 ln \ln ln,于是复杂度 T ( n ) = T ( n / 2 ) + O ( n log n ) = O ( n log n ) T(n)=T(n/2)+O(n\log n)=O(n\log n) T(n)=T(n/2)+O(nlogn)=O(nlogn)。
//落谷 P4726
//https://www.luogu.com.cn/problem/P4726
#include<bits/stdc++.h>
#define fr(x) freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);
using namespace std;
const int mod=998244353,N=4e5+5;
int n,a[N],b[N],w[N],mmax;
inline int rd()
{
int x=0,zf=1;
char ch=getchar();
while(ch<'0'||ch>'9') (ch=='-')and(zf=-1),ch=getchar();
while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
return x*zf;
}
inline void wr(int x)
{
if(x==0) return putchar('0'),putchar(' '),void();
int num[35],len=0;
while(x) num[++len]=x%10,x/=10;
for(int i=len;i>=1;i--) putchar(num[i]+'0');
putchar(' ');
}
inline int bger(int x){return x|=x>>1,x|=x>>2,x|=x>>4,x|=x>>8,x|=x>>16,x+1;}
inline int md(int x){return x>=mod?x-mod:x;}
inline int ksm(int x,int p){int s=1;for(;p;(p&1)&&(s=1ll*s*x%mod),x=1ll*x*x%mod,p>>=1);return s;}
inline void dao(int *a,int n){for(int i=1;i<n;i++) a[i-1]=1ll*i*a[i]%mod;a[n-1]=0;}
inline void ji(int *a,int n){for(int i=n-1;i>=1;i--) a[i]=1ll*ksm(i,mod-2)*a[i-1]%mod;a[0]=0;}
inline void init(int mmax)
{
for(int i=1,j,k;i<mmax;i<<=1)
for(w[j=i]=1,k=ksm(3,(mod-1)/(i<<1)),j++;j<(i<<1);j++)
w[j]=1ll*w[j-1]*k%mod;
}
inline void DNT(int *a,int mmax)
{
for(int i,j,k=mmax>>1,L,*W,*x,*y,z;k;k>>=1)
for(L=k<<1,i=0;i<mmax;i+=L)
for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
*y=1ll*(*x+mod-(z=*y))* *W%mod,*x=md(*x+z);
}
inline void IDNT(int *a,int mmax)
{
for(int i,j,k=1,L,*W,*x,*y,z;k<mmax;k<<=1)
for(L=k<<1,i=0;i<mmax;i+=L)
for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
z=1ll* *W* *y%mod,*y=md(*x+mod-z),*x=md(*x+z);
reverse(a+1,a+mmax);
for(int inv=ksm(mmax,mod-2),i=0;i<mmax;i++) a[i]=1ll*a[i]*inv%mod;
}
inline void NTT(int *a,int *b,int n,int m)
{
mmax=bger(n+m);init(mmax);
DNT(a,mmax);DNT(b,mmax);
for(int i=0;i<mmax;i++) a[i]=1ll*a[i]*b[i]%mod;
IDNT(a,mmax);
}
void INV(int num,int *a,int *b)
{
if(num==1) return b[0]=ksm(a[0],mod-2),void();
INV((num+1)>>1,a,b);
int mmax=bger(num<<1);init(mmax);
static int c[N];
for(int i=0;i<num;i++) c[i]=a[i];for(int i=num;i<mmax;i++) c[i]=0;
DNT(c,mmax);DNT(b,mmax);
for(int i=0;i<mmax;i++) b[i]=1ll*(2-1ll*c[i]*b[i]%mod+mod)%mod*b[i]%mod;
IDNT(b,mmax);
for(int i=num;i<mmax;i++) b[i]=0;
}
inline void Ln(int *a,int n){static int b[N];for(int i=0;i<bger(n<<1);i++) b[i]=0;INV(n,a,b);dao(a,n);NTT(a,b,n,n);ji(a,n);for(int i=n;i<bger(n<<1);i++) a[i]=0;}
inline void Exp(int *a,int *b,int n)
{
if(n==1) return b[0]=1,void();
Exp(a,b,(n+1)>>1);static int c[N];for(int i=0;i<bger(n<<1);i++) c[i]=0;
for(int i=0;i<n;i++) c[i]=b[i];Ln(c,n);
for(int i=0;i<n;i++) c[i]=md(mod-c[i]+a[i]);c[0]=md(c[0]+1);
NTT(b,c,n,n);for(int i=n;i<bger(n<<1);i++) b[i]=0;
}
int main()
{
n=rd();for(int i=0;i<n;i++) a[i]=rd();
Exp(a,b,n);for(int i=0;i<n;i++) wr(b[i]);
return 0;
}
多项式快速幂/加强版
给定一个 n − 1 n-1 n−1 次多项式 A ( x ) A(x) A(x),求一个在 m o d x n \bmod\ x^n mod xn 意义下的多项式 B ( x ) B(x) B(x),使得 B ( x ) ≡ ( A ( x ) ) k ( m o d x n ) B(x) \equiv (A(x))^k \ (\bmod\ x^n) B(x)≡(A(x))k (mod xn)。
多项式的系数在 m o d 998244353 \bmod\ 998244353 mod 998244353 的意义下进行运算。
普通版保证 a 0 = 1 a_0=1 a0=1,加强版则不保证。
$\texttt{code}$注意到 ( A ( x ) ) k = exp ( k ln ( A ( x ) ) ) (A(x))^k=\exp(k\ln(A(x))) (A(x))k=exp(kln(A(x))),于是一次 exp 一次 ln 即可,复杂度 O ( n log n ) O(n\log n) O(nlogn)。
//洛谷 P4726
//https://www.luogu.com.cn/problem/P4726
#include<bits/stdc++.h>
#define fr(x) freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);
using namespace std;
const int mod=998244353,N=4e6+5;
int n,m,a[N],b[N],w[N],mmax;
inline int rd()
{
int x=0,zf=1;
char ch=getchar();
while(ch<'0'||ch>'9') (ch=='-')and(zf=-1),ch=getchar();
while(ch>='0'&&ch<='9') x=(10ll*x%mod+ch-'0')%mod,ch=getchar();
return x*zf;
}
inline void wr(int x)
{
if(x==0) return putchar('0'),putchar(' '),void();
int num[35],len=0;
while(x) num[++len]=x%10,x/=10;
for(int i=len;i>=1;i--) putchar(num[i]+'0');
putchar(' ');
}
inline int bger(int x){return x|=x>>1,x|=x>>2,x|=x>>4,x|=x>>8,x|=x>>16,x+1;}
inline int md(int x){return x>=mod?x-mod:x;}
inline int ksm(int x,int p){int s=1;for(;p;(p&1)&&(s=1ll*s*x%mod),x=1ll*x*x%mod,p>>=1);return s;}
inline void dao(int *a,int n){for(int i=1;i<n;i++) a[i-1]=1ll*i*a[i]%mod;a[n-1]=0;}
inline void ji(int *a,int n){for(int i=n-1;i>=1;i--) a[i]=1ll*ksm(i,mod-2)*a[i-1]%mod;a[0]=0;}
inline void init(int mmax)
{
for(int i=1,j,k;i<mmax;i<<=1)
for(w[j=i]=1,k=ksm(3,(mod-1)/(i<<1)),j++;j<(i<<1);j++)
w[j]=1ll*w[j-1]*k%mod;
}
inline void DNT(int *a,int mmax)
{
for(int i,j,k=mmax>>1,L,*W,*x,*y,z;k;k>>=1)
for(L=k<<1,i=0;i<mmax;i+=L)
for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
*y=1ll*(*x+mod-(z=*y))* *W%mod,*x=md(*x+z);
}
inline void IDNT(int *a,int mmax)
{
for(int i,j,k=1,L,*W,*x,*y,z;k<mmax;k<<=1)
for(L=k<<1,i=0;i<mmax;i+=L)
for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
z=1ll* *W* *y%mod,*y=md(*x+mod-z),*x=md(*x+z);
reverse(a+1,a+mmax);
for(int inv=ksm(mmax,mod-2),i=0;i<mmax;i++) a[i]=1ll*a[i]*inv%mod;
}
inline void NTT(int *a,int *b,int n,int m)
{
mmax=bger(n+m);init(mmax);
DNT(a,mmax);DNT(b,mmax);
for(int i=0;i<mmax;i++) a[i]=1ll*a[i]*b[i]%mod;
IDNT(a,mmax);
}
void INV(int num,int *a,int *b)
{
if(num==1) return b[0]=ksm(a[0],mod-2),void();
INV((num+1)>>1,a,b);
int mmax=bger(num<<1);init(mmax);
static int c[N];
for(int i=0;i<num;i++) c[i]=a[i];for(int i=num;i<mmax;i++) c[i]=0;
DNT(c,mmax);DNT(b,mmax);
for(int i=0;i<mmax;i++) b[i]=1ll*(2-1ll*c[i]*b[i]%mod+mod)%mod*b[i]%mod;
IDNT(b,mmax);
for(int i=num;i<mmax;i++) b[i]=0;
}
inline void Ln(int *a,int n){static int b[N];for(int i=0;i<bger(n<<1);i++) b[i]=0;INV(n,a,b);dao(a,n);NTT(a,b,n,n);ji(a,n);for(int i=n;i<bger(n<<1);i++) a[i]=0;}
inline void Exp(int *a,int *b,int n)
{
if(n==1) return b[0]=1,void();
Exp(a,b,(n+1)>>1);static int c[N];for(int i=0;i<bger(n<<1);i++) c[i]=0;
for(int i=0;i<n;i++) c[i]=b[i];Ln(c,n);
for(int i=0;i<n;i++) c[i]=md(mod-c[i]+a[i]);c[0]=md(c[0]+1);
NTT(b,c,n,n);for(int i=n;i<bger(n<<1);i++) b[i]=0;
}
int main()
{
n=rd();m=rd();for(int i=0;i<n;i++) a[i]=rd();
Ln(a,n);for(int i=0;i<n;i++) a[i]=1ll*a[i]*m%mod;
Exp(a,b,n);for(int i=0;i<n;i++) wr(b[i]);
return 0;
}
$\texttt{code}$加强版只需注意到:
由于费马小定理,幂次可以对 998244352 998244352 998244352 取模。
A ( x ) k = c k ( A ( x ) c ) k A(x)^k=c^k\left(\dfrac{A(x)}{c}\right)^k A(x)k=ck(cA(x))k,而后求后面的幂次。
设 A ( x ) = ∑ i = 0 n − 1 a i x i A(x)=\sum\limits_{i=0}^{n-1} a_ix^i A(x)=i=0∑n−1aixi,第一个非零次为 x t x^t xt,则 A ( x ) k = x t k ( ∑ i = t n − 1 a i x i − t ) k A(x)^k=x^{tk}\left(\sum\limits_{i=t}^{n-1}a_ix^{i-t}\right)^k A(x)k=xtk(i=t∑n−1aixi−t)k,而后求后面的幂次。
于是就能转化为 a 0 = 1 a_0=1 a0=1 的情况啦。复杂度不变。当然加强版基本 useless 啦!
//落谷 P4726
//https://www.luogu.com.cn/problem/P4726
#include<bits/stdc++.h>
#define fr(x) freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);
using namespace std;
const int mod=998244353,N=4e5+5;
int n,m,a[N],b[N],w[N],mmax;
char str[N];
inline int rd()
{
int x=0,zf=1;
char ch=getchar();
while(ch<'0'||ch>'9') (ch=='-')and(zf=-1),ch=getchar();
while(ch>='0'&&ch<='9') x=(10*x+ch-'0'),ch=getchar();
return x*zf;
}
inline void wr(int x)
{
if(x==0) return putchar('0'),putchar(' '),void();
int num[35],len=0;
while(x) num[++len]=x%10,x/=10;
for(int i=len;i>=1;i--) putchar(num[i]+'0');
putchar(' ');
}
inline int bger(int x){return x|=x>>1,x|=x>>2,x|=x>>4,x|=x>>8,x|=x>>16,x+1;}
inline int md(int x){return x>=mod?x-mod:x;}
inline int ksm(int x,int p){int s=1;for(;p;(p&1)&&(s=1ll*s*x%mod),x=1ll*x*x%mod,p>>=1);return s;}
inline void dao(int *a,int n){for(int i=1;i<n;i++) a[i-1]=1ll*i*a[i]%mod;a[n-1]=0;}
inline void ji(int *a,int n){for(int i=n-1;i>=1;i--) a[i]=1ll*ksm(i,mod-2)*a[i-1]%mod;a[0]=0;}
inline void init(int mmax)
{
for(int i=1,j,k;i<mmax;i<<=1)
for(w[j=i]=1,k=ksm(3,(mod-1)/(i<<1)),j++;j<(i<<1);j++)
w[j]=1ll*w[j-1]*k%mod;
}
inline void DNT(int *a,int mmax)
{
for(int i,j,k=mmax>>1,L,*W,*x,*y,z;k;k>>=1)
for(L=k<<1,i=0;i<mmax;i+=L)
for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
*y=1ll*(*x+mod-(z=*y))* *W%mod,*x=md(*x+z);
}
inline void IDNT(int *a,int mmax)
{
for(int i,j,k=1,L,*W,*x,*y,z;k<mmax;k<<=1)
for(L=k<<1,i=0;i<mmax;i+=L)
for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
z=1ll* *W* *y%mod,*y=md(*x+mod-z),*x=md(*x+z);
reverse(a+1,a+mmax);
for(int inv=ksm(mmax,mod-2),i=0;i<mmax;i++) a[i]=1ll*a[i]*inv%mod;
}
inline void NTT(int *a,int *b,int n,int m)
{
mmax=bger(n+m);init(mmax);
DNT(a,mmax);DNT(b,mmax);
for(int i=0;i<mmax;i++) a[i]=1ll*a[i]*b[i]%mod;
IDNT(a,mmax);
}
void INV(int num,int *a,int *b)
{
if(num==1) return b[0]=ksm(a[0],mod-2),void();
INV((num+1)>>1,a,b);
int mmax=bger(num<<1);init(mmax);
static int c[N];
for(int i=0;i<num;i++) c[i]=a[i];for(int i=num;i<mmax;i++) c[i]=0;
DNT(c,mmax);DNT(b,mmax);
for(int i=0;i<mmax;i++) b[i]=1ll*(2-1ll*c[i]*b[i]%mod+mod)%mod*b[i]%mod;
IDNT(b,mmax);
for(int i=num;i<mmax;i++) b[i]=0;
}
inline void Ln(int *a,int n){static int b[N];for(int i=0;i<bger(n<<1);i++) b[i]=0;INV(n,a,b);dao(a,n);NTT(a,b,n,n);ji(a,n);for(int i=n;i<bger(n<<1);i++) a[i]=0;}
inline void Exp(int *a,int *b,int n)
{
if(n==1) return b[0]=1,void();
Exp(a,b,(n+1)>>1);static int c[N];for(int i=0;i<bger(n<<1);i++) c[i]=0;
for(int i=0;i<n;i++) c[i]=b[i];Ln(c,n);
for(int i=0;i<n;i++) c[i]=md(mod-c[i]+a[i]);c[0]=md(c[0]+1);
NTT(b,c,n,n);for(int i=n;i<bger(n<<1);i++) b[i]=0;
}
inline void ksm(int *a,char *str,int n)
{
int t=n,t1,t2,t3,L=strlen(str),X=0,XX=0;
static int A[N],b[N];for(int i=0;i<n;i++) A[i]=b[i]=0;
for(int i=0;i<L;i++) m=(10ll*m%mod+str[i]-'0')%mod;
for(int i=0;i<L;i++) X=(10ll*X%(mod-1)+str[i]-'0')%(mod-1);
for(int i=0;i<min(L,7);i++) XX=(10*XX+str[i]-'0');
if(a[0]==0&&XX>=n){for(int i=0;i<n;i++) a[i]=0;return;}
for(int i=0;i<n;i++) if(a[i]!=0){t=i;break;}
for(int i=t;i<n;i++) A[i-t]=a[i];t1=ksm(A[0],mod-2),t2=ksm(A[0],X);
for(int i=0;i<n-t;i++) A[i]=1ll*A[i]*t1%mod;
Ln(A,n-t);for(int i=0;i<n-t;i++) A[i]=1ll*A[i]*m%mod;t3=min(1ll*t*m,1ll*n);
Exp(A,b,n-t);for(int i=0;i<t3;i++) A[i]=0;for(int i=t3;i<n;i++) A[i]=1ll*b[i-t3]*t2%mod;
for(int i=0;i<n;i++) a[i]=A[i];
}
int main()
{
n=rd();scanf("%s",str);
for(int i=0;i<n;i++) a[i]=rd();
ksm(a,str,n);for(int i=0;i<n;i++) wr(a[i]);
return 0;
}
多项式开根/加强版
给定一个 n − 1 n-1 n−1 次多项式 A ( x ) A(x) A(x) ,求一个在 m o d x n \bmod\ {x^n} mod xn 意义下的多项式 B ( x ) B(x) B(x) ,使得 B 2 ( x ) ≡ A ( x ) ( m o d x n ) B^2(x)\equiv A(x)\pmod {x^n} B2(x)≡A(x)(modxn)。
多项式的系数在 m o d 998244353 \bmod\ 998244353 mod 998244353 的意义下进行运算。
普通版保证 a 0 = 1 a_0=1 a0=1,加强版仅保证 a 0 a_0 a0 是 m o d 998244353 \bmod\ 998244353 mod 998244353 的二次剩余。
加强版 $\texttt{code}$直接讲加强版,设我们求出 a 0 a_0 a0 的二次剩余为 r r r。
考虑牛顿迭代。 g ( X ) = X 2 − A ( x ) g(X)=X^2-A(x) g(X)=X2−A(x),则 g ( B ( x ) ) ≡ 0 ( m o d x n ) g(B(x))\equiv 0\pmod {x^n} g(B(x))≡0(modxn)。
设 B 0 ( x ) B_0(x) B0(x) 满足 g ( B 0 ( x ) ) ≡ 0 ( m o d x ⌈ x / 2 ⌉ ) g(B_0(x))\equiv 0\pmod {x^{\lceil x/2\rceil}} g(B0(x))≡0(modx⌈x/2⌉)
则 B ( x ) ≡ B 0 ( x ) − g ( B 0 ( x ) ) g ′ ( B 0 ( x ) ) = B 0 ( x ) − B 0 ( x ) 2 − A ( x ) 2 B 0 ( x ) = B 0 ( x ) 2 + A ( x ) 2 B 0 ( x ) ( m o d x n ) B(x)\equiv B_0(x)-\dfrac{g(B_0(x))}{g'(B_0(x))}=B_0(x)-\dfrac{B_0(x)^2-A(x)}{2B_0(x)}=\dfrac{B_0(x)}{2}+\dfrac{A(x)}{2B_0(x)}\pmod {x^n} B(x)≡B0(x)−g′(B0(x))g(B0(x))=B0(x)−2B0(x)B0(x)2−A(x)=2B0(x)+2B0(x)A(x)(modxn)。
于是这里一次多项式求逆即可,边界时 n = 1 n=1 n=1, B ( x ) = r B(x)=r B(x)=r。
复杂度 T ( n ) = T ( n / 2 ) + O ( n log n ) = O ( n log n ) T(n)=T(n/2)+O(n\log n)=O(n\log n) T(n)=T(n/2)+O(nlogn)=O(nlogn)。
//洛谷 P5277
//https://www.luogu.com.cn/problem/P5277
#include<bits/stdc++.h>
#define fr(x) freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);
#define LL long long
using namespace std;
const int mod=998244353,N=4e6+5;
int n,m,a[N],b[N],c,w[N],jc[N],inv[N],mmax;
inline int rd()
{
int x=0,zf=1;char ch=getchar();
while(ch<'0'||ch>'9') (ch=='-')and(zf=-1),ch=getchar();
while(ch>='0'&&ch<='9') x=((x<<3)+(x<<1)+ch-'0')%mod,ch=getchar();
return x*zf;
}
inline void wr(int x)
{
if(x==0) return putchar('0'),putchar(' '),void();
int num[35],len=0;
while(x) num[++len]=x%10,x/=10;
for(int i=len;i>=1;i--) putchar(num[i]+'0');
putchar(' ');
}
inline int bger(int x){return x|=x>>1,x|=x>>2,x|=x>>4,x|=x>>8,x|=x>>16,x+1;}
inline int md(int x){return x>=mod?x-mod:x;}
inline int yy(int x){if(x&1) return (x+mod)>>1;return x>>1;}
inline int ksm(int x,int p){int s=1;for(;p;(p&1)&&(s=1ll*s*x%mod),x=1ll*x*x%mod,p>>=1);return s;}
inline void dao(int *a,int n){for(int i=1;i<n;i++) a[i-1]=1ll*i*a[i]%mod;a[n-1]=0;}
inline void ji(int *a,int n){for(int i=n-1;i>=1;i--) a[i]=1ll*ksm(i,mod-2)*a[i-1]%mod;a[0]=0;}
namespace er
{
int n,w;
struct comp
{
int x,y;
inline comp operator *(comp X){return {(1ll*x*X.x+1ll*y*X.y%mod*w)%mod,(1ll*x*X.y+1ll*y*X.x)%mod};}
};
inline comp ksm1(comp x,int p){comp s={1,0};for(;p;(p&1)&&(s=s*x,1),x=x*x,p>>=1);return s;}
inline int Find()
{
int x;
while(1)
{
x=rand()%mod;w=((LL)x*x%mod-n+mod)%mod;
if(ksm(w,(mod-1)>>1)==mod-1) break;
}
return ksm1({x,1},(mod+1)>>1).x;
}
}using er::Find;
inline void init(int mmax)
{
for(int i=1,j,k;i<mmax;i<<=1)
for(w[j=i]=1,k=ksm(3,(mod-1)/(i<<1)),j++;j<(i<<1);j++)
w[j]=1ll*w[j-1]*k%mod;
}
inline void DNT(int *a,int mmax)
{
for(int i,j,k=mmax>>1,L,*W,*x,*y,z;k;k>>=1)
for(L=k<<1,i=0;i<mmax;i+=L)
for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
*y=1ll*(*x+mod-(z=*y))* *W%mod,*x=md(*x+z);
}
inline void IDNT(int *a,int mmax)
{
for(int i,j,k=1,L,*W,*x,*y,z;k<mmax;k<<=1)
for(L=k<<1,i=0;i<mmax;i+=L)
for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
z=1ll* *W* *y%mod,*y=md(*x+mod-z),*x=md(*x+z);
reverse(a+1,a+mmax);
for(int inv=ksm(mmax,mod-2),i=0;i<mmax;i++) a[i]=1ll*a[i]*inv%mod;
}
inline void NTT(int *a,int *b,int n,int m)
{
mmax=bger(n+m);init(mmax);
DNT(a,mmax);DNT(b,mmax);
for(int i=0;i<mmax;i++) a[i]=1ll*a[i]*b[i]%mod;
IDNT(a,mmax);
}
void INV(int num,int *a,int *b)
{
if(num==1) return b[0]=ksm(a[0],mod-2),void();INV((num+1)>>1,a,b);
int mmax=bger(num<<1);init(mmax);static int c[N];
for(int i=0;i<num;i++) c[i]=a[i];for(int i=num;i<mmax;i++) c[i]=0;DNT(c,mmax);DNT(b,mmax);
for(int i=0;i<mmax;i++) b[i]=1ll*(2-1ll*c[i]*b[i]%mod+mod)%mod*b[i]%mod;IDNT(b,mmax);
for(int i=num;i<mmax;i++) b[i]=0;
}
void SQ(int n,int *a,int *b)
{
if(n==1) return b[0]=c,void();SQ((n+1)>>1,a,b);static int c[N],d[N];
for(int i=0;i<n;i++) c[i]=b[i];INV(n,c,d);for(int i=n;i<mmax;i++) c[i]=d[i]=0;
for(int i=0;i<n;i++) c[i]=a[i];NTT(c,d,n-1,n-1);
for(int i=0;i<n;i++) b[i]=yy(md(c[i]+b[i]));
for(int i=0;i<mmax;i++) c[i]=d[i]=0;
}
int main()
{
srand(time(NULL));n=rd();for(int i=0;i<n;i++) a[i]=rd();
er::n=a[0];c=Find();SQ(n,a,b);
for(int i=0;i<n;i++) wr((b[0]<=mod-b[0])?b[i]:mod-b[i]);
return 0;
}
挑战多项式
按照惯例需要给出挑战多项式的代码作为小全家桶。代码。
常系数齐次线性递推
求一个满足 k k k 阶齐次线性递推数列 a i {a_i} ai 的第 n n n 项,即给定 f 1 , 2 … k , a 0 , 1 … k − 1 f_{1,2\dots k},a_{0,1\dots k-1} f1,2…k,a0,1…k−1,以及递推式: a n = ∑ i = 1 k f i × a n − i a_n=\sum\limits_{i=1}^{k}f_i \times a_{n-i} an=i=1∑kfi×an−i,求 a n a_n an
k ≤ 3.2 × 1 0 4 , n ≤ 1 0 9 k\le 3.2\times10^4,n\le 10^9 k≤3.2×104,n≤109.
前置知识-- Bostan-mori 算法
$\texttt{function}$我们有一个生成函数 H ( x ) = F ( x ) G ( x ) H(x)=\dfrac{F(x)}{G(x)} H(x)=G(x)F(x) 。其中 deg ( F ) , deg ( G ) ≤ k ≤ 32000 \deg(F),\deg(G)\le k\le32000 deg(F),deg(G)≤k≤32000。要求其第 n ≤ 1 0 9 n\le10^9 n≤109 项。即求 [ x n ] F ( x ) G ( x ) [x^n]\dfrac{F(x)}{G(x)} [xn]G(x)F(x)
。推柿子: [ x n ] F ( x ) G ( x ) = [ x n ] F ( x ) G ( − x ) G ( x ) G ( − x ) = [ x n ] F 0 ( x 2 ) + x F 1 ( x 2 ) G 0 ( x 2 ) = { [ x n / 2 ] F 0 ( x ) G 0 ( x ) ( 2 ∣ n ) [ x n / 2 ] F 1 ( x ) G 0 ( x ) ( 2 ∤ n ) [x^n]\dfrac{F(x)}{G(x)}=[x^n]\dfrac{F(x)G(-x)}{G(x)G(-x)}=[x^n]\dfrac{F_0(x^2)+xF_1(x^2)}{G_0(x^2)}=\begin{cases}[x^{n/2}]\dfrac{F_0(x)}{G_0(x)}(2\mid n)\\\\ [x^{n/2}]\dfrac{F_1(x)}{G_0(x)}(2\nmid n)\end{cases} [xn]G(x)F(x)=[xn]G(x)G(−x)F(x)G(−x)=[xn]G0(x2)F0(x2)+xF1(x2)=⎩ ⎨ ⎧[xn/2]G0(x)F0(x)(2∣n)[xn/2]G0(x)F1(x)(2∤n)
又 [ x 0 ] F ( x ) G ( x ) = [ x 0 ] F ( x ) [ x 0 ] G ( x ) [x^0]\dfrac{F(x)}{G(x)}=\dfrac{[x^0]F(x)}{[x^0]G(x)} [x0]G(x)F(x)=[x0]G(x)[x0]F(x)。于是就可以求啦!
有一个细节:就是 G ( − x ) G(-x) G(−x) 要和两个多项式分别做卷积,于是要开两个数组(下面的 c , d c,d c,d 数组)分别和两个做卷积,不能只用一个数组分别和两个做卷积。
inline int genf(int *a,int *b,int n,int m,int k)
{
static int c[N],d[N];
for(;k;k>>=1)
{
memset(c,0,sizeof(c));for(int i=0;i<=m;i++) c[i]=b[i];
for(int i=1;i<=m;i+=2) c[i]=(mod-c[i])%mod;
memset(d,0,sizeof(d));for(int i=0;i<=m;i++) d[i]=c[i];
NTT(a,c,n,m);NTT(b,d,m,m);int j=0;
for(int i=(k&1);i<=n+m;i+=2,j++) a[i/2]=a[i];n=j-1;j=0;
for(int i=0;i<=2*m;i+=2,j++) b[i/2]=b[i];m=j-1;
for(int i=n+1;i<=N-5;i++) a[i]=0;
for(int i=m+1;i<=N-5;i++) b[i]=0;
}
return 1ll*a[0]*ksm(b[0],mod-2)%mod;
}
计算递推部分
$\texttt{code}$其实内核不难发现。只是名字高端一点。
已知: a 0 , 1 , . . . , k − 1 , f 1 , 2 , . . . , k , a n = ∑ i = 1 k f i × a n − i a_{0,1,...,k-1},f_{1,2,...,k},a_n=\sum\limits_{i=1}^{k}f_i \times a_{n-i} a0,1,...,k−1,f1,2,...,k,an=i=1∑kfi×an−i。要求 a a a 的生成函数。
设 a a a 表示 ∑ i ≥ 0 a i x i \sum\limits_{i\ge 0} a_ix^i i≥0∑aixi,由递推式可得: a = a × ∑ i = 1 k f i x i + t ( x ) ( deg ( t ) ≤ k − 1 ) a=a\times\sum\limits_{i=1}^k f_ix^i +t(x)(\deg(t)\le k-1) a=a×i=1∑kfixi+t(x)(deg(t)≤k−1)。其中 a a a 是原来的数列, t ( x ) = ∑ i = 0 k − 1 a i x i − ∑ i = 0 k − 2 a i ∑ j = 1 k − 1 − i f j x i + j t(x)=\sum\limits_{i=0}^{k-1}a_ix^i-\sum\limits_{i=0}^{k-2}a_i\sum\limits_{j=1}^{k-1-i}f_jx^{i+j} t(x)=i=0∑k−1aixi−i=0∑k−2aij=1∑k−1−ifjxi+j。是一个幂次小于 k − 1 k-1 k−1 的一个多项式(很难求,于是我们不需要知道它是多少)。
于是 a = t ( x ) T ( x ) ( T ( x ) = 1 − ∑ i = 1 k f i x i ) a=\dfrac{t(x)}{T(x)}(T(x)=1-\sum\limits_{i=1}^k f_ix^i) a=T(x)t(x)(T(x)=1−i=1∑kfixi)。
由于我们知道 a a a 的前 k − 1 k-1 k−1 项 , T ( x ) T(x) T(x) 又有足足 k k k 次,于是 t ( x ) t(x) t(x) 的前 k − 1 k-1 k−1 项可以通过 A × T A\times T A×T 做一个卷积(其中 A A A 为 a a a 的前 k k k 项构成的多项式)模 x k x^k xk 唯一确定。
最后套用 Bostan-mori 算法计算这个生成函数的第 n n n 项即可,复杂度 O ( k log k log n ) O(k\log k\log n) O(klogklogn)。
#include<bits/stdc++.h>
#define LL long long
#define fr(x) freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);
using namespace std;
const int mod=998244353,N=4e5+5;
int n,k,a[N],b[N],f[N],w[N],mmax;
inline int rd()
{
int x=0,zf=1;
char ch=getchar();
while(ch<'0'||ch>'9') (ch=='-')and(zf=-1),ch=getchar();
while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
return x*zf;
}
inline void wr(int x)
{
if(x==0) return putchar('0'),putchar(' '),void();
int num[35],len=0;
while(x) num[++len]=x%10,x/=10;
for(int i=len;i>=1;i--) putchar(num[i]+'0');
putchar(' ');
}
inline int bger(int x){return x|=x>>1,x|=x>>2,x|=x>>4,x|=x>>8,x|=x>>16,x+1;}
inline int md(int x){return x>=mod?x-mod:x;}
inline int ksm(int x,int p){int s=1;for(;p;(p&1)&&(s=1ll*s*x%mod),x=1ll*x*x%mod,p>>=1);return s;}
inline void init(int mmax)
{
for(int i=1,j,k;i<mmax;i<<=1)
for(w[j=i]=1,k=ksm(3,(mod-1)/(i<<1)),j++;j<(i<<1);j++)
w[j]=1ll*w[j-1]*k%mod;
}
inline void DNT(int *a,int mmax)
{
for(int i,j,k=mmax>>1,L,*W,*x,*y,z;k;k>>=1)
for(L=k<<1,i=0;i<mmax;i+=L)
for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
*y=1ll*(*x+mod-(z=*y))* *W%mod,*x=md(*x+z);
}
inline void IDNT(int *a,int mmax)
{
for(int i,j,k=1,L,*W,*x,*y,z;k<mmax;k<<=1)
for(L=k<<1,i=0;i<mmax;i+=L)
for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
z=1ll* *W* *y%mod,*y=md(*x+mod-z),*x=md(*x+z);
reverse(a+1,a+mmax);
for(int inv=ksm(mmax,mod-2),i=0;i<mmax;i++) a[i]=1ll*a[i]*inv%mod;
}
inline void NTT(int *a,int *b,int n,int m)
{
mmax=bger(n+m);init(mmax);
DNT(a,mmax);DNT(b,mmax);
for(int i=0;i<mmax;i++) a[i]=1ll*a[i]*b[i]%mod;
IDNT(a,mmax);
}
inline int genf(int *a,int *b,int n,int m,int k)
{
static int c[N],d[N];
for(;k;k>>=1)
{
memset(c,0,sizeof(c));for(int i=0;i<=m;i++) c[i]=b[i];
for(int i=1;i<=m;i+=2) c[i]=md(mod-c[i]);
memset(d,0,sizeof(d));for(int i=0;i<=m;i++) d[i]=c[i];
NTT(a,c,n,m);NTT(b,d,m,m);int j=0;
for(int i=(k&1);i<=n+m;i+=2,j++) a[i/2]=a[i];n=j-1;j=0;
for(int i=0;i<=2*m;i+=2,j++) b[i/2]=b[i];m=j-1;
for(int i=n+1;i<=N-5;i++) a[i]=0;
for(int i=m+1;i<=N-5;i++) b[i]=0;
}
return 1ll*a[0]*ksm(b[0],mod-2)%mod;
}
inline int rec(int *f,int *a,int n,int k)
{
for(int i=1;i<=k;i++) f[i]=md(mod-f[i]);f[0]=1;
static int c[N];memset(c,0,sizeof(c));
for(int i=0;i<=k;i++) c[i]=f[i];
NTT(c,a,k,k);for(int i=k;i<=N-5;i++) c[i]=0;
return genf(c,f,k,k,n);
}
int main()
{
n=rd();k=rd();
for(int i=1;i<=k;i++) f[i]=md(rd()%mod+mod);
for(int i=0;i<k;i++) a[i]=md(rd()%mod+mod);
wr(rec(f,a,n,k));
return 0;
}