下标模意义下的多项式乘法及其应用

已知两个数组\(a,b\)，求一个数组\(c\)，满足\(c_i=\sum_{j+k\equiv i(MOD\ K)}a_jb_k(MOD\ M)\)。这里我们把这个东西称为下标模意义下的多项式乘法。那么这个东西怎么做呢？

先说结论：如果MOD M意义下存在K次单位根，那么把平时用NTT做多项式乘法时的原根G统一换成这个K次单位根就可以了(不能跑\(O(nlogn)\)的NTT，因为把数组长度扩充到2的倍数会导致出错)。

复习一下单位根的定义

这是为什么呢？

举个例子，当K=17，M=998244353的时候。此时MOD M意义下是存在K次单位根的，它的值是503044，令其为\(w\)。假设现在有一个多项式\(a\)，且项数不止17个，令它对\(x^{17}\)取模之后的多项式为\(b\)，\(b\)只有\(17\)项，且\(b_i=a_i+a_{i+17}+a_{i+34}\cdots\)。以503044作为单位根的值对\(a,b\)分别做一次DFT(系数转点值)，得到\(c,d\)。由于\(w^{17}=1\)，所以\(c,d\)都有长度为17的循环节，并且容易发现\(c=d\)。这时候截取\(c\)的前17项做IDFT，就可以得到\(b\)，也就是下标取模之后的多项式。这说明只要对一个多项式先DFT，取前17项再IDFT，就可以完成一次下标的取模。

那么两个多项式相乘也是一样的，最后把点值相乘的结果IDFT回去的时候只取前17项就行了。

但是注意到\(O(nlogn)\)的NTT不能这么搞，因为NTT会把数组长度先扩充到2的倍数，但是这个问题里数组长度不是17就不对了。

放一张毛啸2016年候选队论文的截图，仅供参考：

例题

n 个点 m 条带权无向边，权值是 0 到 16。问多少种方案选择一些边（不选重边），使得图联通，且边权和模 17 正好为 x。对每个 0≤x≤16 都求答案。\(2 \le n \le 17,1\le m\le 10^5\)。

解法

前置知识：

1. 二维FFT(高维FFT)。有两个矩阵\(a,b\)，现在要求一个新的矩阵\(c\)满足\(c_{i,j}=\sum_{u+v=i,x+y=j}a_{u,x}b_{v,y}\)。做法
2. 子集exp和子集ln。一个子集幂级数指的是一个数组\(f_s\)，下标\(s\)是一个"子集"，可以看成是\([0,2^n-1]\)中的一个整数，用二进制的方式表示n个元素中的哪些在子集中。子集exp和ln则类似EGF的exp和ln。子集exp的组合意义：如果\(g=e^f\)，且\(f_s\)表示将子集\(s\)中所有的元素进行操作A的方案是，则\(g_s\)表示把\(s\)进一步分成若干子集，每个子集进行A操作的总方案数。子集ln则是子集exp的逆运算，用来把g变成f的。关于子集exp和ln的实现方法和证明具体见这里。

首先令\(f_{i,s}\)表示子集\(s\)内部的边全部乱选，且不能选重边，和MOD 17为i的方案数。\(s\)内部的边指的是两段都在\(s\)中的边。同样定义\(g_{i,s}\)表示只选\(s\)中的边，图连通，没有重边且和MOD 17为i的方案数。观察得到转移式：令\(l\)表示\(s\)中的最低的1，则\(g_{i,s}=f_{i,s}-\sum_{t\subsetneq s,l\in t,j+k=i(MOD\ 17)}g_{j,t}\cdot f_{k,s \ xor \ t}\)，也就是枚举l所在的连通块。直接写这玩意儿能做到\(O(3^n\cdot17^2)\)。

如果所有边的权值都为0，那就可以把\(f,g\)的第一维去掉，只留下子集那一维。根据子集exp的组合意义，很容易发现\(f=e^g\)，\(g=ln(f)\)。根据上面链接里的子集ln求法，直接\(O(2^n\cdot n^2)\)求\(f\)的ln即可。\(f\)本身可以通过简单的dp求出，复杂度是\(O(2^n\cdot n\cdot 17^2)\)。在所有边权都为0的时候则可以直接\(O(2^n\cdot n)\)求出。

边权值不全为0的情况则需要加上第一维。注意到第一维是一个上面说的"下标模意义下的多项式乘法"。可以根据二维FFT的方法依葫芦画瓢，先把第一维以503044位单位根进行DFT，然后把第二维子集ln，最后再把第一维IDFT。这样套娃的正确性我不会证明，感性理解。

\(n\le 17\)，所以把n用17表示了。总复杂度\(O(17^3\cdot 2^n)\)。

#include <bits/stdc++.h>

#define rep(i,n) for(int i=0;i<n;++i)
#define repn(i,n) for(int i=1;i<=n;++i)
#define LL long long
#define pii pair <int,int>
#define fi first
#define se second
#define mpr make_pair
#define pb push_back

void fileio()
{
  #ifdef LGS
  freopen("in.txt","r",stdin);
  freopen("out.txt","w",stdout);
  #endif
}
void termin()
{
  #ifdef LGS
  std::cout<<"\n\nPROGRAM TERMINATED";
  #endif
  exit(0);
}

using namespace std;

const LL MOD=998244353;
const LL w=503044;

LL qpow(LL x,LL a)
{
	LL res=x,ret=1;
	while(a>0)
	{
		if(a&1) (ret*=res)%=MOD;
		a>>=1;
		(res*=res)%=MOD;
	}
	return ret;
}

LL n,m,es[20][20][20],f[20][140000],g[20][140000],dp[20][20],modv[20][20],inv17,inv[20];
//f：乱选(但不选重边)；g：连通

vector <LL> DFT(vector <LL> v)
{
  LL cur=1;
  vector <LL> ret;
  rep(i,v.size())
  {
    LL curr=1,res=0;
    rep(j,v.size())
    {
      (res+=v[j]*curr)%=MOD;
      (curr*=cur)%=MOD;
    }
    (cur*=w)%=MOD;
    ret.pb(res);
  }
  return ret;
}
vector <LL> IDFT(vector <LL> v)
{
  LL cur=1,ww=qpow(w,MOD-2);
  vector <LL> ret;
  rep(i,v.size())
  {
    LL curr=1,res=0;
    rep(j,v.size())
    {
      (res+=v[j]*curr)%=MOD;
      (curr*=cur)%=MOD;
    }
    (cur*=ww)%=MOD;
    ret.pb(res);
  }
  rep(i,ret.size()) (ret[i]*=inv17)%=MOD;
  return ret;
}

void add(LL &x,LL y){x+=y;if(x>=MOD) x-=MOD;}

LL cur[140000],h[20][140000],fr[20],to[20];
void subsetLn()
{
  rep(i,18) rep(j,1<<n) h[i][j]=0;
  rep(i,1<<n) h[__builtin_popcount(i)][i]=cur[i];
  //第三维FMT
  rep(i,n+1)
  {
    for(int j=n-1;j>=0;--j) rep(k,1<<n) if((k&(1<<j))==0)
      add(h[i][k|(1<<j)],h[i][k]);
  }
  //第二维ln
  rep(i,1<<n)
  {
    rep(j,n+1) fr[j]=h[j][i],to[j]=0;
    repn(j,n)
    {
      to[j]=fr[j];
      LL sum=0;
      rep(k,j) (sum+=to[k]*fr[j-k]%MOD*k)%=MOD;
      (sum*=inv[j])%=MOD;
      add(to[j],MOD-sum); 
    }
    rep(j,n+1) h[j][i]=to[j];
  }
  //第三维IFMT
  rep(i,n+1)
  {
    rep(j,n) rep(k,1<<n) if((k&(1<<j))==0)
      add(h[i][k|(1<<j)],MOD-h[i][k]);
  }

  rep(i,1<<n) cur[i]=h[__builtin_popcount(i)][i];
}

int main()
{
	fileio();
  //freopen("roads.in","r",stdin);
  //freopen("roads.out","w",stdout);

  inv17=qpow(17,MOD-2);
  rep(i,19) inv[i]=qpow(i,MOD-2);
  rep(i,18) rep(j,18) modv[i][j]=(i+j)%17;
  cin>>n>>m;
  LL x,y,z;
  rep(i,m)
  {
    scanf("%lld%lld%lld",&x,&y,&z);--x;--y;
    if(x>y) swap(x,y);
    ++es[x][y][z];
  }
  rep(i,n) for(int j=i+1;j<n;++j) ++es[i][j][0];
  rep(msk,1<<n) if(msk>0)
  {
    int lowbit=msk&(-msk),lft=msk^lowbit,id=__builtin_ctz(msk);
    if(lft==0) f[0][msk]=1;
    else
    {
      vector <int> nds;rep(j,n) if(lft&(1<<j)) nds.pb(j);
      rep(i,18) rep(j,18) dp[i][j]=0;
      rep(i,17) dp[0][i]=f[i][lft];
      rep(i,nds.size())
      {
        int ii=nds[i];
        rep(j,17) rep(k,17) (dp[i+1][modv[j][k]]+=dp[i][j]*es[id][ii][k])%=MOD;
      }
      rep(i,17) f[i][msk]=dp[nds.size()][i];
    }
  }

  //第一维DFT
  rep(i,1<<n)
  {
    vector <LL> vv;
    rep(j,17) vv.pb(f[j][i]);
    vv=DFT(vv);
    rep(j,17) g[j][i]=vv[j];
  }
  //子集ln
  rep(i,17)
  {
    rep(j,1<<n) cur[j]=g[i][j];
    subsetLn();
    rep(j,1<<n) g[i][j]=cur[j];
  }
  //第二维IDFT
  rep(i,1<<n)
  {
    vector <LL> vv;
    rep(j,17) vv.pb(g[j][i]);
    vv=IDFT(vv);
    rep(j,17) g[j][i]=vv[j];
  }

  rep(i,17) printf("%lld ",g[i][(1<<n)-1]);
  puts("");

	termin();
}