command_block 的博客

command_block 的博客

半在线卷积小记

posted on 2020-05-05 22:29:14 | under 算法 |

广告 : 多项式计数杂谈

题意 : 给出多项式$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)$

$G(mid,r-l)$ 还没有被正确计算出来。直接这样做是错误的。

正确的策略是 :

  • $l=0$ 时

    $F[0,mid]*G[0,mid]\rightarrow F(mid,r)$

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

  • $l>0$ 时

    $\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$ ,然后并行计算即可。