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

来源：互联网发布：java acos函数编辑：程序博客网时间：2024/06/08 00:46

原文链接：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 longusing 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);        }    }}

0 0