zoukankan      html  css  js  c++  java
  • hdu 5279 YJC plays Minecraft——生成函数

    题目:http://acm.hdu.edu.cn/showproblem.php?pid=5279

    令 n 个点的树的 EGF 是 g(x) ,则 ( g(x) = sumlimits_{i=0}^{infty} frac{i^{i-2}}{i!} x^i )

    令 n 个点的森林的 EGF 是 f(x) ,则 ( f(x) = sumlimits_{i=0}^{infty} frac{g(x)^i}{i!} = e^{g(x)} )

    这道题里,每个团调用 f(x) 的对应项系数,n 条边就给答案乘上 2n

    但要减去所有团、团之间的边构成一个环的情况。所以考虑 n 个点、1号点与 n 号点在同一棵树里的森林个数的 EGF h(x) 。

    有这样的递推式: ( h[n]=sumlimits_{i=1}^{n-1} g[n-i]*f[i]*inom{n-2}{i-1} ) ,就是考虑与 1 号点在同一棵树里有 (n-i) 个点。

    即 ( frac{h[n]}{(n-2)!} = sumlimits_{i=1}^{n-1} frac{g[n-i]}{(n-i-1)!} * frac{f[i]}{(i-1)!} )

    到这里,自己本来想的是后面的就是 ( x*f'(x) ) 和 ( x*g'(x) ) ,然后卷积。不过看看题解,发现前面是 h(x) 的二阶导,后面是 f(x) 和 g(x) 的一阶导。

    那么 ( h(x) = int int f'(x)*g'(x) dx dx )

    注意数组长度是 219 而不是 218 ,因为倍增的 %t 的 t 可能是 218 ,此时的 len 应该是 219

    注意预处理逆元要到 219 ,而不是和 jc[ ] ,jcn[ ] 一样只是 105

    注意积分的时候要倒序遍历,以防 a[ ] 和 b[ ] 是同一个数组。

    注意 get_exp( ) 的时候,本层的 len 要在 get_inv( ) 之后再赋值。并注意 ( 1-ln(f_0(x))+a(x) ) 的那个 1 是常数项,不是每项都有一个+1。

    注意 ( i^{i-2} ) 里的 i 不能是 1 。

    注意最后算 h(x) 的时候做积分,随便找了两个数组存积分结果,在乘起来之前要先把无关的项清空!!!

    注意最后乘上 ( i! ) 才是真正的系数。

    #include<cstdio>
    #include<cstring>
    #include<algorithm>
    #define ll long long
    using namespace std;
    int rdn()
    {
      int ret=0;bool fx=1;char ch=getchar();
      while(ch>'9'||ch<'0'){if(ch=='-')fx=0;ch=getchar();}
      while(ch>='0'&&ch<='9')ret=ret*10+ch-'0',ch=getchar();
      return fx?ret:-ret;
    }
    const int N=(1<<19)+5,mod=998244353;//1<<19 not 18 for t=2^18
    int upt(int x){while(x>=mod)x-=mod;while(x<0)x+=mod;return x;}
    int pw(int x,int k)
    {int ret=1;while(k){if(k&1)ret=(ll)ret*x%mod;x=(ll)x*x%mod;k>>=1;}return ret;}
    
    int n,ct[N],f[N],g[N],h[N],jc[N],jcn[N],bin[N];
    int la[N],ia[N],tp[N],A[N],len,r[N],inv[N],wn[N],wn2[N];
    void ntt_pre(int len)
    {
      for(int R=2;R<=len;R++)
        {
          wn[R]=pw(3,(mod-1)/R);
          wn2[R]=pw(3,(mod-1)-(mod-1)/R);
        }
      inv[1]=1;
      for(int i=2;i<=len;i++)/////here not in init()
        inv[i]=(ll)upt(-mod/i)*inv[mod%i]%mod;
    }
    void ntt_len()
    {
      for(int i=0,j=len>>1;i<len;i++)
        r[i]=(r[i>>1]>>1)+((i&1)?j:0);
    }
    void ntt(int *a,bool fx)
    {
      for(int i=0;i<len;i++)
        if(i<r[i])swap(a[i],a[r[i]]);
      for(int R=2;R<=len;R<<=1)
        {
          int Wn=fx?wn2[R]:wn[R];
          for(int i=0,m=R>>1;i<len;i+=R)
        for(int j=0,w=1;j<m;j++,w=(ll)w*Wn%mod)
          {
            int x=a[i+j],y=(ll)w*a[i+m+j]%mod;
            a[i+j]=upt(x+y); a[i+m+j]=upt(x-y);
          }
        }
      if(!fx)return; int iv=inv[len];
      for(int i=0;i<len;i++)a[i]=(ll)a[i]*iv%mod;
    }
    void get_dao(int n,int *a,int *b)
    {
      for(int i=0;i<n;i++)b[i]=(ll)(i+1)*a[i+1]%mod;
      b[n]=0;
    }
    void get_jf(int n,int *a,int *b)
    {
      for(int i=n;i;i--)b[i]=(ll)inv[i]*a[i-1]%mod;//
      b[0]=0;
    }
    void get_inv(int n,int *a,int *b)
    {
      b[0]=pw(a[0],mod-2);
      for(int t=2,yt=1,i,j;yt<n;yt=t,t=len)
        {
          len=t<<1; ntt_len();
          for(i=0;i<t;i++)tp[i]=a[i];
          for(;i<len;i++)tp[i]=0; for(i=yt;i<len;i++)b[i]=0;
          //b[i]=0 here for other get_inv pollute b[]//<len or t?
          ntt(tp,0); ntt(b,0);
          for(i=0;i<len;i++)b[i]=upt((ll)b[i]*(2-(ll)b[i]*tp[i]%mod)%mod);
          ntt(b,1);
        }
      for(int i=n;i<len;i++)b[i]=0;
    }
    void get_ln(int n,int *a,int *b)
    {
      for(int i=0;i<n;i++)ia[i]=0;
      get_dao(n,a,b); get_inv(n,a,ia);
      len=n<<1; ntt_len();//len=n<<1 is ok
      ntt(b,0); ntt(ia,0);
      for(int i=0;i<len;i++)b[i]=(ll)b[i]*ia[i]%mod;
      ntt(b,1); for(int i=n;i<len;i++)b[i]=0;
      get_jf(n-1,b,b);
    }
    void get_exp(int n,int *a,int *b)
    {
      b[0]=1;
      for(int t=2,yt=1,i,j;yt<n;yt=t,t=len)
        {
          for(i=yt;i<t;i++)b[i]=0; get_ln(t,b,la);//b[i]=0 before
          for(i=0;i<t;i++)la[i]=upt(-la[i]+a[i]);
          la[0]=upt(la[0]+1);/////not la[i]=1-la[i]+a[i]!!!
          len=t<<1; ntt_len();//////after get_ln!!!!!
          ntt(la,0); ntt(b,0);
          for(i=0;i<len;i++)b[i]=(ll)la[i]*b[i]%mod;
          ntt(b,1);
        }
      for(int i=n;i<len;i++)b[i]=0;
    }
    void init()
    {
      n=1e5; for(len=1;len<=n;len<<=1);//len=mx_t
      ntt_pre(len<<1);
      jc[0]=1;for(int i=1;i<=n;i++)jc[i]=(ll)jc[i-1]*i%mod;
      jcn[n]=pw(jc[n],mod-2);
      for(int i=n-1;i>=0;i--)jcn[i]=(ll)jcn[i+1]*(i+1)%mod;
      bin[0]=1;for(int i=1;i<=n;i++)bin[i]=upt(bin[i-1]<<1);
    
      g[1]=1;
      for(int i=2;i<=n;i++)//i=2 not i=1!!
        g[i]=(ll)pw(i,i-2)*jcn[i]%mod;
      get_exp(n+1,g,f);
      get_dao(n,g,ia); get_dao(n,f,la);
      for(len=1;len<n<<1;len<<=1); ntt_len();
      for(int i=n;i<len;i++)ia[i]=la[i]=0;//////!!!!!
      ntt(ia,0); ntt(la,0);
      for(int i=0;i<len;i++)ia[i]=(ll)ia[i]*la[i]%mod;
      ntt(ia,1);
      get_jf(n,ia,h); get_jf(n,h,h);
      for(int i=1;i<=n;i++)///
        {
          f[i]=(ll)f[i]*jc[i]%mod;
          g[i]=(ll)g[i]*jc[i]%mod;
          h[i]=(ll)h[i]*jc[i]%mod;
        }
    }
    int main()
    {
      int T=rdn(); init();
      while(T--)
        {
          n=rdn();
          for(int i=1;i<=n;i++)ct[i]=rdn();
          int m1=bin[n],m2=1;
          for(int i=1;i<=n;i++) m1=(ll)m1*f[ct[i]]%mod;
          for(int i=1;i<=n;i++) m2=(ll)m2*upt(f[ct[i]]-h[ct[i]])%mod;
          printf("%d
    ",upt(m1-m2));
        }
      return 0;
    }
  • 相关阅读:
    MATLAB 中sparse函数使用及full函数用法简单介绍(转)
    稀疏矩阵加减,乘除, 逆 (转)
    拟合方法求直线方程系数
    matlab filtfilt 函数
    Typora 精美而强大的Markdown编辑器 转
    MATLAB生成exe脱离matlab运行可执行程序
    matlab 生成.exe文件 转
    C#排序 转
    C# 进制转换(二进制、十六进制、十进制互转) 转载 https://www.cnblogs.com/icebutterfly/p/8884023.html
    一维高斯滤波 转
  • 原文地址:https://www.cnblogs.com/Narh/p/10829343.html
Copyright © 2011-2022 走看看