半在线卷积小记

广告 : 多项式计数杂谈

• P4721 【模板】分治 FFT

题意 : 给出多项式F(x),求G使得:

• G[0]=1
• G[n]=\sum\limits_{i=1}^nG[n-i]F[i]

常识数据范围。

这个题可以直接求逆做，由于和主题无关我们就不介绍了。

此题中，某个位置的系数依赖于前面的计算，而又会向后面贡献。

这和一些数据结构优化 DP 题有点类似，实际上,这就是强制中序遍历转移的 CDQ 分治。

考虑每次计算跨越隔板的贡献。

式子是卷积的形式，且此时没有依赖性，可以直接用 FFT 计算，复杂度是 O(n\log^2n) 的。

来模拟一下，令F=<0,1,3,2>

G=< 1, 0, 0, 0 >

G=<[1 |0] 0, 0 >

使用<1>卷上F的一部分 : <0,1> 得到<*,1>,将其加到对应的位置。

* 表示我们不在意那个位置的具体值。

G=<[1 |1] 0, 0 >

G=<[1, 1 |0, 0]>

使用 <1,1> 卷上 F 的一部分 : <0,1,3,2> 得到 <*,*,4,5> ,将其加到对应的位置。

G=<[1, 1 |4, 5]>

G=< 1, 1 [4 |5]>

使用 <4> 卷上 F 的一部分 : <0,1> 得到 <*,4> ,将其加到对应的位置。

注意，并非和后半部的 <3,2> 进行卷积。

G=< 1, 1, 4, 9 >

这就是答案了。

可以先填充成 2 的幂,这样会比较好写。

由于我们不关心前半部分的结果，我们可以把循环卷积的次数定为 O(len) 而非 O(1.5len) ,这样就算溢出也只会损坏我们不关心的前半部分。这能够较为可观地优化常数，耗时为原来的60~75%。

又能发现，我们只会使用 F 的某个前缀的 DFT 结果，可以预处理出来。耗时进一步减少到原来的70~80%。

#include<algorithm>
#include<cstring>
#include<cstdio>
#define ll long long
#define ull unsigned long long
#define clr(f,n) memset(f,0,sizeof(int)*(n))
#define cpy(f,g,n) memcpy(f,g,sizeof(int)*(n))
using namespace std;
const int _G=3,mod=998244353,Maxn=135000;
inline int read(){
  int X=0;char ch=0;
  while(ch<48||ch>57)ch=getchar();
  while(ch>=48&&ch<=57)X=X*10+(ch^48),ch=getchar();
  return X;
}
ll powM(ll a,int t=mod-2){
  ll ans=1;
  while(t){
    if(t&1)ans=ans*a%mod;
    a=a*a%mod;t>>=1;
  }return ans;
}
const int invG=powM(_G);
int tr[Maxn],tf;
void tpre(int n){
  if (tf==n)return ;tf=n;
  for(int i=0;i<n;i++)
    tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
}
void NTT(int *g,bool op,int n)
{
  tpre(n);
  static ull f[Maxn],w[Maxn]={1};
  for (int i=0;i<n;i++)f[i]=g[tr[i]];
  for(int l=1;l<n;l<<=1){
    ull tG=powM(op?_G:invG,(mod-1)/(l+l));
    for (int i=1;i<l;i++)w[i]=w[i-1]*tG%mod;
    for(int k=0;k<n;k+=l+l)
      for(int p=0;p<l;p++){
        int tt=w[p]*f[k|l|p]%mod;
        f[k|l|p]=f[k|p]+mod-tt;
        f[k|p]+=tt;
      }      
    if (l==(1<<17))
      for (int i=0;i<n;i++)f[i]%=mod;
  }if (!op){
    ull invn=powM(n);
    for(int i=0;i<n;++i)
      g[i]=f[i]%mod*invn%mod;
  }else for(int i=0;i<n;++i)g[i]=f[i]%mod;
}
void px(int *f,int *g,int n)
{for(int i=0;i<n;++i)f[i]=1ll*f[i]*g[i]%mod;}
int G[Maxn],s1[Maxn],s2[Maxn<<1];
void cdq(int l,int r)
{
  if (l+1==r)return ;
  int mid=(l+r)>>1,n=r-l;
  cdq(l,mid);
  cpy(s1,G+l,n/2);clr(s1+(n/2),n/2);NTT(s1,1,n);
  px(s1,s2+n,n);NTT(s1,0,n);
  for (int i=n/2;i<n;i++)
    G[l+i]=(G[l+i]+s1[i])%mod;
  cdq(mid,r);
}
int n,m,F[Maxn];
int main()
{
  m=read();G[0]=1;
  for (int i=1;i<m;i++)F[i]=read();
  for (n=1;(n>>1)<m;n<<=1){
    cpy(s1,F,n);NTT(s1,1,n);
    cpy(s2+n,s1,n);
  }cdq(0,n>>=1);
  for (int i=0;i<m;i++)
    printf("%d ",G[i]);
  return 0;
}

• 多项式 \exp

在有标号计数的推导中我们经常要用到多项式 \exp ，然而在考场上写 O(n\log n) 的牛顿迭代 \exp 性价比并不高。

这时就要请我们的半在线卷积 Exp 上场了,复杂度是 O(n\log^2 n) 的，但由于常数小并不会比牛顿迭代慢。

( 我才不会告诉你们这玩意比我使劲卡常的牛顿迭代Exp快好多 )

设 B(x)=e^{F(x)} ，两边求导可得 B'(x)=B(x)F'(x)

提取第 n 次项系数可得 : (n+1)B[n+1]=\sum\limits_{i=0}^{n}(i+1)F[i+1]B[n-i]

B[n]=\frac{1}{n}\sum\limits_{i=1}^{n}iF[i]B[n-i]

这就回到了经典问题,常数\frac{1}{n}可以在分治边界乘上。

#include<algorithm>
#include<cstring>
#include<cstdio>
#define ll long long
#define ull unsigned long long
#define clr(f,n) memset(f,0,sizeof(int)*(n))
#define cpy(f,g,n) memcpy(f,g,sizeof(int)*(n))
using namespace std;
const int _G=3,mod=998244353,Maxn=1050000;
inline int read(){
  int X=0;char ch=0;
  while(ch<48||ch>57)ch=getchar();
  while(ch>=48&&ch<=57)X=X*10+(ch^48),ch=getchar();
  return X;
}
ll powM(ll a,int t=mod-2){
  ll ans=1;
  while(t){
    if(t&1)ans=ans*a%mod;
    a=a*a%mod;t>>=1;
  }return ans;
}
const int invG=powM(_G);
int tr[Maxn<<1],tf;
void tpre(int n){
  if (tf==n)return ;tf=n;
  for(int i=0;i<n;i++)
    tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
}
void NTT(int *g,bool op,int n)
{
  tpre(n);
  static ull f[Maxn<<1],w[Maxn<<1]={1};
  for (int i=0;i<n;i++)f[i]=g[tr[i]];
  for(int l=1;l<n;l<<=1){
    ull tG=powM(op?_G:invG,(mod-1)/(l+l));
    for (int i=1;i<l;i++)w[i]=w[i-1]*tG%mod;
    for(int k=0;k<n;k+=l+l)
      for(int p=0;p<l;p++){
        int tt=w[p]*f[k|l|p]%mod;
        f[k|l|p]=f[k|p]+mod-tt;
        f[k|p]+=tt;
      }      
    if (l==(1<<17))
      for (int i=0;i<n;i++)f[i]%=mod;
  }if (!op){
    ull invn=powM(n);
    for(int i=0;i<n;++i)
      g[i]=f[i]%mod*invn%mod;
  }else for(int i=0;i<n;++i)g[i]=f[i]%mod;
}
void px(int *f,int *g,int n)
{for(int i=0;i<n;++i)f[i]=1ll*f[i]*g[i]%mod;}
int F[Maxn],G[Maxn],s1[Maxn],s2[Maxn];
void cdq(int l,int r)
{
  if (l+1==r){
    if (l>0)G[l]=powM(l)*G[l]%mod;
    return ;
  }int mid=(l+r)>>1,n=r-l;
  cdq(l,mid);
  cpy(s1,G+l,n/2);clr(s1+(n/2),n/2);NTT(s1,1,n);
  cpy(s2,F,n);NTT(s2,1,n);
  px(s1,s2,n);NTT(s1,0,n);
  for (int i=n/2;i<n;i++)
    G[l+i]=(G[l+i]+s1[i])%mod;
  cdq(mid,r);
}
int n,m;
int main()
{
  m=read();G[0]=1;
  for (int i=0;i<m;i++)F[i]=1ll*read()*i%mod;
  for (n=1;n<m;n<<=1);
  cdq(0,n);
  for (int i=0;i<m;i++)
    printf("%d ",G[i]);
  return 0;
}

• 更加真实的半在线卷积

  已知 F[0]=G[0]=0,\ F[1]=1

  令 G[n]=F[0...n] 的异或和。

  F[n]=\sum\limits_{i=1}^{n-1}G[i]F[n-i]

前面的问题都是事先给定 G ,此时 G[n] 要在 F[0...n] 都得到之后才能算出。

前面我们的策略是 : F[l,mid]*G[0,r-l)\rightarrow F(mid,r)

当 l=0 时是这样的 : F[0,mid]*G[0,r-l)\rightarrow F(mid,r)

正确的策略是 : - $l=0$ 时 $F[0,mid]*G[0,mid]\rightarrow F(mid,r)

这样,暂时只有 F[mid+1] 是正确的。余者没有考虑 F[0,mid]*G(mid,r-l) 的贡献。

• $\left.\begin{matrix}G[l,mid]*F[0,r-l) \\ F[l,mid]*G[0,r-l) \end{matrix}\right\}\rightarrow F(mid,r)

  容易发现，此时左侧涉及到的所有系数都已经计算完毕。

  下半部是经典的，上半部就是在补充 l=0 时漏掉的贡献。

• 应用 : [数学记录]P5900 无标号无根树计数

• 更加真实的半在线卷积 II

  已知 F[0]=1 ,给出 G。

  F[n]=\sum\limits_{i=1}^{n}G[i]F^2[n-i]

此时，当 F 多算得 F[n] 一位时，并不能马上快速更新 F^2[n]。

使用半在线卷积计算 F^2 ，然后并行计算即可。

• 应用 : [数学记录]Uoj#50. 【UR #3】链式反应