对于非负整数 n,m，设 B(n,m) 为 n 个点 m 条边的二分图个数。给定非负整数 N 和 \textbf{\textit{A}}\in \mathbb F_{2^{64}}，对于每个 n=0,1,\cdots,N，试求 \sum\limits_{m}B(n,m)\textbf{\textit{A}}^m。
## 题解
### 思路
先考虑在模 $998\,244\,353$ 意义下怎么做。这是经典的，设 $F(z)$ 为答案的 EGF，而 $G(z)$ 是每个二分图染色的权值之和对应的 EGF（即，每个二分图会被算 $2$ 的连通块个数次方次），则显然有 $F(z)^2=G(z)$。（可以构造双射）
但是这个做法不能照搬到 $\mathbb F_{2^{64}}$ 上。一个浅显的原因是**除以** $2$ 做不了所以没有 EGF。更深刻的原因是即使你定义两个**序列**的 EGF 卷积为
$$
\textbf{\textit{c}}_k=\sum_{i+j=k}\dbinom{k}{i}\textbf{\textit{a}}_i\textbf{\textit{b}}_j
$$
也会发现任何**序列** $\{\textbf{\textit{a}}_n\}$ 与自身的卷积为 $\textbf{1},\textbf{0},\textbf{0},\cdots$，或者说单位元。因此你无法开根号。
那怎么办？肯定还是要用染色来考虑，因为二分图只有这一个容易考虑的了。但是二分图的染色方案数一定是偶数啊？考虑钦定 $1$ 号点的颜色为黑色，这时**贡献**非 $\textbf{0}$ 的二分图一定是连通二分图。换言之，我们可以容易求出连通二分图对应的**答案**：
$$
\textbf{\textit{f}}_n=\sum_{i=1}^n\dbinom{n-1}{i-1}\textbf{(\textit{A}+1)}^{i(n-i)}
$$
其中有经典的 $\textbf{\textit{x}}^{i(n-i)}=\dfrac{\textbf{\textit{x}}^{n(n-1)/2}}{\textbf{\textit{x}}^{i(i-1)/2}\textbf{\textit{x}}^{(n-i)(n-i-1)/2}}$（特判 $\textbf{\textit{A}}=\textbf{1}$），故 $\{\textbf{\textit{f}}_n\}$ 为 **EGF 卷积**形式，做一次**卷积**即可。
然后做一次 **exp** 即可。时间复杂度 $O(n\log ^2 np(\mathbb F_{V}))$，其中 $p(V)$ 表示在 $V$ 内做 **Nim 积**的复杂度。
### Nimber EGF 怎么做？
先考虑卷积。
$$
\textbf{\textit{c}}_k=\sum_{i+j=k}\dbinom{k}{i}\textbf{\textit{a}}_i\textbf{\textit{b}}_j
$$
根据 Lucas 定理，你发现这个实际上就是要求了 $i$ 和 $j$ 的按位与为 $0$，求卷积。这就是子集卷积，可以在 $O(n\log ^2n)$ 的复杂度内完成。
至于 exp，直接套用牛顿迭代的公式即可。注意我们上面说过一个 **Nimber 序列** 的 **EGF 卷积逆**是它本身。EGF 的意义下，一个**多项式**求导就是**序列的值**平移。
### 实现细节
你发现你需要很快地求 Nim 积……吗？
一个很常见的 Nim 积求法是选择 $\mathbb F_{2^{64}}$ 的一个**本原元** $\textbf{\textit{g}}$。（或者说，其在 $\mathbb F_2$ 上的**最小多项式** $\textbf{\textit{p}}(\textbf{\textit{g}})=\textbf{0}$ 的次数为 $64$）则把 Nimber 同构到 $\mathbb F_2[X]/p(X)$ 上去，算出来乘法再同构回来。
这个求法是针对每次给定两个数求的。本题中，其实可以一开始就打到这个同构，最后输出的时候再同构回来！
那么怎么在 $\mathbb F_2[X]/p(X)$ 中快速运算呢？需要支持：
- 求两个 $<2^{64}$ 的数的不进位乘法。
- 求一个 $<2^{127}$ 的数对某个 $p$ 的不进位取模。
对于前者，可以发现指令集里有 `_mm_clmulepi64_si128` 就是用来算这个的！而后者，考虑取一个长得好一点的 $p(X)$（注意它需要在 $\mathbb F_2$ 上不可约），这里取了 $\textbf{\textit{p}}(\textbf{\textit{x}})=\textbf{\textit{x}}^{64}+\textbf{\textit{x}}^4+\textbf{\textit{x}}^3+\textbf{\textit{x}}+\textbf{1}$，这样取模只需要几次位运算就解决了。
具体的代码实现：
```cpp
inline u64 mult_gf264(u64 x,u64 y)
{
__m128i a=_mm_set_epi64x(0,x),b=_mm_set_epi64x(0,y);
__m128i c=_mm_clmulepi64_si128(a,b,0x00);
x=c[0];y=c[1];
u64 q=y^(y<<1),t=q>>61;q^=t^(t<<1);
return x^q^(q<<3);
}
```
等等我们怎么找这个同构？换言之我们要找到 $p$ 的根。求根的算法在论文里也讲了——首先 $\textbf{\textit{P}}(\textbf{\textit{x}})=(\textbf{\textit{x}}^{2^{64}}-\textbf{\textit{x}})$ 恰为所有 $\textbf{\textit{x}}-\textbf{\textit{a}},\textbf{\textit{a}}\in \mathbb F_{2^{64}}$ 的积——这是 Fermat 小定理。注意到取 $\textbf{\textit{Q}}(\textbf{\textit{x}})=\sum\limits_{i=0}^{63}\textbf{\textit{x}}^{2^i}$，则有 $\textbf{\textit{Q}}(\textbf{\textit{x}})(\textbf{\textit{Q}}(\textbf{\textit{x}})+\textbf{1})=\textbf{\textit{P}}(\textbf{\textit{x}})$。根据原论文的结论，随机取 $l(x)$ 求 $\gcd(\textbf{\textit{Q}}(\textbf{\textit{l}}(\textbf{\textit{x}})),\textbf{\textit{p}}(\textbf{\textit{x}}))$ 有不少于 $\dfrac12$ 的概率求出一个非平凡因式？我不知道为什么，但是它确实表现很好。
详解：[EI blog](https://www.cnblogs.com/Elegia/p/polynomial-factorization-over-finite-field.html)
贴一下解方程工具的代码：
```cpp
#include <bits/stdc++.h>
using namespace std;
using u64 = uint64_t;
using u128 = unsigned __int128;
const int B = 65536;
u64 g_pow[B * 2], g_lg[B];
u64 mult_q_init(u64 x, int lvl = 64) {
if (lvl == 1) return x;
int v = lvl >> 1;
if (x < (1ull << v)) return mult_q_init(x, v) << v;
u64 a = x >> v, b = x & ((1ull << v) - 1);
return (mult_q_init(a ^ b, v) << v) | mult_q_init(mult_q_init(a, v), v);
}
u64 mult_q(u64 x, int lvl = 64) {
if (x == 0) return 0;
if (lvl == 16) return g_pow[g_lg[x] + 40029];
int v = lvl >> 1;
if (x < (1ull << v)) return mult_q(x, v) << v;
u64 a = x >> v, b = x & ((1ull << v) - 1);
return (mult_q(a ^ b, v) << v) | mult_q(mult_q(a, v), v);
}
u64 mult_blk_init(u64 x, u64 y, int u = 16) {
if (x == 0 || y == 0) return 0;
if (x == 1 || y == 1) return x ^ y ^ 1;
int v = u / 2;
u64 a = x >> v, b = x & ((1ull << v) - 1), c = y >> v,
d = y & ((1ull << v) - 1);
if (a == 0) return mult_blk_init(b, c, v) << v | mult_blk_init(b, d, v);
if (c == 0) return mult_blk_init(a, d, v) << v | mult_blk_init(b, d, v);
u64 ac = mult_blk_init(a, c, v), bd = mult_blk_init(b, d , v),
apb_cpd = mult_blk_init(a ^ b, c ^ d, v);
u64 hi = apb_cpd ^ bd, lo = mult_q_init(ac, v) ^ bd;
return (hi << v) | lo;
}
u64 mult_blk(u64 x, u64 y, int u = 64) {
if (x == 0 || y == 0) return 0;
if (x < B && y < B) return g_pow[g_lg[x] + g_lg[y]];
if (x == 1 || y == 1) return x ^ y ^ 1;
int v = u / 2;
u64 a = x >> v, b = x & ((1ull << v) - 1), c = y >> v,
d = y & ((1ull << v) - 1);
if (a == 0) return mult_blk(b, c, v) << v | mult_blk(b, d, v);
if (c == 0) return mult_blk(a, d, v) << v | mult_blk(b, d, v);
u64 ac = mult_blk(a, c, v), bd = mult_blk(b, d , v),
apb_cpd = mult_blk(a ^ b, c ^ d, v);
u64 hi = apb_cpd ^ bd, lo = mult_q(ac, v) ^ bd;
return (hi << v) | lo;
}
void init(u64 g) {
g_pow[0] = 1;
for (int i = 1; i < B * 2; i++) g_pow[i] = mult_blk_init(g_pow[i - 1], g);
for (int i = 0; i < B - 1; i++)
if (!g_lg[g_pow[i]])
g_lg[g_pow[i]] = i;
else {
cout << g << " " << i << endl;
assert(0);
}
}
u64 qkpw(u64 a,u64 b)
{
u64 r=1;
while(b)
{
if(b&1)r=mult_blk(r,a);
a=mult_blk(a,a);
b>>=1;
}
return r;
}
u64 inv(u64 a)
{
return qkpw(a,u64(-2));
}
void unify(int N,u64 *g)
{
if(g[N]==1) return ;
u64 x=inv(g[N]);
for(int i=0;i<=N;++i) g[i]=mult_blk(g[i],x);
}
void mod(int M,u64 *g,int N,u64 *h)
{
unify(N,h);
for(int i=M;i>=N;--i)
{
if(g[i])
{
u64 cur=g[i];
for(int j=0;j<=N;++j)
g[j+i-N]^=mult_blk(cur,h[j]);
}
}
}
int find_degree(int N,u64 *g)
{
int cur=N;while(cur>=0&&g[cur]==0)--cur;
return g[cur]?cur:(-1);
}
void mult(int M,u64 *g,int N,u64 *h,u64 *ans)
{
for(int i=0;i<=M+N;++i) ans[i]=0;
for(int i=0;i<=M;++i)
{
for(int j=0;j<=N;++j)
{
ans[i+j]^=mult_blk(g[i],h[j]);
}
}
}
const int Bq=5;
mt19937_64 myrand(chrono::steady_clock::now().time_since_epoch().count());
u64 rand64() {
return myrand();
}
int GCD(int N,u64 *g,int M,u64 *h)
{
unify(N,g);unify(M,h);
while(1)
{
if(N>=0&&M>=0)
{
if(N>=M)
{
mod(N,g,M,h);
N=find_degree(M-1,g);
unify(N,g);
}
else
{
mod(M,h,N,g);
M=find_degree(N-1,h);
unify(M,h);
}
}
else
{
if(N==-1)
{
N=M;for(int i=0;i<=M;++i) g[i]=h[i];
}
else
{
M=N;for(int i=0;i<=N;++i) h[i]=g[i];
}
return M;
}
}
}
void solve(int N,u64 *f)
{
if(N==0) return ;
unify(N,f);
if(N==1)
{
cout<<f[0]<<'\n';
return ;
}
u64 l[191],ret[191],qaq[191];
int d;
int cnt=0;
while(1)
{
for(int i=0;i<=Bq;++i) l[i]=rand64();
if(l[Bq]==0) l[Bq]=1;
for(int i=0;i<N;++i) ret[i]=0;
for(int __=0;__<64;++__)
{
for(int i=0;i<N;++i) ret[i]^=l[i];
mult(N-1,l,N-1,l,qaq);
mod(2*N-2,qaq,N,f);
for(int i=0;i<N;++i) l[i]=qaq[i];
}
for(int i=0;i<=N;++i) qaq[i]=f[i];
d=GCD(find_degree(N-1,ret),ret,N,qaq);
if(d!=0&&d!=N) break;
if((++cnt)>100) return ;
}
u64 a[191],b[191];
for(int i=0;i<=d;++i) a[i]=qaq[i];
for(int i=N;i>=d;--i)
{
b[i-d]=f[i];
if(f[i])
{
u64 cur=f[i];
for(int j=0;j<=d;++j)
f[j+i-d]^=mult_blk(cur,a[j]);
}
}
solve(d,a);solve(N-d,b);
}
int n;
u64 f[1145];
int main() {
init(2025);
cin>>n;
for(int i=0;i<=n;++i) cin>>f[i];
solve(n,f);
return 0;
}
```
求出来可以取的一个根是 $\textbf{4\,928\,496\,685\,556\,603\,065}$。
还有一些卡常细节，比如为了避免 cache miss 数组要开大一点，具体看代码吧。
### 代码
```cpp
#include <bits/stdc++.h>
#include <immintrin.h>
#include <emmintrin.h>
using namespace std;
#define popc(i) popcount(unsigned(i))
using u64 = uint64_t;
using u128 = unsigned __int128;
const int B = 65536;
u64 g_pow[B * 2], g_lg[B];
u64 mult_q_init(u64 x, int lvl = 64)
{
if (lvl == 1) return x;
int v = lvl >> 1;
if (x < (1ull << v)) return mult_q_init(x, v) << v;
u64 a = x >> v, b = x & ((1ull << v) - 1);
return (mult_q_init(a ^ b, v) << v) | mult_q_init(mult_q_init(a, v), v);
}
u64 mult_q(u64 x, int lvl = 64)
{
if (x == 0) return 0;
if (lvl == 16) return g_pow[g_lg[x] + 40029];
int v = lvl >> 1;
if (x < (1ull << v)) return mult_q(x, v) << v;
u64 a = x >> v, b = x & ((1ull << v) - 1);
return (mult_q(a ^ b, v) << v) | mult_q(mult_q(a, v), v);
}
u64 mult_blk_init(u64 x, u64 y, int u = 16)
{
if (x == 0 || y == 0) return 0;
if (x == 1 || y == 1) return x ^ y ^ 1;
int v = u / 2;
u64 a = x >> v, b = x & ((1ull << v) - 1), c = y >> v,
d = y & ((1ull << v) - 1);
if (a == 0) return mult_blk_init(b, c, v) << v | mult_blk_init(b, d, v);
if (c == 0) return mult_blk_init(a, d, v) << v | mult_blk_init(b, d, v);
u64 ac = mult_blk_init(a, c, v), bd = mult_blk_init(b, d , v),
apb_cpd = mult_blk_init(a ^ b, c ^ d, v);
u64 hi = apb_cpd ^ bd, lo = mult_q_init(ac, v) ^ bd;
return (hi << v) | lo;
}
u64 mult_blk(u64 x, u64 y, int u = 64)
{
if (x == 0 || y == 0) return 0;
if (x < B && y < B) return g_pow[g_lg[x] + g_lg[y]];
if (x == 1 || y == 1) return x ^ y ^ 1;
int v = u / 2;
u64 a = x >> v, b = x & ((1ull << v) - 1), c = y >> v,
d = y & ((1ull << v) - 1);
if (a == 0) return mult_blk(b, c, v) << v | mult_blk(b, d, v);
if (c == 0) return mult_blk(a, d, v) << v | mult_blk(b, d, v);
u64 ac = mult_blk(a, c, v), bd = mult_blk(b, d , v),
apb_cpd = mult_blk(a ^ b, c ^ d, v);
u64 hi = apb_cpd ^ bd, lo = mult_q(ac, v) ^ bd;
return (hi << v) | lo;
}
void init(u64 g) {
g_pow[0] = 1;
for (int i = 1; i < B * 2; i++) g_pow[i] = mult_blk_init(g_pow[i - 1], g);
for (int i = 0; i < B - 1; i++)
if (!g_lg[g_pow[i]])
g_lg[g_pow[i]] = i;
else {
cout << g << " " << i << endl;
assert(0);
}
}
u64 qkpw(u64 a,u64 b)
{
u64 r=1;
while(b)
{
if(b&1)r=mult_blk(r,a);
a=mult_blk(a,a);
b>>=1;
}
return r;
}
u64 inv(u64 a)
{
return qkpw(a,u64(-2));
}
u64 pw[1145],bas[1145],val[1145];
void init2(u64 g)
{
pw[0]=1;for(int i=1;i<128;++i)pw[i]=mult_blk(pw[i-1],g);
for(int i=0;i<64;++i)
{
u64 cur=pw[i],vl=1ull<<i;
for(int j=63;j>=0;--j)
{
if((cur>>j)&1)
{
if(bas[j])
{
cur^=bas[j];vl^=val[j];
}
else
{
bas[j]=cur;val[j]=vl;
break;
}
}
}
}
for(int j=0;j<=63;++j)
{
for(int k=0;k<j;++k) if((bas[j]>>k)&1)
{
bas[j]^=bas[k];
val[j]^=val[k];
}
}
}
u64 nimber_to_gf264(u64 x)
{
u64 ret=0;
for(int i=0;i<64;++i) if((x>>i)&1) ret^=val[i];
return ret;
}
inline u64 mult_gf264(u64 x,u64 y)
{
__m128i a=_mm_set_epi64x(0,x),b=_mm_set_epi64x(0,y);
__m128i c=_mm_clmulepi64_si128(a,b,0x00);
x=c[0];y=c[1];
u64 q=y^(y<<1),t=q>>61;q^=t^(t<<1);
return x^q^(q<<3);
}
u64 gf264_to_nimber(u64 x)
{
u64 ret=0;
for(int i=0;i<64;++i) if((x>>i)&1) ret^=pw[i];
return ret;
}
const int Bt=1048576;
u64 p[21][Bt+256],q[21][Bt+256],r[21][Bt+256];
void conv(int k,u64 *a,u64 *b,u64 *ans)
{
int n=1<<k;
for(int i=0;i<=k;++i)
{
memset(p[i],0,sizeof(u64)<<k);
memset(q[i],0,sizeof(u64)<<k);
memset(r[i],0,sizeof(u64)<<k);
}
for(int i=0;i<n;++i)
{
p[popc(i)][i]=a[i];
q[popc(i)][i]=b[i];
}
for(int l=0;l<k;++l)
for(int j=0;j<n;++j) if((j>>l)&1)
{
for(int i=0;i<k;++i)
{
p[i][j]^=p[i][j^(1<<l)];
q[i][j]^=q[i][j^(1<<l)];
}
}
for(int i=0;i<=k;++i)
{
for(int j=0,l=i;l>=0;++j,--l)
{
for(int id=0;id<n;++id)
r[i][id]^=mult_gf264(p[j][id],q[l][id]);
}
}
for(int l=0;l<k;++l)
for(int j=0;j<n;++j) if((j>>l)&1)
{
for(int i=0;i<=k;++i)
{
r[i][j]^=r[i][j^(1<<l)];
}
}
for(int i=0;i<n;++i) ans[i]=r[popc(i)][i];
}
void conv_small_large(int k,u64 *b,u64 *ans)
{
int n=1<<k;
for(int i=0;i<=k;++i)
{
memset(q[i],0,sizeof(u64)<<k);
memset(r[i],0,sizeof(u64)<<k);
}
for(int i=0;i<n;++i)
{
q[popc(i)][i]=b[i];
}
for(int l=0;l<k;++l)
for(int j=0;j<n;++j) if((j>>l)&1)
{
for(int i=0;i<k;++i)
{
q[i][j]^=q[i][j^(1<<l)];
}
}
for(int i=0;i<=k;++i)
{
for(int j=0,l=i;l>=0;++j,--l)
{
for(int id=0;id<n/2;++id)
r[i][id]^=mult_gf264(p[j][id],q[l][id]);
for(int id=n/2;id<n;++id)
r[i][id]^=mult_gf264(p[j][id-n/2],q[l][id]);
}
}
for(int l=0;l<k;++l)
for(int j=0;j<n;++j) if((j>>l)&1)
{
for(int i=0;i<=k;++i)
{
r[i][j]^=r[i][j^(1<<l)];
}
}
for(int i=0;i<n;++i) ans[i]=r[popc(i)][i];
}
u64 qaq[Bt+256];
void exp(int k,u64 *f,u64 *ans)
{
ans[0]=1;
for(int i=1;i<=k;++i)
{
ans[1<<(i-1)]=0;
conv(i-1,ans,ans+1,qaq+1);
qaq[0]=1;
for(int j=0;j<(1<<i);++j) qaq[j]^=f[j];
conv_small_large(i,qaq,ans);
}
}
int paranoia,k,n;
u64 a,iva;
u64 pia[Bt+256],piab[Bt+256];
u64 x[Bt+256],y[Bt+256],z[Bt+256];
u64 f[Bt+256];
u64 ans[Bt+256];
u64 pwinva(u64 t)
{
return mult_gf264(pia[t&1048575],piab[t>>20]);
}
int main()
{
init(2025);
init2(4928496685556603065);
#ifndef DEBUG
ios_base::sync_with_stdio(false);cin.tie(0);
#endif
cin>>paranoia>>a;
a^=1;iva=inv(a);
k=1;
while((1<<k)<=paranoia) ++k;
n=1<<k;
a=nimber_to_gf264(a);iva=nimber_to_gf264(iva);
if(a==0)
{
for(int i=1;i<n;++i) f[i]=1;
}
else
{
pia[0]=1;for(int i=1;i<=Bt;++i) pia[i]=mult_gf264(pia[i-1],iva);
piab[0]=1;for(int i=1;i<=Bt;++i) piab[i]=mult_gf264(piab[i-1],pia[Bt]);
for(int i=0;i<n;++i)
{
x[i]=pwinva(i*(i+1ull)/2);
y[i]=pwinva(i*(i-1ull)/2);
}
conv(k,x,y,z);
u64 mult=1,qaaq=1;
for(int i=1;i<n;++i)
{
f[i]=mult_gf264(z[i-1],qaaq);
mult=mult_gf264(mult,a);
qaaq=mult_gf264(qaaq,mult);
}
}
exp(k,f,ans);
for(int i=0;i<=paranoia;++i)
{
cout<<gf264_to_nimber(ans[i])<<" \n"[i==paranoia];
}
return 0;
}
```
## 致谢
感谢 @[zhuzhu2891](https://www.luogu.com.cn/user/515385) 和 @[Galois_Field_1048576](https://www.luogu.com.cn/user/360265) 对我解决本题提供的帮助。
感谢 UOJ 群群友告诉我 AtCoder 的 exit code 9 不是 RE 是 TLE。
参考文献：
《浅谈 Nimber 和多项式算法》宁波市镇海中学 罗煜翔