zoukankan      html  css  js  c++  java
  • 【BZOJ4555】[TJOI2016&HEOI2016] 求和(NTT)

    点此看题面

    大致题意: 计算(sum_{i=0}^nsum_{j=0}^iS(i,j)*2^j*(j!)),其中(S)为第二类斯特林数。

    推式子

    首先让我们来推一波式子:

    因为当(i<j)时,(S(i,j)=0),所以,为了方便式子的化简,我们可以先将第二个(sum)的上限全部改成(n),即:

    [sum_{i=0}^nsum_{j=0}^nS(i,j)*2^j*(j!) ]

    这样一来,(sum_{j=0}^n2^j*(j!))这个式子就与(i)无关,则我们可以将其提前:

    [sum_{j=0}^n2^j*(j!)sum_{i=0}^nS(i,j) ]

    套用第二类斯特林数的通项公式(S(n,m)=sum_{i=0}^mfrac{(-1)^i*(m-i)^n}{i!(m-i)!})可以得到:

    [sum_{j=0}^n2^j*(j!)sum_{i=0}^nsum_{k=0}^jfrac{(-1)^k*(j-k)^i}{k!(j-k)!} ]

    接下来,我们可以考虑再将所有与(i)无关的项提前,得到下面这个式子:

    [sum_{j=0}^n2^j*(j!)sum_{k=0}^jfrac{(-1)^k}{k!(j-k)!}*sum_{i=0}^n(j-k)^i ]

    我们可以枚举(j),这样前面的(sum_{j=0}^n2^j*(j!))一项就可以轻松搞定。

    那么接下来的问题就是如何对于给定的(j),求出下面这个式子的值:

    [sum_{k=0}^jfrac{(-1)^k}{k!(j-k)!}*sum_{i=0}^n(j-k)^i ]

    仔细观察,可以发现似乎这个式子前半部分的分母和后半部分都有((j-k))这个因式,则容易想到将它们移到一起,即:

    [sum_{k=0}^jfrac{(-1)^k}{k!}*sum_{i=0}^nfrac{(j-k)^i}{(j-k)!} ]

    这样一移的效果是显著的,因为此时这个式子的前半部分只与(k)有关,后半部分只与((j-k))有关。

    则可以设(f)(g)如下:

    [f_x=frac{(-1)^x}{x!} ]

    [g_x=sum_{i=0}^nfrac{x^i}{x!}=frac{sum_{i=0}^nx^i}{x!}=frac{frac{x^{n+1}-1}{x-1}}{x!}=frac{x^{n+1}-1}{(x-1)x!} ]

    带入原式就可以得到:

    [sum_{k=0}^jf_{k}*g_{j-k} ]

    而这个式子实际上就相当于:

    [(f*g)(j) ]

    至此化简完毕,最终的式子就是:

    [sum_{j=0}^n2^j*(j!)*(f*g)(j) ]

    具体实现

    我们可以先(O(n))(O(nlogn))?)预处理出(f)(g)两个数组,然后(NTT)即可。

    这应该还是比较简单的吧。

    代码

    #include<bits/stdc++.h>
    #define Tp template<typename Ty>
    #define Ts template<typename Ty,typename... Ar>
    #define Reg register
    #define RI Reg int
    #define Con const
    #define CI Con int&
    #define I inline
    #define W while
    #define N 100000
    #define X 998244353
    #define Qinv(x) (Qpow(x,X-2))
    #define swap(x,y) (x^=y^=x^=y)
    #define Inc(x,y) ((x+=(y))>=X&&(x-=X))
    #define XSum(x,y) ((x)+(y)>=X?(x)+(y)-X:(x)+(y))
    #define XSub(x,y) ((x)-(y)<0?(x)-(y)+X:(x)-(y))
    using namespace std;
    int n,a[(N<<2)+5],b[(N<<2)+5],Fac[N+5],Inv[N+5];
    I int Qpow(RI x,RI y) {RI res=1;W(y) y&1&&(res=1LL*res*x%X),x=1LL*x*x%X,y>>=1;return res;}
    class NTT//NTT求模意义下的卷积
    {
        private:
            static const int SZ=N,PR=3,IP=(X+1)/3;int P,L,R[(SZ<<2)+5];
            I void Transform(int* s,CI op)
            {
                RI i,j,k,U,S,tx,ty;for(i=0;i^P;++i) i<R[i]&&swap(s[i],s[R[i]]);
                for(i=1;i^P;i<<=1) for(U=Qpow(~op?PR:IP,(X-1)/(i<<1)),j=0;j^P;j+=(i<<1))
                    for(S=1,k=0;k^i;++k,S=1LL*S*U%X) tx=s[j+k],ty=1LL*S*s[i+j+k]%X,s[j+k]=XSum(tx,ty),s[i+j+k]=XSub(tx,ty);
            }
        public:
            I void Solve(CI n,CI m,int* a,int* b)
            {
                RI i,t;P=1,L=0,memset(R,0,sizeof(R));
                W(P<=n+m) P<<=1,++L;for(i=0;i^P;++i) R[i]=(R[i>>1]>>1)|((i&1)<<L-1);
                for(Transform(a,1),Transform(b,1),i=0;i^P;++i) a[i]=1LL*a[i]*b[i]%X;
                for(Transform(a,-1),t=Qinv(P),i=0;i<=n+m;++i) a[i]=1LL*a[i]*t%X;
            }
    }NTT;
    int main()
    {
        RI i,t,ans=0;for(scanf("%d",&n),Fac[0]=1,i=1;i<=n;++i) Fac[i]=1LL*Fac[i-1]*i%X;//预处理阶乘
        for(Inv[n]=Qpow(Fac[n],X-2),i=n-1;~i;--i) Inv[i]=1LL*Inv[i+1]*(i+1)%X;//预处理阶乘的逆元
        for(t=1,i=0;i<=n;++i,t=1LL*t*(X-1)%X) a[i]=1LL*t*Inv[i]%X;//预处理第一个数组
        for(b[0]=1,b[1]=n+1,i=2;i<=n;++i) b[i]=1LL*(Qpow(i,n+1)-1)*Qinv(i-1)%X*Inv[i]%X;//预处理第二个数组
        for(NTT.Solve(n,n,a,b),t=1,i=0;i<=n;++i,(t<<=1)>=X&&(t-=X)) Inc(ans,1LL*a[i]*t%X*Fac[i]%X);//做一遍NTT,然后统计答案
        return printf("%d",ans),0;//输出答案
    }
    
  • 相关阅读:
    Asp.net获取客户端的IP地址排除::1
    EF 筛选列包含NULL会报错
    layUI关于table编辑列支持方向键功能
    .NET CORE 发布IIS问题收集
    VS2019最新源代码管理工具+附下载地址
    关于Mysql可视化工具Navicat Premium12激活使用【亲测】
    经典SQL 语句
    事务的四种隔离级别 [转载]
    HTML 特殊符号编码对照表
    github本地文件Push到仓库
  • 原文地址:https://www.cnblogs.com/chenxiaoran666/p/BZOJ4555.html
Copyright © 2011-2022 走看看