HDU5519 Kykneion asma (指数生成函数+快速数论变换模任意数+启发式合并思想)

原文链接：HDU5519 Kykneion asma (指数生成函数+快速数论变换模任意数+启发式合并思想)

先说一下，这个不是正解。但是也可以过。
题意：有5个数字——0,1,2,3,4，每个数字分别有a0,a1,a2,a3,a4个。问这些数字能组成多少个n位数？
数据范围：a<=30000,n<=15000
时限：6s
分析：
首先n位数肯定是排列，每种数字有很多个，就是多重集，这个就是多重集的排列问题。显然是指数生成函数可以做的。
关于指数生成函数可以看看我前面的生成函数这个课件。
现在有个问题，就是怎么消除前导0？
设W(n,a0,a1,a2,a3,a4) 代表可以有前导0的n位5进制数字，至多ai个数字为i的方案数，那么答案等于W(n,a0,a1,a2,a3,a4)−W(n−1,a0−1,a1,a2,a3,a4)，相当于减去一定有前导0的方案数。
这下问题转换成了求多项式第n次项的系数乘以n的阶乘.
这个显然直接ntt就完了。
但是这里有个问题，就是模的数不是费马质数，这个显然直接上ntt不能做。
所以我们可以嘿嘿嘿
找3个1e9左右的费马质数，比如998244353,1005060097,950009857，然后对每个质数都做一次ntt，用CRT合并解。
注意到CRT的过程要把模数乘起来，这样long long不是很资瓷，于是需要手写一个高精度或者int128.
这个题我学会了一个东西就是：在写ntt模任意数的时候，多项式必须两两相乘，因为CRT这里限制了不能把多项式全部ntt然后乘起来变回去。
所以这里还有一个思想，就是启发式合并。对于t个多项式，我们开一个集合维护它们的长度。每次选取最短的两个进行ntt，然后将得到的新多项式塞集合。这样显然比较快。
对于本题，因为t只有5，比较小，所以启发式合并就不需要了，但是也不应该用一个数组保存解，一直乘下去，这样明显不优。应该两个两个地乘，再把每一对的乘积拿来乘。
下面是代码

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define ll long long
using namespace std;
const ll M=131073ll,MOD=1000000007ll,P[3]={998244353ll,1005060097ll,950009857ll},G[3]={3ll,5ll,7ll},Inv[3]={644348675ll,675933219ll,647895261ll};
int n,m,cas,a0,a1,a2,a3,a4;
ll inv[M],fac[M];
inline ll pow(ll x,ll y,ll P)
{
    ll re=1;
    while(y)
    {
        if(y&1)re=re*x%P;
        x=x*x%P;
        y>>=1;
    }
    return re;
}
struct Int_128{
    unsigned long long a,b;
    Int_128(ll x){a=0,b=x;}
    friend bool operator < (Int_128 x,Int_128 y)
    {
        return x.a<y.a||x.a==y.a&&x.b<y.b;
    }
    friend Int_128 operator + (Int_128 x,Int_128 y)
    {
        Int_128 re(0);
        re.a=x.a+y.a+(x.b+y.b<x.b);
        re.b=x.b+y.b;
        return re;
    }
    friend Int_128 operator - (Int_128 x,Int_128 y)
    {
        y.a=~y.a;y.b=~y.b;
        return x+y+1;
    }
    void Div2()
    {
        b>>=1;b|=(a&1ll)<<63;a>>=1;
    }
    friend Int_128 operator * (Int_128 x,Int_128 y)
    {
        Int_128 re=0;
        while(y.a||y.b)
        {
            if(y.b&1)re=re+x;
            x=x+x;y.Div2();
        }
        return re;
    }
    friend Int_128 operator % (Int_128 x,Int_128 y)
    {
        Int_128 temp=y;
        int cnt=0;
        while(temp<x)temp=temp+temp,++cnt;
        for(;cnt>=0;cnt--)
        {
            if(temp<x)x=x-temp;
            temp.Div2();
        }
        return x;
    }
};
void ntt(ll*a,int n,int f,int P,int G)
{
    for(int i=1,j=0;i<n-1;++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=pow(G,(P-1)/(i<<1),P);
        for(int j=0;j<n;j+=i<<1)
        {
            ll w=1;
            for(int k=0;k<i;++k,w=w*wn%P)
            {
                ll x=a[j+k],y=w*a[j+k+i]%P;
                a[j+k]=(x+y)%P;
                a[j+k+i]=(x-y+P)%P;
            }
        }
    }
    if(f==-1)
    {
        for(int i=1;i<(n>>1);++i)swap(a[i],a[n-i]);
        ll inv=pow(n,P-2,P);
        for(int i=0;i<n;++i)a[i]=a[i]*inv%P;
    }
}
inline void Polynomial_Multiplication(ll a[],ll b[],ll c[],int n)
{
    static ll A[3][M],B[3][M];
    for(int j=0;j<3;++j)
    {
        for(int i=0;i<n;++i)
        {
            A[j][i]=a[i]%P[j];
            B[j][i]=b[i]%P[j];
        }
        ntt(A[j],n,1,P[j],G[j]);
        ntt(B[j],n,1,P[j],G[j]);
        for(int i=0;i<n;i++)A[j][i]=A[j][i]*B[j][i]%P[j];
        ntt(A[j],n,-1,P[j],G[j]);
    }
    Int_128 _MOD=Int_128(P[0]*P[1])*P[2];
    for(int i=0;i<n;++i)
    {
        Int_128 temp=
        Int_128(P[1]*P[2])*Int_128(Inv[0]*A[0][i])+
        Int_128(P[0]*P[2])*Int_128(Inv[1]*A[1][i])+
        Int_128(P[0]*P[1])*Int_128(Inv[2]*A[2][i]);
        c[i]=(temp%_MOD%MOD).b;
    }
}
int main()
{
    scanf("%d",&cas);
    inv[0]=inv[1]=fac[0]=fac[1]=1;
    for(int i=2;i<=15000;++i)inv[i]=(MOD-MOD/i)*inv[MOD%i]%MOD,fac[i]=fac[i-1]*i%MOD;
    for(int i=2;i<=15000;++i)inv[i]=inv[i]*inv[i-1]%MOD;
    static ll A[M],B[M],C[M],D[M],E[M],F[M];
    for(int t=1;t<=cas;++t)
    {
        scanf("%d%d%d%d%d%d",&n,&a0,&a1,&a2,&a3,&a4);
        for(int i=0;i<=min(n,a0);++i)A[i]=inv[i];
        for(int i=min(n,a0)+1;i<=m;++i)A[i]=0;
        for(int i=0;i<=min(n,a1);++i)B[i]=inv[i];
        for(int i=min(n,a1)+1;i<=m;++i)B[i]=0;
        for(int i=0;i<=min(n,a2);++i)C[i]=inv[i];
        for(int i=min(n,a2)+1;i<=m;++i)C[i]=0;
        for(int i=0;i<=min(n,a3);++i)D[i]=inv[i];
        for(int i=min(n,a3)+1;i<=m;++i)D[i]=0;
        for(int i=0;i<=min(n,a4);++i)E[i]=inv[i];
        for(int i=min(n,a4)+1;i<=m;++i)E[i]=0;
        /*运用了启发式合并的思想，跑了5.1s*/
        for(m=1;m<min(a1,n)+min(a2,n);m<<=1);
        Polynomial_Multiplication(B,C,B,m);
        for(m=1;m<min(a3,n)+min(a4,n);m<<=1);
        Polynomial_Multiplication(D,E,D,m);
        for(m=1;m<min(a1,n)+min(a2,n)+min(a3,n)+min(a4,n);m<<=1);
        Polynomial_Multiplication(B,D,B,m);

    /*
        如果不用启发式合并，直接B数组乘到底，5.8s，已经是卡着过了
        for(m=1;m<min(a1,n)+min(a2,n);m<<=1);
        Polynomial_Multiplication(B,C,B,m);
        for(m=1;m<min(a1,n)+min(a2,n)+min(a3,n);m<<=1);
        Polynomial_Multiplication(B,D,B,m);
        for(m=1;m<min(a1,n)+min(a2,n)+min(a3,n)+min(a4,n);m<<=1);
        Polynomial_Multiplication(B,E,B,m);
    */  for(m=1;m<min(a0,n)+min(a1,n)+min(a2,n)+min(a3,n)+min(a4,n);m<<=1);
        if(!a0)
        {
            Polynomial_Multiplication(A,B,F,m);
            printf("Case #%d: %I64d\n",t,F[n]*fac[n]%MOD);
        }
        else
        {
            Polynomial_Multiplication(A,B,F,m);
            ll ans=F[n]*fac[n]%MOD;
            A[min(a0,n)]=0;
            Polynomial_Multiplication(A,B,F,m);
            printf("Case #%d: %I64d\n",t,((ans-F[n-1]*fac[n-1]%MOD)%MOD+MOD)%MOD);
        }
    }
}