51nod 1258：序列求和 V4（伯努利数+多项式求逆+NTT+CRT）

题面
题意：求$S_k(n)$
n≤1e18，k≤5e4。

用多项式逆元求出伯努利数，代公式就可以了。
上维基查了下伯努利数，$B^+$与$B^-$的差别只是$B_1$的正负而已.

但我依然不知道什么时候用哪个，两个都试一下，发现这题用$B^+$。

所以我做这题是为了学CRT的。

假设题意要求模P，三模数为p1,p2,p3。

我起初还天真地以为可以算出分别模p1,p2,p3的伯努利数，然后合并
以至于我三模数没有开数组，代码量*3。

根据我粗鄙地理解，若x不是模P意义下的，则它的实际值不能太大。
比如做A*B*C的时候应先做出A*B%P的系数表达式，再乘C。
和拆系数FFT很类似。

通过这题，我对拆系数FFT更加膜拜了。

#include <iostream>
#include <fstream>
#include <algorithm>
#include <cmath>
#include <ctime>
#include <cstdio>
#include <cstdlib>
#include <cstring>

using namespace std;
#define mmst(a, b) memset(a, b, sizeof(a))
#define mmcp(a, b) memcpy(a, b, sizeof(b))

typedef long long LL;

const int N=200200;
const LL p1=998244353,p2=1004535809,p3=469762049,P=1e9+7;

int T,k,m=65536,rev[N];
LL n;
LL B1[N],B2[N],B3[N];
LL X1[N],X2[N],X3[N];
LL I[N],jc[N],Ijc[N],cf[N];
LL ans;

LL cheng(LL a,LL b,LL p)
{
    LL res=1;
    for(;b;b>>=1,a=a*a%p)
    if(b&1)
    res=res*a%p;
    return res;
}

LL mul(LL a,LL b,LL p)
{
    a%=p; b%=p;
    return ((a*b-(LL)((LL)((long double)a/p*b+1e-3)*p))%p+p)%p;
}

const LL M=p1*p2,inv1=cheng(p1%p2,p2-2,p2),inv2=cheng(p2%p1,p1-2,p1),inv3=cheng(M%p3,p3-2,p3);

LL CRT(LL a1,LL a2,LL a3)
{
    LL A=(mul(a1*p2%M,inv2,M)+mul(a2*p1%M,inv1,M))%M;
    LL K=(a3+p3-A%p3)*inv3%p3;
    return (K*(M%P)+A)%P;
}

void ntt(LL *a,LL ops,LL p)
{
    for(int i=0;i<n;i++)
    if(i<rev[i])
    swap(a[i],a[rev[i]]);
    for(int l=2,m=1;l<=n;l<<=1,m<<=1)
    {
        LL wn=(ops) ? cheng(3,(p-1)/l,p) : cheng(3,p-1-(p-1)/l,p);
        for(int i=0;i<n;i+=l)
        {
            LL w=1;
            for(int k=0;k<m;k++)
            {
                LL t=a[i+k+m]*w%p;
                a[i+k+m]=(a[i+k]-t+p)%p;
                a[i+k]=(a[i+k]+t)%p;
                w=w*wn%p;
            }
        }
    }
    if(!ops)
    {
        LL inv=cheng(n,p-2,p);
        for(int i=0;i<n;i++)
        a[i]=a[i]*inv%p;
    }
}

void ny(LL *a,int len)
{
    if(len==1)
    return;
    ny(a,len/2);

    for(int i=0;i<len;i++)
    X1[i]=a[i]%p1,X2[i]=a[i]%p2,X3[i]=a[i]%p3;

    n=1;
    int k=-1;
    while(n<=len)
    n<<=1,k++;
    for(int i=0;i<n;i++)
    rev[i]=(rev[i>>1]>>1) | ((i&1)<<k);

    ntt(X1,1,p1);
    ntt(B1,1,p1);
    ntt(X2,1,p2);
    ntt(B2,1,p2);
    ntt(X3,1,p3);
    ntt(B3,1,p3);

    for(int i=0;i<n;i++)
    {
        X1[i]=X1[i]*B1[i]%p1;
        X2[i]=X2[i]*B2[i]%p2;
        X3[i]=X3[i]*B3[i]%p3;
    }

    ntt(X1,0,p1);
    ntt(X2,0,p2);
    ntt(X3,0,p3);

    for(int i=0;i<n;i++)
    {
        X1[i]=CRT(X1[i],X2[i],X3[i]);
        X1[i]=(P-X1[i])%P;
    }

    X1[0]=(2+X1[0])%P;

    for(int i=0;i<n;i++)
    {
        X2[i]=X1[i]%p2;
        X3[i]=X1[i]%p3;
        X1[i]=X1[i]%p1;
    }

    ntt(X1,1,p1);
    ntt(X2,1,p2);
    ntt(X3,1,p3);

    for(int i=0;i<n;i++)
    {
        B1[i]=B1[i]*X1[i]%p1;
        B2[i]=B2[i]*X2[i]%p2;
        B3[i]=B3[i]*X3[i]%p3;
    }

    ntt(B1,0,p1);
    ntt(B2,0,p2);
    ntt(B3,0,p3);

    for(int i=0;i<n;i++)
    B1[i]=B2[i]=B3[i]=CRT(B1[i],B2[i],B3[i]);
    for(int i=len;i<n;i++)
    B1[i]=B2[i]=B3[i]=0;
}

int main()
{
    B1[0]=B2[0]=B3[0]=I[1]=Ijc[0]=jc[0]=1;
    for(int i=2;i<=m;i++)
    I[i]=(P-P/i)*I[P%i]%P;

    for(int i=1;i<m;i++)
    jc[i]=jc[i-1]*i%P,Ijc[i]=Ijc[i-1]*I[i+1]%P;

    ny(Ijc,m);
    for(int i=0;i<m;i++)
    B1[i]=B1[i]*jc[i]%P;
    B1[1]=(P+1)/2;

    cin>>T;
    while(T--)
    {
        ans=0;
        cin>>n>>k;
        n%=P;
        Ijc[0]=cf[0]=1;

        for(int i=1;i<=k+1;i++)
        Ijc[i]=Ijc[i-1]*I[i]%P,cf[i]=cf[i-1]*n%P;

        for(int i=0;i<=k;i++)
        ans=(ans+jc[k+1]*Ijc[i]%P*Ijc[k+1-i]%P*B1[i]%P*cf[k+1-i]%P)%P;
        ans=ans*I[k+1]%P;

        cout<<ans<<endl;
    }

    return 0;
}