BZOJ 4555 [Tjoi2016&Heoi2016]求和

最新推荐文章于 2019-07-04 22:58:56 发布

wang3312362136

最新推荐文章于 2019-07-04 22:58:56 发布

阅读量245

点赞数

CC 4.0 BY-SA版权

分类专栏： FFT

本文链接：https://blog.youkuaiyun.com/wang3312362136/article/details/81699757

FFT 专栏收录该内容

8 篇文章

订阅专栏

本文详细解析了一道算法竞赛题目，该题目要求计算一个包含斯特林数的复杂序列的和。通过数学推导简化了问题，并利用快速傅立叶变换（NTT）优化了计算过程。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

题目链接

https://www.lydsy.com/JudgeOnline/problem.php?id=4555

题解

f (n) = \sum i = 0 n \sum j = 0 i S (i, j) \times 2 j \times j! = \sum i = 0 n \sum j = 0 n S (i, j) \times 2 j \times j! = \sum i = 0 n \sum j = 0 n 2 j \sum k = 0 j (- 1) k (j k) (j - k) i = \sum i = 0 n \sum j = 0 n 2 j \times j! \sum k = 0 j ( - 1 ) k k ! \times ( j - k ) i ( j - k ) ! = \sum j = 0 n 2 j \times j! \sum k = 0 j ( - 1 ) k k ! \times \sum n i = 0 ( j - k ) i ( j - k ) ! (1) (2) (1) (3) (2)

$\begin{align} f(n)&=\sum_{i=0}^n \sum_{j=0}^i S(i,j)\times 2^j\times j!\\ &=\sum_{i=0}^n \sum_{j=0}^n S(i,j)\times 2^j\times j!\\ &=\sum_{i=0}^n \sum_{j=0}^n 2^j\sum_{k=0}^j (-1)^k\binom{j}{k}(j-k)^i\tag{1}\\ &=\sum_{i=0}^n \sum_{j=0}^n 2^j\times j!\sum_{k=0}^j\frac{(-1)^k}{k!}\times \frac{(j-k)^i}{(j-k)!}\\ &=\sum_{j=0}^n 2^j\times j!\sum_{k=0}^j\frac{(-1)^k}{k!}\times \frac{\sum_{i=0}^n (j-k)^i}{(j-k)!}\tag{2} \end{align}$

其中 $(1)$ 步骤使用了第二类斯特林数的展开式 $S(i,j)=\frac{1}{j!}\sum_{k=0}^j (-1)^k \binom{j}{k}(j-k)^i$ ， $(2)$ 步骤是一个卷积形式，模数比较特殊可以用NTT优化， $\sum_{i=0}^n (j-k)^i$ 很明显是一个等比数列求和。

代码

#include <cstdio>
#include <algorithm>

int read()
{
  int x=0,f=1;
  char ch=getchar();
  while((ch<'0')||(ch>'9'))
    {
      if(ch=='-')
        {
          f=-f;
        }
      ch=getchar();
    }
  while((ch>='0')&&(ch<='9'))
    {
      x=x*10+ch-'0';
      ch=getchar();
    }
  return x*f;
}

const int maxn=200000;
const int maxm=340000;
const int mod=998244353;
const int G=3;

int quickpow(int a,int b,int m)
{
  int res=1;
  while(b)
    {
      if(b&1)
        {
          res=1ll*res*a%m;
        }
      a=1ll*a*a%m;
      b>>=1;
    }
  return res;
}

int add(int x,int y,int m)
{
  int res=x+y;
  if(res>=m)
    {
      res-=m;
    }
  return res;
}

int minus(int x,int y,int m)
{
  int res=x-y;
  if(res<0)
    {
      res+=m;
    }
  return res;
}

int rev[maxm+10],a[maxm+10],b[maxm+10],ans[maxm+10];

int getrev(int n)
{
  int m=1,len=0;
  while(m<=n)
    {
      m<<=1;
      ++len;
    }
  for(int i=1; i<m; ++i)
    {
      rev[i]=(rev[i>>1]>>1)+((i&1)<<(len-1));
    }
  return m;
}

int fft(int *s,int len)
{
  for(int i=0; i<len; ++i)
    {
      if(rev[i]<i)
        {
          std::swap(s[rev[i]],s[i]);
        }
    }
  for(int i=2; i<=len; i<<=1)
    {
      int gn=quickpow(G,(mod-1)/i,mod);
      for(int j=0; j<len; j+=i)
        {
          int g=1;
          for(int k=0; k<(i>>1); ++k)
            {
              int x=s[j+k],y=1ll*g*s[j+k+(i>>1)]%mod;
              s[j+k]=add(x,y,mod);
              s[j+k+(i>>1)]=minus(x,y,mod);
              g=1ll*g*gn%mod;
            }
        }
    }
  return 0;
}

int main()
{
  int n=read();
  a[0]=1;
  int v=1;
  for(int i=1; i<=n; ++i)
    {
      v=1ll*minus(0,v,mod)*quickpow(i,mod-2,mod)%mod;
      a[i]=v;
    }
  b[0]=1;
  v=1;
  for(int i=1; i<=n; ++i)
    {
      v=1ll*v*quickpow(i,mod-2,mod)%mod;
      b[i]=1ll*minus(1,quickpow(i,n+1,mod),mod)*quickpow(minus(1,i,mod),mod-2,mod)%mod*v%mod;
    }
  b[1]=n+1;
  int m=getrev(n<<1);
  fft(a,m);
  fft(b,m);
  for(int i=0; i<m; ++i)
    {
      ans[i]=1ll*a[i]*b[i]%mod;
    }
  fft(ans,m);
  std::reverse(ans,ans+m+1);
  v=quickpow(m,mod-2,mod);
  for(int i=0; i<m; ++i)
    {
      ans[i]=1ll*ans[i]*v%mod;
    }
  v=1;
  int u=1;
  for(int i=0; i<=n; ++i)
    {
      v=add(v,1ll*quickpow(2,i,mod)*u%mod*ans[i]%mod,mod);
      u=1ll*u*(i+1)%mod;
    }
  printf("%d\n",v);
  return 0;
}