FFT & NTT学习心得

基础知识浅谈

FFT—快速傅里叶变换

基本功能

在O( (n+m)log(n+m) )的时间复杂度内计算：n次多项式乘m次多项式.

实现方式

欲求多项式A*多项式B：
1.对多项式A，B分别进行快速傅里叶变换,得到A1,B1;
2.将A1,B1的对应项相乘得到多项式C1.即：C1[i]=A1[i]*B1[i] 其中C1[i]表示C1的第i项系数.
3.将多项式C1进行逆变换得到多项式C.则多项式C即是多项式A与多项式B相乘的结果.
具体证明这里不多说了.

NTT—快速数论变换

基本功能

与FFT类似，但没有精度差.

实现方式

与FFT类似，只是用数论原根的幂代替FFT中的w[n].

原根的含义及求法

根据欧拉定理，若a,p互质则有：

aϕ(p)≡1  (mod  p)

如果当且仅当x=ϕ(p) 时

ax≡1  (mod  p)

则a是p的一个原根.
用更专业的语言来描述即：a模p的阶等于ϕ(p) 则a是p的原根.
原根的求法：没有什么特别的方法，只能暴力枚举，验证时略有技巧，这里就不细说了.

关于NTT算法的一些说明

优点

相比起FFT来说，NTT最显著的优势在于没有精度差.由于FFT会用到complex double ,在大数据下不排除出现精度差的可能.在某些评测机上，效率可能不如NTT.(NTT虽然不用复数运算，但是取模很多).

限制

相比起FFT来说,NTT的限制很多.

• 所求的多项式要求是整系数
• 如果题目要求结果对质数p取模,这个质数往往只能是998244353，否则会有很多麻烦，这个会在后面谈到.
• 所求多项式的项数应在223 之内,因为998244353=7∗17∗223+1
• 结果的系数不应超过质数P.(P是自己选择的质数,一般定为P=998244353)

具体代码实现

FFT模板

const double pi=3.1415926535897932384626433832795;complex<double>Wi[MAXN];char s[MAXN],t[MAXN];void FFT(complex<double> A[],int nn,int ty){    int i,j,k,m;    complex<double> t0,t1;    for(i=0;i<nn;i++)    {        for(j=0,k=i,m=1;m<nn;m<<=1,j=(j<<1)|(k&1),k>>=1);        if(i<j)t0=A[i],A[i]=A[j],A[j]=t0;    }//这段for循环不建议擅自改动,极易出错.    Wi[0]=1;    for(m=1;m<nn;m<<=1)    {        t0=exp(complex<double>(0,ty*pi/m));        for(i=1;i<m;i++)            Wi[i]=Wi[i-1]*t0;        for(k=0;k<nn;k+=m<<1)            for(i=k;i<k+m;i++)            {                t0=A[i];                t1=A[i+m]*Wi[i-k];                A[i]=t0+t1;                A[i+m]=t0-t1;            }    }    if(ty==1)        return;//ty==-1时为逆变换.    t0=1.0/nn;    for(i=0;i<nn;i++)        A[i]*=t0;}

NTT模板

const int P=998244353;const int g=3;//P的原根int W[MAXN];int exp(int a,int k){    ll A=1LL*a,ANS=1LL;    for(;k;k>>=1,A=A*A%P)    {        if(k&1)        {            ANS=ANS*A%P;        }    }    return (int)ANS%P;}//其实这里的函数名取为exp是不恰当的，不过无伤大雅.void NTT(int A[],int nn,int ty){    int t1,t2,i,j,k,m;    for(i=0;i<nn;i++)    {        for(j=0,k=i,m=1;m<nn;m<<=1,j=(j<<1)|(k&1),k>>=1);        if(i<j)        {            t1=A[i];            A[i]=A[j];            A[j]=t1;        }    }    W[0]=1;    for(m=1;m<nn;m<<=1)    {        t1=exp(g,P-1+ty*(P-1)/(m<<1));        for(i=1;i<m;i++)        {            W[i]=1LL*W[i-1]*t1%P;        }        for(k=0;k<nn;k+=m<<1)        {            for(i=k;i<k+m;i++)            {                t1=A[i];                t2=1LL*A[i+m]*W[i-k]%P;                A[i]=t1+t2;                A[i]-=A[i]>P?P:0;                A[i+m]=t1-t2;                A[i+m]+=A[i+m]<0?P:0;            }        }    }    if(ty==1)    {        return ;    }    t1=exp(nn,P-2);    for(i=0;i<nn;i++)    {        A[i]=1LL*A[i]*t1%P;    }    return ;}

FFT及NTT的简单应用

最基本应用–高精度乘法

将十进制数看成一个多项式,然后利用FFT或NTT求解.(注意进位).
比如：
123*456=56088可以看成:

123=1∗102+2∗101+3∗100

456=4∗102+5∗101+6∗100

123∗456=4∗104+13∗103+28∗102+27∗101+18∗100

考虑进位则有：

123∗456=5∗104+6∗103+0∗102+8∗101+8∗100=56088

稍微高端一点的应用

特殊的计数问题

题目大意

给出两个长度分别为n,m的01序列A、B，有Q次询问，每次询问要求回答：当把A的第i位和B的第j位对齐时，A,B公共部分有多少对对齐且相同的数.
如：
输入n,m,A串,B串,Q,每次询问(i,j)
6 6
000111
111100
3
1 1
1 2
4 2
输出每次询问的答案:
1
0
3

数据范围:

n<=100000,Q<=1000000

解题方法

注意到所有询问总共只有2n-1种本质不同的询问,我们希望预处理出所有这些本质不同的询问的答案.
比如询问为(i,j)我们可以等价转化为(1,j-i)，因此我们不妨将询问设为(1,k)
我们先考虑有多少位上的1是对齐的：
其实两个位置上的1对齐,当且仅当两个位置都是1,且两个位置相差k-1.
如果将其中一个序列倒序排列的话，那么两个之前对齐的位置现在的位置标号之和就是定值了.于是就可以构造多项式进行计数.
这之后，我们将所有0，1取反,就可以同样计算出有多少位0是对齐的.

参考代码

#include<cstdio>#include<cstring> #define MAXN 425000const int P=998244353;const int G=3;char sa[MAXN],sb[MAXN];int N,W[MAXN*2];void _r(int& x){    char c=getchar();    while(c<'0'||c>'9')    {        c=getchar();    }    for(x=0;c>='0'&&c<='9';c=getchar())    {        x=(x<<1)+(x<<3)+c-'0';    }    return ;}bool check(int x[],int nn){    for(int i=0;i<nn;i++)    {        if(x[i])        {            return false;        }    }    return true;}int CNT[MAXN*2];int exp(int a,int k){    int ans=1;    for(;k;k>>=1)    {        if(k&1)        {            ans=1ll*ans*a%P;        }        a=1ll*a*a%P;    }    return ans;}void fft(int A[],int nn,int ty){    int t1,t2,i,j,k,m;    for(i=0;i<nn;i++)    {        for(j=0,k=i,m=1;m<nn;m<<=1,j=(j<<1)|(k&1),k>>=1);        if(i<j)        {            t1=A[i];            A[i]=A[j];            A[j]=t1;        }    }    W[0]=1;    for(m=1;m<nn;m<<=1)    {        t1=exp(G,P-1+ty*(P-1)/(m<<1));        for(i=1;i<m;i++)        {            W[i]=1ll*W[i-1]*t1%P;        }        for(k=0;k<nn;k+=m<<1)        {            for(i=k;i<k+m;i++)            {                t1=A[i];                t2=1ll*A[i+m]*W[i-k]%P;                A[i]=t1+t2;                A[i]-=A[i]>P?P:0;                A[i+m]=t1-t2;                A[i+m]+=A[i+m]<0?P:0;            }        }    }    if(ty==1)    {        return ;    }    t1=exp(nn,P-2);    for(i=0;i<N;i++)    {        A[i]=1ll*A[i]*t1%P;    }    return ;}//其实是NTT只是函数名取成了fftint A[MAXN*2],B[MAXN*2],n,m,q;int o[12],tot=0;void PUT(int x){    if(x==0)    {        putchar('0');        return ;    }    for(tot=0;x;x/=10)    {        o[++tot]=x%10;    }    for(;tot;tot--)    {        putchar(o[tot]+'0');    }    return ;}int main(){    _r(n);    _r(m);    scanf("%s%s",sa,sb);    for(int i=0;i<n;i++)    {        A[i]=sa[i]-'0';    }    for(int i=0;i<m;i++)    {        B[i]=sb[m-i-1]-'0';    }    for(N=1;N<n+m+1;N<<=1);    if(check(A,N)||check(B,N));//特判两个多项式是否存在为0的,否则用NTT会有诈.    else    {        fft(A,N,1);        fft(B,N,1);        for(int i=0;i<N;i++)        {            A[i]=1ll*A[i]*B[i]%P;        }        fft(A,N,-1);        for(int i=0;i<N;i++)        {            CNT[i]+=A[i];        }    }    memset(A,0,sizeof(A));    memset(B,0,sizeof(B));    for(int i=0;i<n;i++)    {        A[i]=sa[i]-'0';        A[i]^=1;    }    for(int i=0;i<m;i++)    {        B[i]=sb[m-i-1]-'0';        B[i]^=1;    }    if(check(A,N)||check(B,N));    else    {         fft(A,N,1);        fft(B,N,1);        for(int i=0;i<N;i++)        {            A[i]=1ll*A[i]*B[i]%P;        }        fft(A,N,-1);        for(int i=0;i<N;i++)        {            CNT[i]+=A[i];        }    }     _r(q);    for(int i=1,u,v;i<=q;i++)    {        _r(u);        _r(v);        int p=u-1+m-v;        PUT(CNT[p]);        putchar('\n');    }    return 0;} 

还有一道很牛逼的题目—万径人踪灭,也用到了这个方法.

NTT&FFT与多项式问题

相关知识

• 多项式求逆
• 多项式除法
• 多项式取模
• 多项式开方
• 多项式快速幂
• 多项式取对数
• ………………

这些东西本人也是一知半解，具体算法的证明过程，代码实现可以参考网上的其他博客。这里给一个参考的链接(加载非常慢，但内容非常好，如有需要请耐心等待):http://picks.logdown.com/posts/197262-polynomial-division

简单地谈一谈

多项式求逆

已知多项式A,欲求多项式B,使A∗B≡1 (mod xn).
首先,常数项取B[0]=A[0]−1
假如我们已经得到B’，使得A∗B′≡1 (mod xn),欲求B使得A∗B≡1 (mod x2n)
那么有

A∗(B−B′)≡0 (mod xn)

于是

B−B′≡0 (mod xn)

平方

B2+B′2−2B∗B′≡0 (mod x2n)

同时乘A得

A∗B2+A∗B′2−2A∗B∗B′≡0 (mod x2n)

即

B+A∗B′2−2B′≡0 (mod x2n)

于是得到

B≡2B′−A∗B′2 (mod x2n)

这里用到了倍增的思想.
多项式开方与之类似.
时间复杂度

T(n)=T(n/2)+O(nlog n),即O(nlog n)

一些典型题目

贝壳串

【问题描述】
海边市场有长度分别为1到n的贝壳串出售，其中长度为i的贝壳串有a[i]种，每种贝壳串有无限个，问用这些贝壳串链接成长度为n的串有多少种方案?

【输入格式】
第一行，一整数n，
第二行，n个整数ai表示长度为i的贝壳串的种类数

【输出格式】
输出方案数，结果模313

【输入输出样例】
in
3
1 3 7
out
14

in
4
2 2 2 2
out
54

int
7
86 58 87 145 510 32 263
out
152
【数据范围】
对于50%的数据n<=1000
对于100%的数据n<=100000，0<=ai<=10000000

官方题解

FFT+分治

另解

FFT+多项式求逆

然而我是一个异端

题目要求对313取模,本是断绝了NTT的路,但我硬是强用NTT把这个题过了！
做法非常的“淫荡”.为此我专门写了一篇题解：

另类题解

【贝壳串shell】的另类解法—-强行NTT.

1.DP方程
设f[i]为配成长度为i的项链的方案数,a[i]表示长度为i的贝壳的种类数.
枚举最后一串贝壳的长度,可得转移方程：

f[i]=∑i=1nf[n−i]⋅a[i]

并规定f[0]=1.
直接暴力DP可以拿到50分.

2.构造母函数

F(x)=∑i=0+∞f[i]⋅xi        A(x)=∑i=0+∞a[i]⋅xi

于是根据dp方程有:

F(x)=F(x)∗A(x)+1

等价于:

F(x)=11−A(x)

到此为止用FFT就足以解决问题了.

3.强用NTT
我们知道，多项式求逆需要多次使用NTT，这样不可避免地会使系数对质数P(P=998244353)取模.但是题目要求对313取模,这样得到的答案就不对了.

举个例子吧：
用NTT计算(1+x)^n的一次项系数,答案对9973取模.如果我们直接把多项式进行(n-1)次乘法,那么得到的结果将是n%998244353%9973，这与n%9973是不等价的.
但是如果在每次乘以(x+1)后,将结果的各项系数都对9973取模,得到的将是正确的结果.因为NTT实际上是在大质数的剩余系下进行乘法运算,如果单次NTT变换的结果不会超出大质数,所得结果将与正常的乘法一样.

所以，贝壳串这道题也可以使用这样的办法.

很遗憾的是,就贝壳串一题而言,单次运算的上界=313*313*100000,大约100亿左右,显然超出了998244353之类的数的范围，怎么办？

我们可以在P=998244353*1004535809的剩余系下完成乘法。(998244353,1004535809都是常见的用于NTT的大质数,本人更喜欢用前者.）我们可以在每次NTT之后都把系数对313取模,这样中间过程一直都不会超过P.

但是P=998244353*1004535809不能直接用于NTT，于是我们可以把多项式分别在998244353，1004535809的剩余系下进行NTT变换,之后用中国剩余定理求出系数模P的结果.

不过还有一个坑点,这样得到的结果不一定是真实的结果,不能简单地对313取模.按照以上做法，得到的结果除了可能是真实结果之外,还有可能是(真实结果+P),这时的真实结果是负数.因为在做NTT变换时,负数一般是+P转为正数的.好在负数结果+P的值要远大于正数结果,我们可以很简单区分开二者,然后加以特判即可.其实这个问题非常容易被忽视！考场上我起初是没有发现这个问题的，调了将近一个小时才发现这个问题。

至此这个问题终于讲清楚了,贴上我那丑陋的代码！

//1004535809 998244353#include<cstdio>#define MAXN 450000#define LL long long#define mod 313const int P=998244353,Q=1004535809,G=3;const LL PQ=1ll*998244353*1004535809,lim=1ll*mod*mod*400000;int exp(int a,int k,int MOD){    int an=1;    for(;k;k>>=1)    {        if(k&1)        {            an=1ll*an*a%MOD;        }        a=1ll*a*a%MOD;    }    return an;}const int ip=exp(P,Q-2,Q),iq=exp(Q,P-2,P);int A[MAXN],F[MAXN];int n;void _r(int& x){    char c=getchar();    while(c<'0'||c>'9')    {        c=getchar();    }    for(x=0;c>='0'&&c<='9';c=getchar())    {        x=(x<<1)+(x<<3)+c-'0';    }    return ;}void work1(){    F[0]=1;    for(int i=1;i<=n;i++)    {        for(int j=1;j<=i;j++)        {            F[i]=(F[i]+F[i-j]*A[j]%mod)%mod;        }    }    printf("%d\n",F[n]);    return ;}int W[MAXN];void NTT(int A[],int nn,int ty,int M){    int i,j,t1,t2,k,m;    for(i=0;i<nn;i++)    {        for(j=0,k=i,m=1;m<nn;m<<=1,j=(j<<1)|(k&1),k>>=1);        if(i<j)        {            t1=A[i];            A[i]=A[j];            A[j]=t1;        }    }    W[0]=1;    for(m=1;m<nn;m<<=1)    {        t1=exp(G,M-1+ty*(M-1)/(m<<1),M);        for(i=1;i<m;i++)        {            W[i]=1ll*W[i-1]*t1%M;        }        for(k=0;k<nn;k+=m<<1)        {            for(i=k;i<k+m;i++)            {                t1=A[i];                t2=1ll*A[i+m]*W[i-k]%M;                A[i]=t1+t2;                A[i]-=A[i]>=M?M:0;                A[i+m]=t1-t2;                A[i+m]+=A[i+m]<0?M:0;            }        }    }    if(ty==1)    {        return ;    }    t1=exp(nn,M-2,M);    for(i=0;i<nn;i++)    {        A[i]=1ll*A[i]*t1%M;    }    return ;}LL mul(LL a,LL b){    LL an=0;    for(;b;b>>=1)    {        if(b&1)        {            an=an+a;            an-=an>=PQ?PQ:0;        }        a=a+a;        a-=a>=PQ?PQ:0;    }    return an;}int T1[MAXN],T2[MAXN],tmp[MAXN];void inv(int A[],int X[],int nn){    int i,j,k,t,m,n;    LL t5;    X[0]=exp(A[0],mod-2,mod);    for(m=1;m<nn;m<<=1)    {        n=m<<1;        for(i=0;i<n;i++)        {            T1[i]=A[i];            T2[i]=A[i];        }        for(n<<=1;i<n;i++)        {            T1[i]=0;            T2[i]=0;        }        for(i=0;i<n;i++)        {            tmp[i]=X[i];        }        NTT(T1,n,1,P);        NTT(tmp,n,1,P);        for(i=0;i<n;i++)        {            T1[i]=2+P-1ll*T1[i]*tmp[i]%P;            T1[i]-=T1[i]>=P?P:0;            T1[i]+=T1[i]<0?P:0;            tmp[i]=1ll*tmp[i]*T1[i]%P;        }        NTT(tmp,n,-1,P);        NTT(T2,n,1,Q);        NTT(X,n,1,Q);        for(i=0;i<n;i++)        {            T2[i]=2+Q-1ll*T2[i]*X[i]%Q;            T2[i]-=T2[i]>=Q?Q:0;            T2[i]+=T2[i]<0?Q:0;            X[i]=1ll*X[i]*T2[i]%Q;        }        NTT(X,n,-1,Q);        for(i=0;i<n;i++)        {            t5=(mul(mul(1ll*X[i],1ll*P),1ll*ip)+mul(mul(1ll*tmp[i],1ll*Q),1ll*iq))%PQ;            if(t5>lim)            {                X[i]=t5%mod-PQ%mod;                X[i]+=X[i]<0?mod:0;            }            else            {                X[i]=t5%mod;            }        }        for(i=n>>1;i<n;i++)        {            X[i]=0;        }    }    for(i=nn;i<n;i++)    {        X[i]=0;    }    for(i=0;i<n;i++)    {        tmp[i]=T1[i]=T2[i]=0;    }    return ;}       void work2(){    A[0]=1;    for(int i=1;i<=n;i++)    {        A[i]%=mod;        A[i]=mod-A[i];    }    inv(A,F,n+1);    printf("%d\n",(F[n]%mod+mod)%mod);    return ;}int main(){    //freopen("shell.in","r",stdin);    //freopen("shell.out","w",stdout);    _r(n);    for(int i=1;i<=n;i++)    {        _r(A[i]);    }    if(n<=1100)    {        work1();    }    else    {        work2();    }    return 0;}

上面的代码很乱，是考场上写的，有很多冗余的变量，将就了看吧.
对了,忘记说了，P取得相当大，用中国剩余定理的时候，乘法要用快速乘，不然会超long long.

---

小朋友与二叉树

题目描述参见bzoj 3625 / Codeforces Round #250

题解

这道题是一个计数问题，容易想到dp的做法.
设f[n]表示点权和为n的二叉树有多少种.考虑n的组成：
左子树点权和+右子树点权和+根节点权值=n.
设权值为k的数给定序列中出现次数为c[k],则根据分步计数原理：

f[n]=∑i+j+k=nf[i]⋅f[j]⋅c[k]

构造多项式：

F(x)=∑i=0+∞f[i]⋅xi      C(x)=∑i=0+∞c[i]⋅xi

则有:

F(x)=F(x)∗F(x)∗C(x)

???有诈!!!
我们在递推式中默认F(0)=f[0]=1.这样才可以包含子树为空的情形.但是C(0)=c[0]=0.
于是F(0)=F(0)*F(0)*C(0)=C(0)=0,推出矛盾！
所以正确的式子：

F(x)=F(x)∗F(x)∗C(x)+1

解得:

F(x)=1±1−4C(x)−−−−−−−−√2C(x)

这种方程真的能这么解吗？其实只需要知道，把这个解代入原方程，方程肯定成立.我能解释的也只有这么多了.
注意到,分母没有常数项，分子要整除分母，则分子也不能有常数项，即±应取负号.
为了简单，再做一下变形：

F(x)=1±1−4C(x)−−−−−−−−√2C(x)=21+1−4C(x)−−−−−−−−√

用多项式开方即可.

参考代码

算了，这题我就不给代码了！本人2017年3月1日在bzoj提交了这道题，并在当时莫名其妙地超过了各位大牛,拿到了rank 1.如果放出代码，rank 1岂不是保不住了！！！

---

bzoj 4555

第二类斯特林数

这个题目做法很多,NTT+分治或者多项式求逆都可以做.
只要式子化得好,直接上裸的NTT也是可以做的.
网上题解很多，这里就不多说了.

参考代码(不用求逆或分治的常数最小的版本)

这个做法的数学推导比较妙，建议先学习一下第二类斯特林数的展开式.

#include<cstdio>#define MAXN 410000#define P 998244353#define g 3int n,A[MAXN],B[MAXN],N;int exp(int a,int k,int mod=P){    int an=1;    for(;k;k>>=1)    {        if(k&1)        {            an=1ll*an*a%mod;        }        a=1ll*a*a%mod;    }    return an;}int fac[MAXN],inv[MAXN],W[MAXN];void NTT(int A[],int nn,int ty){    int i,j,k,m,t1,t2;    for(i=0;i<nn;i++)    {        for(j=0,k=i,m=1;m<nn;m<<=1,j=(j<<1)|(k&1),k>>=1);        if(i<j)        {            t1=A[i];            A[i]=A[j];            A[j]=t1;        }    }    W[0]=1;    for(m=1;m<nn;m<<=1)    {        t1=exp(g,P-1+ty*(P-1)/(m<<1),P);        for(i=1;i<m;i++)        {            W[i]=1ll*W[i-1]*t1%P;        }        for(k=0;k<nn;k+=m<<1)        {            for(i=k;i<k+m;i++)            {                t1=A[i];                t2=1ll*A[i+m]*W[i-k]%P;                A[i]=t1+t2;                A[i]-=A[i]>=P?P:0;                A[i+m]=t1-t2;                A[i+m]+=A[i+m]<0?P:0;            }        }    }    if(ty==1)    {        return ;    }    t1=exp(nn,P-2,P);    for(i=0;i<nn;i++)    {        A[i]=1ll*A[i]*t1%P;    }    return ;}int main(){    scanf("%d",&n);    fac[0]=inv[0]=1;    for(int i=1;i<=n;i++)    {        inv[i]=1ll*inv[i-1]*i%P;        fac[i]=inv[i];    }    inv[n]=exp(inv[n],P-2,P);    for(int i=n;i;i--)    {        inv[i-1]=1ll*inv[i]*i%P;    }    for(int i=0,p;i<=n;i++)    {        if(i&1)        {            A[i]=P-inv[i];        }        else        {            A[i]=inv[i];        }        if(i==0)        {            B[i]=1;        }        else if(i==1)        {            B[i]=1ll*(n+1)*inv[i]%P;        }        else        {            B[i]=1ll*(exp(i,n+1)-1+P)%P*exp(i-1,P-2,P)%P*inv[i]%P;        }    }    for(N=1;N<=n*2;N<<=1);    NTT(A,N,1);    NTT(B,N,1);    for(int i=0;i<N;i++)    {        A[i]=1ll*A[i]*B[i]%P;    }    NTT(A,N,-1);    int ans=0,tmp;    for(int i=0;i<=n;i++)    {        tmp=1ll*exp(2,i,P)*fac[i]%P;        tmp=1ll*tmp*A[i]%P;        ans=ans+tmp;        ans-=ans>=P?P:0;    }    printf("%d\n",ans);    return 0;}

---

【SDOI2013 R1 Day2】淘金

题目来源：sdoi 2013

题号：bzoj 3131

题目描述

小Z在玩一个叫做《淘金者》的游戏。游戏的世界是一个二维坐标。X轴、Y轴坐标范围均为1..N。初始的时候，所有的整数坐标点上均有一块金子，共N*N块。
一阵风吹过，金子的位置发生了一些变化。细心的小Z发现，初始在(i,j)坐标处的金子会变到(f(i),f(j))坐标处。其中f(x)表示x各位数字的乘积，例如f(99)=81，f(12)=2，f(10)=0。如果金子变化后的坐标不在1..N的范围内，我们认为这块金子已经被移出游戏。同时可以发现，对于变化之后的游戏局面，某些坐标上的金子数量可能不止一块，而另外一些坐标上可能已经没有金子。这次变化之后，游戏将不会再对金子的位置和数量进行改变，玩家可以开始进行采集工作。
小Z很懒，打算只进行K次采集。每次采集可以得到某一个坐标上的所有金子，采集之后，该坐标上的金子数变为0。
现在小Z希望知道，对于变化之后的游戏局面，在采集次数为K的前提下，最多可以采集到多少块金子？
答案可能很大，小Z希望得到对1000000007(10^9+7)取模之后的答案。

数据规模

N≤1012,k≤105.

题解

一些简单的想法：

• 处理出每个位置的金子数
• 如果设g(x)表示所有满足f(i)=x的i的数目，那么(x,y)的金子数为g(x)*g(y)
• 如果我们处理出了g(x)的值，那么这就是一个经典的用堆处理的问题了.
• 处理g(x)应该可以用数位dp.

官方题解：数位dp+堆.

另一些想法：

• 是不是可以暴力求解g(x)?
• g(x)!=0,意味着x=2a∗3b∗5c∗7d，这样的x并不多！
• 如果N≤106,我们暴力求解即可
• N≤1012,则N≤106∗106我们暴力处理出106以内的g(x)，然后用NTT做一次乘法就可以得出1012以内的g(x)了！
• 显然x可以大到912,不过我们可以分解质因数用a,b,c,d存储.为了方便起见，估计a,b,c,d的最大值，然后采取合适的进制就可以了.
• 于是这道题居然就可以用NTT+暴力就水过去了！
• 当然这个乘法并不简单，考虑到数位dp的特点，我们在做乘法时要分有限制和无限制两部分,具体实现看代码吧

参考代码

丑的要命的代码：

#include<cstdio>#include<queue>#include<cstring>#include<iostream>#include<algorithm>using namespace std;#define MAXN 1080000#define mod 1000000007#define LL long long#define G 998244353#define K 1004535809#define gg 3int exp(int a,int k){    int an=1;    for(;k;k>>=1)    {        if(k&1)        {            an=1ll*an*a%G;        }        a=1ll*a*a%G;    }    return an;}int W[MAXN];void NTT(int A[],int nn,int ty){    int i,j,t1,t2,k,m;    for(i=0;i<nn;i++)    {        for(j=0,k=i,m=1;m<nn;m<<=1,j=(j<<1)|(k&1),k>>=1);        if(i<j)        {            t1=A[i];            A[i]=A[j];            A[j]=t1;        }    }    W[0]=1;    for(m=1;m<nn;m<<=1)    {        t1=exp(gg,G-1+ty*(G-1)/(m<<1));        for(i=1;i<m;i++)        {            W[i]=1ll*W[i-1]*t1%G;        }        for(k=0;k<nn;k+=m<<1)        {            for(i=k;i<k+m;i++)            {                t1=A[i];                t2=1ll*A[i+m]*W[i-k]%G;                A[i]=t1+t2;                A[i]-=A[i]>=G?G:0;                A[i+m]=t1-t2;                A[i+m]+=A[i+m]<0?G:0;            }        }    }    if(ty==1)    {        return ;    }    t1=exp(nn,G-2);    for(i=0;i<nn;i++)    {        A[i]=1ll*A[i]*t1%G;    }    return ;}LL n,k;int f[MAXN],g[MAXN],tot=0,h[MAXN],oto=0;int ans=0;void work1(){    for(int i=1,tt;i<=n;i++)    {        tt=1;        for(int x=i;x;x/=10)        {            tt*=x%10;        }        f[tt]++;    }    for(int i=1;i<=n;i++)    {        if(f[i])        {            g[++tot]=f[i];        }    }    sort(g+1,g+1+tot);    for(int i=tot;i>=1&&i>=tot-1000;i--)    {        for(int j=tot;j>=1&&j>=tot-1000;j--)        {            h[++oto]=1ll*g[i]*g[j]%mod;        }    }    sort(h+1,h+1+oto);    for(int i=oto;i>=1&&i>=oto-k+1;i--)    {        ans+=h[i];        ans-=ans>=mod?mod:0;    }    printf("%d\n",ans);    return ;}//分段，前面小数据直接暴力！const int QQ=1000000;LL n1,n2,m1,m2,AA[MAXN*3];int pp[5]={0,2,3,5,7},qq[5]={0,40,26,18,15};int A[MAXN],B[MAXN],C[MAXN],D[MAXN],NN=1048576;struct node{    int x,y;    LL val;    node(int xx=0,int yy=0)    {        x=xx;        y=yy;        val=1ll*D[x]*D[y];    }};bool cmp(LL a,LL b){    return a>b;}bool operator < (node x,node y){    return x.val<y.val;}priority_queue<node>Q;void work2(){    n1=n/QQ;    n2=n1-1;    m1=n%QQ;    m2=QQ-1;    memset(f,0,sizeof(f));     for(int i=QQ/10,tt;i<=m2;i++)    {        tt=1;        for(int x=i;x;x/=10)        {            tt*=x%10;        }        f[tt]++;    }    tot=0;    for(int i=1;i<=m2;i++)    {        if(f[i])        {            int tt=0;            for(int j=4,q=i;j>=1;j--)            {                tt*=qq[j];                while(q%pp[j]==0)                {                    q/=pp[j];                    ++tt;                }            }            D[tt]=f[i];        }    }    memset(f,0,sizeof(f));     for(int i=QQ/10,tt;i<=m1;i++)    {        tt=1;        for(int x=i;x;x/=10)        {            tt*=x%10;        }        f[tt]++;    }    tot=0;    for(int i=1;i<=m1;i++)    {        if(f[i])        {            int tt=0;            for(int j=4,q=i;j>=1;j--)            {                tt*=qq[j];                while(q%pp[j]==0)                {                    q/=pp[j];                    ++tt;                }            }            C[tt]=f[i];        }    }    memset(f,0,sizeof(f));     for(int i=1,tt;i<=n2;i++)    {        tt=1;        for(int x=i;x;x/=10)        {            tt*=x%10;        }        f[tt]++;    }    tot=0;    for(int i=1;i<=n2;i++)    {        if(f[i])        {            int tt=0;            for(int j=4,q=i;j>=1;j--)            {                tt*=qq[j];                while(q%pp[j]==0)                {                    q/=pp[j];                    ++tt;                }            }            B[tt]=f[i];        }    }    int ttt=1;    for(int x=n1;x;x/=10)    {        ttt*=x%10;    }    int tt=0;    if(ttt)    {        for(int j=4,q=ttt;j>=1;j--)        {            tt*=qq[j];            while(q%pp[j]==0)            {                q/=pp[j];                ++tt;            }        }    }    A[tt]=1;    NTT(A,NN,1);    NTT(C,NN,1);    for(int i=0;i<NN;i++)    {        A[i]=1ll*A[i]*C[i]%G;    }    NTT(A,NN,-1);    NTT(B,NN,1);    NTT(D,NN,1);    for(int i=0;i<NN;i++)    {        B[i]=1ll*B[i]*D[i]%G;    }    NTT(B,NN,-1);    for(int i=0;i<NN;i++)    {        C[i]=A[i]+B[i];        C[i]-=C[i]>=mod?mod:0;    }    memset(f,0,sizeof(f));     memset(D,0,sizeof(D));    for(int i=1,tt;i<=m2;i++)    {        tt=1;        for(int x=i;x;x/=10)        {            tt*=x%10;        }        f[tt]++;    }    tot=0;    for(int i=1;i<=m2;i++)    {        if(f[i])        {            int tt=0;            for(int j=4,q=i;j>=1;j--)            {                tt*=qq[j];                while(q%pp[j]==0)                {                    q/=pp[j];                    ++tt;                }            }            D[tt]=f[i];        }    }    for(int i=0;i<NN;i++)    {        C[i]+=D[i];        C[i]-=C[i]>=mod?mod:0;    }    int tttt=0;    for(int i=0;i<NN;i++)    {        if(C[i])        {            D[++tttt]=C[i];        }    }    sort(D+1,D+1+tttt,cmp);    Q.push(node(1,1));    int cnt=0;    node q;    LL ans=0;    while(Q.size())    {        q=Q.top();        Q.pop();        ans+=q.val%mod;        ans-=ans>=mod?mod:0;        ++cnt;        if(cnt==k)        {            break;        }        if(q.x!=q.y)        {            ans+=q.val%mod;            ans-=ans>=mod?mod:0;            ++cnt;            if(cnt==k)            {                break;            }            if(q.x+1<=tttt)            {                Q.push(node(q.x+1,q.y));            }        }        if(q.x==1&&q.y+1<=tttt)        {            Q.push(node(q.x,q.y+1));        }    }    printf("%I64d\n",ans);    return ;}int main(){    //freopen("gold.in","r",stdin);    //freopen("gold.out","w",stdout);    cin>>n>>k;    if(n<=2000000)    {        work1();    }    else    {        work2();    }    return 0;}

代码中的A,B,C,D数组是什么鬼？t,tt,ttt,tttt记的是什么？这些嘛，自己看自己想，反正就是乱整就对了.

至于为什么NTT不会出现结果超出998244353的情况，我只能说，打表发现，A,B,C,D这些东西并不是很大.

这份代码在bzoj上光荣拿到了倒数第2(用户:problem_set)，想想吧，人家数位dp自然是比这个暴力+NTT要快很多的