【总结】MTT算法
前言
51nod1258真是道好题。。。
一道题，学会了3个东西：伯努利数，自然数幂和，MTT…
前置科技（其实学MTT的人估计都会。。。）
CRT(中国剩余定理)
NTT
NTT算法的局限
众所周知，NTT是通过原根的性质来进行快速傅里叶变化
不过，同时也对其模数做出了要求，对于一个模数M，若 φ ( M ) \varphi(M) φ(M)中，2的次数较少（小于要转的序列长度）就不可做了
所以，也因此衍生出了“NTT模数”，即 φ ( M ) \varphi(M) φ(M)中2的次数较多的模数。
常见的有： 998244353 = 2 23 ∗ 119 + 1 998244353=2^{23}*119+1 998244353=223∗119+1， 1004535809 = 2 21 ∗ 479 + 1 1004535809=2^{21}*479+1 1004535809=221∗479+1， 469762049 = 2 26 ∗ 7 + 1 469762049=2^{26}*7+1 469762049=226∗7+1
（请至少记住这三个模数。。。因为后面要用。。。）
但是，有些无良的出题人，由于懒得想好题，于是就把一些比较简单的多项式题目中的模数换掉（如换成1e9+7），以此来强迫选手要么写拆系数FFT，要么写MTT…
拆系数FFT容易卡精然后WA
MTT容易卡常然后T
To wa or to tle,that’s a question.
MTT
MTT说白了就是用三个NTT模数分别算出一个解，然后再用CRT合并出原来的解
为什么要用三个模数呢。。这是因为，根据CRT的使用条件，如果真正的答案比你的模数之积大，那是得不到唯一解的。所以要求你的模数之积超过真正的解得值域上界。然而一般来说，用MTT都是因为模数是1e9+7，然后 ( 1 e 9 + 7 ) 2 ∗ N ( 1 0 5 左 右 ) ≈ 1 0 23 (1e9+7)^2*N(10^5左右)≈10^{23} (1e9+7)2∗N(105左右)≈1023，所以为了超过这个值，一般都用3个模数。
不过很多情况下这是会爆long long的（不爆long long还不如直接去FFT然后强转。。。）
所以，不能直接利用CRT。
但可以通过一些技巧，算出真正的解
已知
a
n
s
≡
c
0
(
m
o
d
m
0
)
ans\equiv c_0(mod \ m_0)
ans≡c0(mod m0)
a
n
s
≡
c
1
(
m
o
d
m
1
)
ans\equiv c_1(mod \ m_1)
ans≡c1(mod m1)
a
n
s
≡
c
2
(
m
o
d
m
2
)
ans\equiv c_2(mod \ m_2)
ans≡c2(mod m2)
首先利用CRT合并前两项
记
I
m
0
=
m
0
−
1
(
m
o
d
m
1
)
,
I
m
1
=
m
1
−
1
(
m
o
d
m
0
)
Im_0=m_0^{-1}(mod \ m_1),Im_1=m_1^{-1}(mod\ m_0)
Im0=m0−1(mod m1),Im1=m1−1(mod m0)
a
n
s
=
(
c
0
∗
m
1
∗
I
m
1
+
c
1
∗
m
0
∗
I
m
0
)
(
m
o
d
m
0
m
1
)
ans=(c_0*m_1*Im_1+c_1*m_0*Im_0)(mod\ m_0m_1)
ans=(c0∗m1∗Im1+c1∗m0∗Im0)(mod m0m1)
注意这一步需要快速乘避免爆long long
所以设
a
n
s
=
x
∗
M
+
C
(
M
=
m
0
m
1
)
ans= x*M+C(M=m_0m_1)
ans=x∗M+C(M=m0m1)
然后带入到
a
n
s
≡
c
2
(
m
o
d
m
2
)
ans\equiv c_2(mod\ m_2)
ans≡c2(mod m2)中，得
x
≡
c
2
−
C
M
(
m
o
d
m
2
)
x\equiv \frac {c_2-C} {M}(mod\ m_2)
x≡Mc2−C(mod m2)
所以将
c
2
−
C
M
\frac {c_2-C} {M}
Mc2−C设为q
得
x
=
q
+
k
m
2
x=q+km_2
x=q+km2，再代回
a
n
s
=
x
∗
M
+
C
ans= x*M+C
ans=x∗M+C中，有
a
n
s
=
k
m
0
m
1
m
2
+
q
M
+
C
ans=km_0m_1m_2+qM+C
ans=km0m1m2+qM+C
由于
a
n
s
<
m
0
m
1
m
2
ans< m_0m_1m_2
ans<m0m1m2所以
k
=
0
k=0
k=0
所以得到
a
n
s
=
q
M
+
C
ans=qM+C
ans=qM+C。。。
#include<cstdio>
#include<algorithm>
#include<vector>
#define SF scanf
#define PF printf
#define MAXN 300010
#define MOD 1000000007
using namespace std;
typedef long long ll;
const int G=3;
ll g[MAXN],B[MAXN],ifac[MAXN],fac[MAXN];
ll fsa(ll x,ll y,ll mod){
// x%=mod,y%=mod;
// ll res=0;
// while(y){
// if(y&1ll)
// res=(res+x)%mod;
// x=(x+x)%mod;
// y>>=1ll;
// }
// return res;
return (x*y-(ll)((long double)x/mod*y)*mod+mod)%mod;
}
ll fsp(ll x,ll y,ll mod){
x%=mod;
ll res=1;
while(y){
if(y&1ll)
res=res*x%mod;
x=x*x%mod;
y>>=1ll;
}
return res;
}
void NTT(ll A[],int N,int flag,ll mod){
for(int i=1,j=0;i<N;i++){
for(int d=N;j^=d>>=1,~j&d;);
if(i<j)
swap(A[i],A[j]);
}
for(int i=1;i<N;i<<=1){
ll wn=fsp(G,(mod-1)/(i<<1),mod);
if(flag)
wn=fsp(wn,mod-2,mod);
for(int j=0;j<N;j+=(i<<1)){
ll w=1;
for(int k=0;k<i;k++,w=w*wn%mod){
ll x=A[j+k],y=A[i+j+k]*w%mod;
A[j+k]=(x+y)%mod;
A[i+j+k]=(x-y+mod)%mod;
}
}
}
if(flag){
ll invN=fsp(N,mod-2,mod);
for(int i=0;i<N;i++)
A[i]=A[i]*invN%mod;
}
}
void Mul(ll A[],ll B[],int N,int M,ll res[],ll mod){
static ll tmp1[MAXN],tmp2[MAXN];
for(int i=0;i<N;i++)
tmp1[i]=A[i];
for(int i=0;i<M;i++)
tmp2[i]=B[i];
int p=1;
while(p<=N+M)
p<<=1;
NTT(tmp1,p,0,mod);
NTT(tmp2,p,0,mod);
for(int i=0;i<p;i++)
res[i]=tmp1[i]*tmp2[i]%mod;
NTT(res,p,1,mod);
for(int i=0;i<p;i++)
tmp1[i]=tmp2[i]=0;
}
ll Mod[3]={998244353,469762049,1004535809};
void MTT(ll A[],ll B[],int N,int M,ll res[]){
static ll res1[MAXN],res2[MAXN],res3[MAXN];
Mul(A,B,N,M,res1,Mod[0]);
Mul(A,B,N,M,res2,Mod[1]);
Mul(A,B,N,M,res3,Mod[2]);
ll M1=Mod[0]*Mod[1];
ll IM0=fsp(Mod[0],Mod[1]-2,Mod[1]);
ll IM1=fsp(Mod[1],Mod[0]-2,Mod[0]);
for(int i=0;i<N+M;i++){
ll c=(fsa(res1[i]*Mod[1],IM1,M1)+fsa(res2[i]*Mod[0],IM0,M1))%M1;
ll x=((res3[i]-c)%Mod[2]+Mod[2])%Mod[2]*fsp(M1,Mod[2]-2,Mod[2])%Mod[2];
res[i]=(x*(M1%MOD)%MOD+c)%MOD;
}
for(int i=0;i<N+M;i++)
res1[i]=res2[i]=res3[i]=0;
}
void Inv(ll A[],ll B[],int N){
if(N==1){
B[0]=fsp(A[0],MOD-2,MOD);
return ;
}
Inv(A,B,(N+1)>>1);
static ll tmp1[MAXN],tmp2[MAXN];
MTT(A,B,N,N,tmp1),MTT(B,tmp1,N,N,tmp2);
// if(N==8){
// for(int i=0;i<N;i++)
// PF("[%lld]",B[i]);
// PF("\n");
// for(int i=0;i<N;i++)
// PF("[%lld]",tmp1[i]);
// PF("\n");
// for(int i=0;i<N;i++)
// PF("[%lld]",tmp2[i]);
// PF("\n");
// }
for(int i=0;i<N;i++) B[i]=(B[i]*2ll-tmp2[i]+MOD)%MOD;
for(int i=0;i<2*N;i++) tmp1[i]=tmp2[i]=0;
}
ll C(int n,int m){
return fac[n]*ifac[m]%MOD*ifac[n-m]%MOD;
}
int main(){
int n=65536;
fac[0]=1;
for(int i=1;i<=n*2;i++)
fac[i]=fac[i-1]*i%MOD;
ifac[n*2]=fsp(fac[n*2],MOD-2,MOD);
for(int i=n*2;i>=1;i--)
ifac[i-1]=ifac[i]*i%MOD;
for(int i=0;i<n;i++)
g[i]=ifac[i+1];
// for(int i=0;i<8;i++)
// PF("%lld ",g[i]);
Inv(g,B,n);
for(int i=0;i<n;i++)
B[i]=B[i]*fac[i]%MOD;
// for(int i=0;i<8;i++)
// PF("%lld ",B[i]);
int T,k;
ll N;
SF("%d",&T);
while(T--){
SF("%lld%d",&N,&k);
ll ans=0;
ll nw=(N+1)%MOD,q=(N+1)%MOD;
for(int i=k;i>=0;i--,nw=nw*q%MOD)
ans=(ans+C(k+1,i)*B[i]%MOD*nw%MOD)%MOD;
ans=ans*fsp(k+1,MOD-2,MOD)%MOD;
PF("%lld\n",ans);
}
}