command_block
2019-01-31 13:47:34
(
(
(
(
(
(
书接上回 : 傅里叶变换(FFT)学习笔记
学完了FFT乘法之后这个坑就鸽了半年……
这里的全家桶指 {求逆+牛顿迭代+带余除法+Ln+Exp+快速幂+MTT+高维多项式}
此外,插值求值请见 从拉插到快速插值求值
如果对生成函数运算的定义方式与合法性感兴趣,可见 : rqy : 浅谈 OI 中常用的一些生成函数运算的合法与正确性
FFT 对单位根的依赖,决定该算法必须采用浮点数,进而引发精度问题。
糟糕的是,数学家们已经证明,在复数域
大部分的计数题都是在模意义下完成的,我们自然希望为 FFT 找一个模意义下的替代品。
考虑引入模意义同余系
可能有帮助的文章 : 原根&离散对数相关
我们先罗列一下 FFT 中用到的
①
这正是单位根定义的一部分。
由此,我们只需要找出一个满足要求的
②
否则点值就会有重复,插值就会不正确。
③
可以导出对称引理 :
当
综合①②③,等价于 :
④
等价于
我们来看看原根的定义。
众所周知,对于素数
对于
显然,
但是,能够发现,
所以,当
当
这说明,
④ :
所以,
前文提到过,当且仅当
比如
熟练的选手会记下一些常用模数的原根。
懒人福利 : 安利大佬的原根表。
至于如何寻找原根,可见上面拉链的那篇文章。
但是,一定非要原根不可吗?
摆脱原根的方法来自 Uoj : hly1204
假如我们已经知道
随机出一个二次非剩余
那么有
则有
好的,我们找到单位根的替代品了,现在我们来改代码。
上次讲的 FFT 的最终版本(没有“三次变两次优化”):
Old评测记录
#include<algorithm>
#include<cstdio>
#include<cmath>
#define Maxn 1350000
using namespace std;
const double Pi=acos(-1);
int n,m;
struct CP
{
CP (double xx=0,double yy=0){x=xx,y=yy;}
double x,y;
CP operator + (CP const &B) const
{return CP(x+B.x,y+B.y);}
CP operator - (CP const &B) const
{return CP(x-B.x,y-B.y);}
CP operator * (CP const &B) const
{return CP(x*B.x-y*B.y,x*B.y+y*B.x);}
}f[Maxn<<1],p[Maxn<<1];
int tr[Maxn<<1];
void fft(CP *f,bool flag)
{
for (int i=0;i<n;i++)
if (i<tr[i])swap(f[i],f[tr[i]]);
for(int p=2;p<=n;p<<=1){
int len=p>>1;
CP tG(cos(2*Pi/p),sin(2*Pi/p));
if(!flag)tG.y*=-1;
for(int k=0;k<n;k+=p){
CP buf(1,0);
for(int l=k;l<k+len;l++){
CP tt=buf*f[len+l];
f[len+l]=f[l]-tt;
f[l]=f[l]+tt;
buf=buf*tG;
}
}
}
}
int main()
{
scanf("%d%d",&n,&m);
for (int i=0;i<=n;i++)scanf("%lf",&f[i].x);
for (int i=0;i<=m;i++)scanf("%lf",&p[i].x);
for(m+=n,n=1;n<=m;n<<=1);
for(int i=0;i<n;i++)
tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
fft(f,1);fft(p,1);//DFT
for(int i=0;i<n;++i)f[i]=f[i]*p[i];
fft(f,0);
for(int i=0;i<=m;++i)printf("%d ",(int)(f[i].x/n+0.49));
return 0;
}
摇身一变!NTT出现!
还是原来的配方,还是熟悉的味道。
Code
#include<algorithm>
#include<cstdio>
#define ll long long
#define mod 998244353
#define G 3
#define Maxn 1350000
using namespace std;
inline int read()
{
register char ch=0;
while(ch<48||ch>57)ch=getchar();
return ch-'0';
}
ll powM(ll a,ll t=mod-2)
{
ll ans=1;
while(t){
if(t&1)ans=ans*a%mod;
a=a*a%mod;t>>=1;
}return ans;
}
int n,m,tr[Maxn<<1];
ll f[Maxn<<1],g[Maxn<<1],invn;
const ll invG=powM(G);
void NTT(ll *f,bool op)
{
for (int i=0;i<n;i++)
if (i<tr[i])swap(f[i],f[tr[i]]);
for(int p=2;p<=n;p<<=1){
int len=p>>1,
tG=powM(op ? G:invG,(mod-1)/p);
for(int k=0;k<n;k+=p){
ll buf=1;
for(int l=k;l<k+len;l++){
int tt=buf*f[len+l]%mod;
f[len+l]=f[l]-tt;
if (f[len+l]<0)f[len+l]+=mod;
f[l]=f[l]+tt;
if (f[l]>mod)f[l]-=mod;
buf=buf*tG%mod;
}
}
}
}
int main()
{
scanf("%d%d",&n,&m);n++;m++;
for (int i=0;i<n;i++)f[i]=read();
for (int i=0;i<m;i++)g[i]=read();
for(m+=n,n=1;n<m;n<<=1);
for(int i=0;i<n;i++)
tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
NTT(f,1);NTT(g,1);
for(int i=0;i<n;++i)f[i]=f[i]*g[i]%mod;
NTT(f,0);invn=powM(n);
for(int i=0;i<m-1;++i)
printf("%d ",(int)(f[i]*invn%mod));
return 0;
}
评测记录
未经毒瘤卡常时,单个 NTT
变换的时间为 FFT
的 75%~80%
.
其中这一句话改变得很大:
ll tG=powM(op==1 ? G:invG,(mod-1)/p);
即
一些额外的技巧:
预处理单位根,这样可以省去部分乘法操作,并访问一段连续内存。
注意到只会有加减法操作,可以使用 ull
存储,能耐受 18*mod*mod
的范围,在常规范围下可以不取模。
模板题范围较大,在中间取一次模即可。
在不同的题目中,可以加速 15%~40%
不等。
#include<algorithm>
#include<cstdio>
#define ll long long
#define ull unsigned long long
using namespace std;
const int _G=3,mod=998244353,Maxn=1050000;
inline int read(){
int X=0;char ch=0;
while(ch<48||ch>57)ch=getchar();
while(ch>=48&&ch<=57)X=X*10+(ch^48),ch=getchar();
return X;
}
ll powM(ll a,int t=mod-2){
ll ans=1;
while(t){
if(t&1)ans=ans*a%mod;
a=a*a%mod;t>>=1;
}return ans;
}
const int invG=powM(_G);
int tr[Maxn<<1];
void NTT(int *g,bool op,int n)
{
static ull f[Maxn<<1],w[Maxn<<1]={1};
for (int i=0;i<n;i++)f[i]=g[tr[i]];
for(int l=1;l<n;l<<=1){
ull tG=powM(op?_G:invG,(mod-1)/(l+l));
for (int i=1;i<l;i++)w[i]=w[i-1]*tG%mod;
for(int k=0;k<n;k+=l+l)
for(int p=0;p<l;p++){
int tt=w[p]*f[k|l|p]%mod;
f[k|l|p]=f[k|p]+mod-tt;
f[k|p]+=tt;
}
if (l==(1<<17))
for (int i=0;i<n;i++)f[i]%=mod;
}if (!op){
ull invn=powM(n);
for(int i=0;i<n;++i)
g[i]=f[i]%mod*invn%mod;
}else for(int i=0;i<n;++i)g[i]=f[i]%mod;
}
int n,m,F[Maxn<<1],G[Maxn<<1];
int main()
{
n=read()+1;m=read()+1;
for (int i=0;i<n;i++)F[i]=read();
for (int i=0;i<m;i++)G[i]=read();
int len=1;for (;len<n+m;len<<=1);
for(int i=0;i<len;i++)
tr[i]=(tr[i>>1]>>1)|((i&1)?len>>1:0);
NTT(F,1,len);NTT(G,1,len);
for(int i=0;i<len;++i)F[i]=1ll*F[i]*G[i]%mod;
NTT(F,0,len);
for (int i=0;i<n+m-1;i++)
printf("%d ",(int)F[i]);
return 0;
}
多项式模板丰富之后,如何封装成为令人头痛的大难题。
前面讲过,利用 FFT(NTT)
变换的线性性可以减少运算。
这样就需要适度封装以利用性质卡常数,下面给出比较智能的封装:
inline int read(){
int X=0;char ch=0;
while(ch<48||ch>57)ch=getchar();
while(ch>=48&&ch<=57)X=X*10+(ch^48),ch=getchar();
return X;
}
inline void print(ll *f,int len){
for (int i=0;i<len;i++)
printf("%lld ",f[i]);
puts("");
}
#define ll long long
#define ull unsigned long long
#define clr(f,n) memset(f,0,sizeof(int)*(n))
#define cpy(f,g,n) memcpy(f,g,sizeof(int)*(n))
const int _G=3,mod=998244353,Maxn=;
ll powM(ll a,ll t=mod-2){
ll ans=1;
while(t){
if(t&1)ans=ans*a%mod;
a=a*a%mod;t>>=1;
}return ans;
}
const int invG=powM(_G);
int tr[Maxn<<1],tf;
void tpre(int n){
if (tf==n)return ;tf=n;
for(int i=0;i<n;i++)
tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
}
void NTT(int *g,bool op,int n)
{
tpre(n);
static ull f[Maxn<<1],w[Maxn<<1]={1};
for (int i=0;i<n;i++)f[i]=(((ll)mod<<5)+g[tr[i]])%mod;
for(int l=1;l<n;l<<=1){
ull tG=powM(op?_G:invG,(mod-1)/(l+l));
for (int i=1;i<l;i++)w[i]=w[i-1]*tG%mod;
for(int k=0;k<n;k+=l+l)
for(int p=0;p<l;p++){
int tt=w[p]*f[k|l|p]%mod;
f[k|l|p]=f[k|p]+mod-tt;
f[k|p]+=tt;
}
if (l==(1<<10))
for (int i=0;i<n;i++)f[i]%=mod;
}if (!op){
ull invn=powM(n);
for(int i=0;i<n;++i)
g[i]=f[i]%mod*invn%mod;
}else for(int i=0;i<n;++i)g[i]=f[i]%mod;
}
void px(int *f,int *g,int n)
{for(int i=0;i<n;++i)f[i]=1ll*f[i]*g[i]%mod;}
void times(int *f,int *g,int len,int lim)
{
static int sav[Maxn<<1];
int n=1;for(n;n<len+len;n<<=1);
clr(sav,n);cpy(sav,g,n);
NTT(f,1,n);NTT(sav,1,n);
px(f,sav,n);NTT(f,0,n);
clr(f+lim,n-lim);clr(sav,n);
}
熟练的选手能够在 8~15min
内默完,保险起见建议拿 times
跟暴力卷积拍一下。
这样我们就能够肆无忌惮的使用 NTT
和一系列缩写了。
小心地分析清零带来的影响,宁愿多清不能漏网,而且这种问题一般很难检查和拍出来。
善用宏可以有效避免重名问题。
本文所讨论的多项式均为形式幂级数,但限于篇幅,并没有给出其严格定义和应用。
它只是辅助分析的数学工具,所谓多项式中的“未知数”并没有实际意义,只是一个代数对象。
如果想了解更多相关内容,可见 rqy : 浅谈 OI 中常用的一些生成函数运算的合法与正确性。
简化地讲,我们将以多项式加法和加法卷积来定义所有运算。
在大多数题目中,我们只对多项式的前若干项感兴趣,此时我们为所有运算设定一个次数上界,即
由于多项式加法和乘法转移的单向性(只会从低次向高次贡献),不难发现下列两条性质:
前面提到过,所有运算都是利用加法和乘法定义的,所以这能说明,我们可以在每一步运算中都保持次数界,利用其去除无用的高次项。
这样就能避免对无穷幂级数的处理,以保障复杂度。
终于从乘法向多项式全家桶迈进了!
约定 :
也就是说,
多项式的加减定义为 :
其实就是简单地将对应项相加减,
多项式乘法(加法卷积)定义为 :
使用
目前为止,加法,减法,乘法都和我们平时多项式运算的经验保持一致。
接下来就是多项式乘法逆元,定义为 :
满足
类比同余系下数字的逆元,多项式逆元可以把除法转换成乘法。
怎么求一个多项式的乘法逆元呢?
方法一 : 递推
由
可以
这个递推算法并没有用到 NTT,所以不可能达到
方法二 : 倍增
我们考虑采用倍增的思想,不断扩大次数上界。
当界为
假设我们已经得到了
设
明显有
如何把
得
展开得到
等式两边同乘
移项得到
根据这个即可从
复杂度为
Code : (配合前面的封装食用)
void invp(int *f,int m)
{
int n;for (n=1;n<m;n<<=1);
static int w[Maxn<<1],r[Maxn<<1],sav[Maxn<<1];
w[0]=powM(f[0]);
for (int len=2;len<=n;len<<=1){
for (int i=0;i<(len>>1);i++)
r[i]=(w[i]<<1)%mod;
cpy(sav,f,len);
NTT(w,1,len<<1);px(w,w,len<<1);
NTT(sav,1,len<<1);px(w,sav,len<<1);
NTT(w,0,len<<1);clr(w+len,len);
for (int i=0;i<len;i++)
w[i]=(r[i]-w[i]+mod)%mod;
}cpy(f,w,m);clr(sav,n+n);clr(w,n+n);clr(r,n+n);
}
int n,F[Maxn<<1];
int main()
{
n=read();
for (int i=0;i<n;i++)F[i]=read();
invp(F,n);print(F,n);
return 0;
}
注意,由于中间 NTT
是两倍长度,最后清空数组要清空两倍的 n
。
检查的方法 : 随机一个多项式求逆两次看看是不是自己。
评测记录
众所周知
观察倍增公式
我们最终只需要得到答案的后
先计算
这部分结果的前
再乘上一个
void invp(int *f,int m)
{
int n;for (n=1;n<m;n<<=1);
static int w[Maxn<<1],r[Maxn<<1],sav[Maxn<<1];
w[0]=powM(f[0]);
for (int len=2;len<=n;len<<=1){
for (int i=0;i<(len>>1);i++)r[i]=w[i];
cpy(sav,f,len);NTT(sav,1,len);
NTT(r,1,len);px(r,sav,len);
NTT(r,0,len);clr(r,len>>1);
cpy(sav,w,len);NTT(sav,1,len);
NTT(r,1,len);px(r,sav,len);
NTT(r,0,len);
for (int i=len>>1;i<len;i++)
w[i]=(w[i]*2ll-r[i]+mod)%mod;
}cpy(f,w,m);clr(sav,n);clr(w,n);clr(r,n);
}
评测记录 (可以加速两倍多)
继续学习之前,建议读者先自行了解简单微积分。
首先定义多项式的求导 :
定义多项式的积分 :
这和经典的求导、积分是一致的。
Code:
void dao(int *f,int m){
for (int i=1;i<m;i++)
f[i-1]=1ll*f[i]*i%mod;
f[m-1]=0;
}
int inv[Maxn];
void Init(int lim){
inv[1]=1;
for (int i=2;i<=lim;i++)
inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod;
}
void jifen(int *f,int m){
for (int i=m;i;i--)
f[i]=1ll*f[i-1]*inv[i]%mod;
f[0]=0;
}
定义多项式复合为 :
将所有
用于解决下列问题:
已知函数
实际应用中函数
证明可在 多项式计数杂谈 中查看。
结论 :
也就是说我们采用这个式子倍增地求解上述问题。
这个式子很牛逼,可以帮我们省下思维难度(肝),一定好好记住。
需要注意的是,
比如说上面的 多项式求逆
根据
设
所以
(下文默认
得
注意
得
这与我们上面采用平方法推出的式子是等价的。
多项式开方
P5205 【模板】多项式开根
得
令
则
根据牛顿迭代得到 :
根据这个倍增即可。
多项式只有常数时,开方的结果是二次剩余,题目中保证常数为
Code:
void sqrtp(int *f,int m)
{
int n;for (n=1;n<m;n<<=1);
static int b1[Maxn<<1],b2[Maxn<<1];
b1[0]=1;
for (int len=2;len<=n;len<<=1){
for (int i=0;i<(len>>1);i++)b2[i]=(b1[i]<<1)%mod;
invp(b2,len);
NTT(b1,1,len);px(b1,b1,len);NTT(b1,0,len);
for (int i=0;i<len;i++)b1[i]=(f[i]+b1[i])%mod;
times(b1,b2,len,len);
}cpy(f,b1,m);clr(b1,n+n);clr(b2,n+n);
}
评测记录。
P4512 【模板】多项式除法
如果保证没有余多项式(即整除),直接使用普通多项式除法就能得到商多项式。
我们考虑用一种奥妙重重的方式把余多项式去掉。
前面我们经常利用
然而此处余多项式的次数低,考虑反转系数,将原来次数高的变成次数低的。
我们定义
手玩一下就会发现
比如说
首先有
换元得到
两边同乘
写成反转的形式,得 :
此时我们机智地
变形得
这就可以直接用多项式求逆做了。求出
两者均由麦克劳林级数定义 :
为啥非得要用不友好的麦克劳林级数?因为这样就能做到仅使用加法和乘法来定义运算了。
经典的性质在这个定义下仍然成立,如
P4725 【模板】多项式对数函数(多项式 ln)
当然,可以直接使用定义式来计算答案,但是复杂度会十分糟糕。
考虑利用更加简明的性质,如
题目有
同时积分得到
也就是说,一个求导,一个求逆,一个乘法,一个积分就好了。
Code:
void lnp(int *f,int m)
{
static int g[Maxn<<1];
cpy(g,f,m);
invp(g,m);dao(f,m);
times(f,g,m,m);
jifen(f,m-1);
clr(g,m);
}
评测记录
P4726 【模板】多项式指数函数(多项式 exp)
递推
有
提取系数可得
即
每次可以
其实,由
由
倍增
和
回忆式子
根据
求导得
代入
根据上式倍增即可。
注意,必须保证原多项式
Code:
void exp(int *f,int m)
{
static int s[Maxn<<1],s2[Maxn<<1];
int n=1;for(;n<m;n<<=1);
s2[0]=1;
for (int len=2;len<=n;len<<=1){
cpy(s,s2,len>>1);lnp(s,len);
for (int i=0;i<len;i++)
s[i]=(f[i]-s[i]+mod)%mod;
s[0]=(s[0]+1)%mod;
times(s2,s,len,len);
}cpy(f,s2,m);clr(s,n);clr(s2,n);
}
评测记录
观察上面的递推式,不难发现是个半在线卷积,分治FFT
即可。
半在线卷积具体怎么做可见 : 半在线卷积小记
虽然复杂度是
P5245 【模板】多项式快速幂
这里是
公式十分简单:
( 也就是把乘法转化成指数上的加法 )
( 常数高达一个
Code:
可以当做全套模板来用。
#include<algorithm>
#include<cstring>
#include<cstdio>
#define ll long long
#define ull unsigned long long
#define clr(f,n) memset(f,0,sizeof(int)*(n))
#define cpy(f,g,n) memcpy(f,g,sizeof(int)*(n))
const int _G=3,mod=998244353,Maxn=135000;
int read(){
int X=0;char ch=0;
while(ch<48||ch>57)ch=getchar();
while(ch>=48&&ch<=57)X=X*10+(ch^48),ch=getchar();
return X;
}
inline void print(int *f,int len){
for (int i=0;i<len;i++)
printf("%d ",f[i]);
puts("");
}
ll powM(ll a,ll t=mod-2){
ll ans=1;
while(t){
if(t&1)ans=ans*a%mod;
a=a*a%mod;t>>=1;
}return ans;
}
const int invG=powM(_G);
int tr[Maxn<<1],tf;
void tpre(int n){
if (tf==n)return ;tf=n;
for(int i=0;i<n;i++)
tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
}
void NTT(int *g,bool op,int n)
{
static ull f[Maxn<<1],w[Maxn]={1};
tpre(n);
for (int i=0;i<n;i++)f[i]=(((ll)mod<<5)+g[tr[i]])%mod;
for(int l=1;l<n;l<<=1){
ull tG=powM(op?_G:invG,(mod-1)/(l+l));
for (int i=1;i<l;i++)w[i]=w[i-1]*tG%mod;
for(int k=0;k<n;k+=l+l)
for(int p=0;p<l;p++){
int tt=w[p]*f[k|l|p]%mod;
f[k|l|p]=f[k|p]+mod-tt;
f[k|p]+=tt;
}
if (l==(1<<10))
for (int i=0;i<n;i++)f[i]%=mod;
}if (!op){
ull invn=powM(n);
for(int i=0;i<n;++i)
g[i]=f[i]%mod*invn%mod;
}else for(int i=0;i<n;++i)g[i]=f[i]%mod;
}
void px(int *f,int *g,int n)
{for(int i=0;i<n;++i)f[i]=1ll*f[i]*g[i]%mod;}
int _g1[Maxn<<1];
#define sav _g1
void times(int *f,int *g,int len,int lim)
{
int n=1;while(n<len+len)n<<=1;
cpy(sav,g,n);clr(sav+len,n-len);
NTT(f,1,n);NTT(sav,1,n);
px(f,sav,n);NTT(f,0,n);
clr(f+lim,n-lim);clr(sav,n);
}
void invp(int *f,int m)
{
int n;for (n=1;n<m;n<<=1);
static int w[Maxn<<1],r[Maxn<<1];
w[0]=powM(f[0]);
for (int len=2;len<=n;len<<=1){
for (int i=0;i<(len>>1);i++)r[i]=w[i];
cpy(sav,f,len);NTT(sav,1,len);
NTT(r,1,len);px(r,sav,len);
NTT(r,0,len);clr(r,len>>1);
cpy(sav,w,len);NTT(sav,1,len);
NTT(r,1,len);px(r,sav,len);
NTT(r,0,len);
for (int i=len>>1;i<len;i++)
w[i]=(w[i]*2ll-r[i]+mod)%mod;
}cpy(f,w,m);clr(sav,n);clr(w,n);clr(r,n);
}
void dao(int *f,int m){
for (int i=1;i<m;i++)
f[i-1]=1ll*f[i]*i%mod;
f[m-1]=0;
}
int inv[Maxn];
void Init(int lim){
inv[1]=1;
for (int i=2;i<=lim;i++)
inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod;
}
void jifen(int *f,int m){
for (int i=m;i;i--)
f[i]=1ll*f[i-1]*inv[i]%mod;
f[0]=0;
}
void lnp(int *f,int m)
{
static int g[Maxn<<1];
cpy(g,f,m);
invp(g,m);dao(f,m);
times(f,g,m,m);
jifen(f,m-1);
clr(g,m);
}
void exp(int *f,int m)
{
static int s[Maxn<<1],s2[Maxn<<1];
int n=1;for(;n<m;n<<=1);
s2[0]=1;
for (int len=2;len<=n;len<<=1){
cpy(s,s2,len>>1);lnp(s,len);
for (int i=0;i<len;i++)
s[i]=(f[i]-s[i]+mod)%mod;
s[0]=(s[0]+1)%mod;
times(s2,s,len,len);
}cpy(f,s2,m);clr(s,n);clr(s2,n);
}
char str[Maxn];
int k,f[Maxn<<1],m;
int main()
{
scanf("%d%s",&m,str);Init(m);
for (int i=0;'0'<=str[i]&&str[i]<='9';i++)
k=(10ll*k+str[i]-'0')%mod;
for (int i=0;i<m;i++)f[i]=read();
lnp(f,m);
for(int i=0;i<m;i++)f[i]=1ll*f[i]*k%mod;
exp(f,m);print(f,m);
return 0;
}
评测记录
设原多项式
对
将
拿这两个式子搞搞能得到
提取系数可得 :
挖出次数最高的项 :
把
前面的和式枚举复杂度都是
//默认f[0]=1
void powp(int *f,int k,int n,int p)
{
static int g[Maxn]={1};
for (int m=0;m<n-1;m++){
for (int i=1;i<k;i++)
g[m+1]=g[m+1]+(1ll*p*i-m+i-1)%mod*f[i]*g[m-i+1];
g[m+1]=(ll)(g[m+1]+mod)*inv[m+1]%mod;
}cpy(f,g,n);
}
这里也要求
( 令
注意,当 Ln+EXP
快,后者只能跑
注意到其不保证满足
考虑把
不难发现
注意,常数
求出
int main()
{
scanf("%d%s",&m,str);
for (int i=0;i<m;i++)f[i]=read();
int p=0,c;while(!f[p])p++;c=powM(f[p]);
for (int i=0;'0'<=str[i]&&str[i]<='9';i++){
k=(10ll*k+str[i]-'0')%mod;
k2=(10ll*k2+str[i]-'0')%(mod-1);
if (1ll*k*p>m){
for (int i=0;i<m;i++)printf("0 ");
return 0;
}
}Init(m=m-p*k);
for (int i=0;i<m;i++)f[i]=1ll*f[i+p]*c%mod;
clr(f+m,p*k);lnp(f,m);
for(int i=0;i<m;i++)f[i]=1ll*f[i]*k%mod;
exp(f,m);
for (int i=0;i<p*k;i++)printf("0 ");
c=powM(c,mod-1-k2);
for (int i=0;i<m;i++)f[i]=1ll*f[i]*c%mod;
print(f,m);
return 0;
}
评测记录
讲了这么多,我们把几个算法的核心式子都集合一下。
(其实主要是背代码)
可以通过ull
以及预处理单位根优化。
(带*的表示
已知
倍增公式 :
定义一个
得
变化为反转的形式,得
要求常数项为
倍增公式
要求常数项为
附送一个 vector
板子 :
#include<algorithm>
#include<cstring>
#include<cstdio>
#include<vector>
#define ll long long
#define ull unsigned long long
#define clr(f,n) memset(f,0,sizeof(int)*(n))
#define cpy(f,g,n) memcpy(f,g,sizeof(int)*(n))
const int _G=3,mod=998244353,Maxn=1<<17|500;
using namespace std;
ll powM(ll a,ll t=mod-2){
ll ans=1;
while(t){
if(t&1)ans=ans*a%mod;
a=a*a%mod;t>>=1;
}return ans;
}
const int invG=powM(_G);
int tr[Maxn<<1],tf;
void tpre(int n){
if (tf==n)return ;tf=n;
for(int i=0;i<n;i++)
tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
}
void NTT(int *g,bool op,int n)
{
tpre(n);
static ull f[Maxn<<1],w[Maxn<<1];w[0]=1;
for (int i=0;i<n;i++)f[i]=(((ll)mod<<5)+g[tr[i]])%mod;
for(int l=1;l<n;l<<=1){
ull tG=powM(op?_G:invG,(mod-1)/(l+l));
for (int i=1;i<l;i++)w[i]=w[i-1]*tG%mod;
for(int k=0;k<n;k+=l+l)
for(int p=0;p<l;p++){
int tt=w[p]*f[k|l|p]%mod;
f[k|l|p]=f[k|p]+mod-tt;
f[k|p]+=tt;
}
if (l==(1<<10))
for (int i=0;i<n;i++)f[i]%=mod;
}if (!op){
ull invn=powM(n);
for(int i=0;i<n;++i)
g[i]=f[i]%mod*invn%mod;
}else for(int i=0;i<n;++i)g[i]=f[i]%mod;
}
void px(int *f,int *g,int n)
{for(int i=0;i<n;++i)f[i]=1ll*f[i]*g[i]%mod;}
#define Poly vector<int>
Poly operator + (const Poly &A,const Poly &B)
{
Poly C=A;C.resize(max(A.size(),B.size()));
for (int i=0;i<B.size();i++)C[i]=(C[i]+B[i])%mod;
return C;
}
Poly operator - (const Poly &A,const Poly &B)
{
Poly C=A;C.resize(max(A.size(),B.size()));
for (int i=0;i<B.size();i++)C[i]=(C[i]+mod-B[i])%mod;
return C;
}
Poly operator * (const int c,const Poly &A)
{
Poly C;C.resize(A.size());
for (int i=0;i<A.size();i++)C[i]=1ll*c*A[i]%mod;
return C;
}
int lim;
Poly operator * (const Poly &A,const Poly &B)
{
static int a[Maxn<<1],b[Maxn<<1];
for (int i=0;i<A.size();i++)a[i]=A[i];
for (int i=0;i<A.size();i++)a[i]=A[i];
cpy(a,&A[0],A.size());
cpy(b,&B[0],B.size());
Poly C;C.resize(min(lim,(int)(A.size()+B.size()-1)));
int n=1;for(n;n<A.size()+B.size()-1;n<<=1);
NTT(a,1,n);NTT(b,1,n);
px(a,b,n);NTT(a,0,n);
cpy(&C[0],a,C.size());
clr(a,n);clr(b,n);
return C;
}
void pinv(const Poly &A,Poly &B,int n)
{
if (n==1)B.push_back(powM(A[0]));
else if (n&1){
pinv(A,B,--n);
int sav=0;
for (int i=0;i<n;i++)sav=(sav+1ll*B[i]*A[n-i])%mod;
B.push_back(1ll*sav*powM(mod-A[0])%mod);
}else {
pinv(A,B,n/2);
Poly sA;sA.resize(n);
cpy(&sA[0],&A[0],n);
B=2*B-B*B*sA;
B.resize(n);
}
}
Poly pinv(const Poly &A)
{
Poly C;pinv(A,C,A.size());
return C;
}
int inv[Maxn];
void Init()
{
inv[1]=1;
for (int i=2;i<=lim;i++)
inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod;
}
Poly dao(const Poly &A)
{
Poly C=A;
for (int i=1;i<C.size();i++)
C[i-1]=1ll*C[i]*i%mod;
C.pop_back();
return C;
}
Poly ints(const Poly &A)
{
Poly C=A;
for (int i=C.size()-1;i;i--)
C[i]=1ll*C[i-1]*inv[i]%mod;
C[0]=0;
return C;
}
Poly ln(const Poly &A)
{return ints(dao(A)*pinv(A));}
void pexp(const Poly &A,Poly &B,int n)
{
if (n==1)B.push_back(1);
else if (n&1){
pexp(A,B,n-1);n-=2;
int sav=0;
for (int i=0;i<=n;i++)sav=(sav+1ll*(i+1)*A[i+1]%mod*B[n-i])%mod;
B.push_back(1ll*sav*inv[n+1]%mod);
}else {
pexp(A,B,n/2);
Poly lnB=B;
lnB.resize(n);lnB=ln(lnB);
for (int i=0;i<lnB.size();i++)
lnB[i]=(mod+A[i]-lnB[i])%mod;
lnB[0]++;
B=B*lnB;
B.resize(n);
}
}
Poly pexp(const Poly &A)
{
Poly C;pexp(A,C,A.size());
return C;
}
Poly F;
int main()
{
int n;scanf("%d",&n);
lim=n;Init();
F.resize(n);
for (int i=0;i<F.size();i++)scanf("%d",&F[i]);
F=pexp(F);
for (int i=0;i<F.size();i++)printf("%d ",F[i]);
return 0;
}
本部分未完。
先来介绍二元幂级数,即形如
二元幂级数的乘法仍然定义为在各个维度上的加法卷积。
有
二元幂级数的“系数”可以仍然是幂级数,如
也有
可以先对
然后上式中
这可以得到二元幂级数卷积的过程 :
先对行做
也可以通过二元函数插值来解释。
有了卷积,我们再来推导多项式全家桶。
求逆仍然可以使用倍增公式
令
这样,虽然单个
这样,只需记录模
实现时,第二维可以单独封装成一个结构体。这样也可以直接套用一元牛顿迭代。
板子 :[数学记录]Uoj#596. 【集训队互测2021】三维立体混元劲
P4245 【模板】任意模数NTT
有的时候题目给的膜数并不是 看似平淡无奇实则暗藏杀机)。
那么,我们苦苦研究的 NTT 就都没有用武之地了吗?
任意模数多项式乘法原理 : 假设卷积前,每个数都小于
卷积后,能产生的最大的数也就是
我们只要弄一个精度够高的多项式乘法,就能做到不丢精度。
一种想法是 : 跑九次 NTT (三次卷积),把答案用中国剩余定理合并,精度可达
这种做法常数非常大,而且 long long
存不下,还需要一些黑科技辅助,所以不推荐。
另一种想法是 : 拆系数 FFT ,又称为 MTT 。
把一个数拆成
把
那么
这样做需要
冷静分析 : 每个多项式只用插值一次,共
点乘复杂度忽略,最后得到的四个多项式各插值一次,共
这样就是
我们考虑推推式子来优化:
我们有四个多项式
考虑到
那么我们设复多项式
那么
我们又设
那么