跳转至

Welcome to Liam's Blog

Pollard-Rho 算法

推导

对于一个已知的合数 \(n>3\),我们的目标是快速找到它的一个非平凡因子(除 \(1\)\(n\) 外)。

显然问题可以转化为找到一个 \(k\) 使得 \(1 \lt \gcd(k,n) \lt n\)

考虑一个序列 \(f\),满足 \(f_0=0\)\(f_i = (f_{i-1}^2 + c) \bmod n\),其中 \(c\) 是常数。由生日悖论1\(f\) 中不同的数量的期望为 \(O(\sqrt n)\)

\(m\)\(n\) 的最小质因子,显然有 \(m \le \sqrt n\)

\(g_i = f_i \bmod m\),则 \(g\) 中不同的数量的期望为 \(O(\sqrt m)\)

于是我们可以在期望 \(O(\sqrt m) \le O(n ^ {\frac 1 4})\) 的时间内找到两个位置 \(i,j\) 使得 \(f_i \not = f_j \land g_i = g_j\),即 \(n \nmid (f_i - f_j) \land m \mid (f_i - f_j)\),于是有 \(1 \lt m \le \gcd(|f_i - f_j|, n) \lt n\)

于是我们可以在期望 \(O(n ^ {\frac 1 4})\) 的时间内找到 \(k\) 满足 \(1 \lt gcd(k,n) \lt n\),也就找到了一个 \(n\) 的非平凡因子。

实现

我们随机 \(c \in [1,n)\)

Floyd 判环

由于 \(f\) 的值域在 \([0,n)\) 之间,于是其必定有循环节,因为其值的排列图像形似 \(\rho\),于是这个算法被称为 Pollard-Rho

我们考虑枚举 \(i\),判断是否 \(1 \lt \gcd(|f_i - f_{2i}|,n) \lt n\),并在当 \(f_i = f_{2i}\) 时终止。

code
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
ll pollard_rho(ll n) {
    const ll c = randint(1,n-1);
    const auto f = [&](ll x) -> ll {
        return (mul(x,x,n)+c)%n;
    };
    ll s=f(0),t=f(f(0));
    while(true) {
        if(s==t) break;
        ll d=gcd(abs(s-t),n);
        if(d>1) return d;
        s=f(s);
        t=f(f(t));
    }
    return n;
}

倍增思想

考虑设两个变量 \(s,t\),我们让 \(t\) 记下 \(s\) 的目前值,并让 \(s\) 跑一定距离,在这过程中检验 \(|s-t|\),并每次都把距离翻倍,显然当 \(s\) 跑到 \(t\),即出现环时,\(s\) 并没有多跑多远。

code
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
ll pollard_rho(ll n) {
    const ll c = randint(1,n-1);
    const auto f = [&](ll x) -> ll {
        return (mul(x,x,n)+c)%n;
    };
    ll s=0,t=0;
    for(int j=1; ; j<<=1,t=s)
        for(int i=1; i<=j; ++i) {
            s=f(s);
            if(s==t) return n;
            ll d=gcd(abs(s-t),n);
            if(d>1) return d;
        }
    return n;
}

基于倍增的二次优化

我们发现上面两种做法时间复杂度都不是单纯的 \(O(n ^ {\frac 1 4})\),而是 \(O(n ^ {\frac 1 4} \log n)\),因为还要调用 \(\gcd\)

为了减小调用 \(gcd\) 的次数,我们考虑若 \(\gcd(a,b) \gt 1\),显然也有 \(\gcd(ac \bmod b,b) \gt 1\),于是我们把一段路程放到一起检验。

根据前人的智慧,取段长为 \(127\) 效果最好。

code
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
ll pollard_rho(ll n) {
    const ll c = randint(1,n-1);
    const auto f = [&](ll x) -> ll {
        return (mul(x,x,n)+c)%n;
    };
    ll s=0,t=0,v=1;
    for(int j=1; ; j<<=1,t=s,v=1) {
        for(int i=1; i<=j; ++i) {
            s=f(s);
            v=mul(v,abs(s-t),n);
            if(!v) return n;
            if(i%127==0) {
                ll d=gcd(v,n);
                if(d>1) return d;
            }
        }
        ll d=gcd(v,n);
        if(d>1) return d;
    }
    return n;
}

完整代码

P4718 【模板】Pollard-Rho
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
#include <cstdio>
#include <random>
#include <chrono>
#include <algorithm>
#include <cmath>
using namespace std;
using ll = long long;
int T;
ll n,maxf;
mt19937 rnd(chrono::steady_clock::now().time_since_epoch().count());
ll randint(ll l,ll r) {
    return uniform_int_distribution<ll>(l,r)(rnd);
}
ll mul(ll a,ll b,ll p) {
    ll c=1.0l*a/p*b;
    return (ll(1ull*a*b-1ull*c*p)+p)%p;
}
ll qpow(ll a,ll n,ll p) {
    a%=p;
    ll x=1;
    for(; n; n>>=1,a=mul(a,a,p)) if(n&1) x=mul(x,a,p);
    return x;
}
bool miller_rabin(ll n) {
    if(n<3) return n==2;
    if(n%2==0) return false;
    ll s=n-1,t=0;
    while(s%2==0) s/=2,++t;
    for(int i=0; i<8; ++i) {
        ll x=randint(2,n-1),u=qpow(x,s,n);
        if(u==1) continue;
        {
            int j;
            for(j=0; j<t; ++j) {
                if(u==n-1) break;
                u=mul(u,u,n);
            }
            if(j==t) return false;
        }
    }
    return true;
}
ll gcd(ll a,ll b) {
    return !b?a:gcd(b,a%b);
}
ll pollard_rho(ll n) {
    const ll c = randint(1,n-1);
    const auto f = [&](ll x) -> ll {
        return (mul(x,x,n)+c)%n;
    };
    ll s=0,t=0,v=1;
    for(int j=1; ; j<<=1,t=s,v=1) {
        for(int i=1; i<=j; ++i) {
            s=f(s);
            v=mul(v,abs(s-t),n);
            if(!v) return n;
            if(i%127==0) {
                ll d=gcd(v,n);
                if(d>1) return d;
            }
        }
        ll d=gcd(v,n);
        if(d>1) return d;
    }
    return n;
}
void split(ll x) {
    if(x<=maxf||x<2||miller_rabin(x)) {
        maxf=max(maxf,x);
        return ;
    }
    ll d=x;
    while(d==x) d=pollard_rho(x);
    while(x%d==0) x/=d;
    split(x);
    split(d);
}
int main() {
    scanf("%d",&T);
    while(T--) {
        scanf("%lld",&n);
        maxf=0;
        split(n);
        if(maxf==n) printf("Prime\n");
        else printf("%lld\n",maxf);
    }
    return 0;
}

  1. 请读者自行阅读其他材料 

CF1267H Help BerLine

题意

已知有 \(n \le 8500\) 个基地,以及开启它们的顺序排列 \(p\),初始所有基地都关闭。

求一个整数序列 \(f\) 使得在任意时刻,已开启的基地的 \(f\) 值按顺序排成的序列满足:

  • 任意子段都有一个只出现在其中过一次的值(以下记为合法

  • \(\forall i, f_i \in [1,24]\)

思考

发现 \(\lceil \log _{\frac 3 2} 8500 \rceil = 23\),于是考虑缩小问题规模递归求解。

考虑用一种颜色 \(c\) 尽可能多且合法地填。

首先合法的必要条件显然是任意时刻序列中相邻的值都不同。

发现如果填 \(c\) 的位置在任意时刻都不相邻,且剩下的位置提取出来后任意时刻合法,那么必有当前序列任意时刻合法。

形式化地,设 \(S \subseteq \{1,2,\dots,n\}\) 为填 \(c\) 的位置集合,即 \(\forall i \in [1,n], f_i = c \iff i \in S\)

那么如果任意时刻均有 \(\forall i,j \in S \land i \not = j, i,j\) 不相邻,便可递归解决剩下的位置集合 \(\{1,2,\dots,n\} \setminus S\)

感性理解:对于任意时刻,显然任意值为 \(c\) 的位置不相邻,那么便如下图:

显然任意子段都属于上面三种之一,且都合法。

解法

于是我们便可以尽可能地在满足任意时刻填 \(c\) 的位置都不相邻的前提下,填尽可能多的 \(c\)

于是问题就简单了,任意时刻填 \(c\) 的位置都不相邻可以转化为当它刚好被加入序列的时刻它的前驱后继都不为 \(c\),于是从后往前扫一遍就可以了。

且显然有一个点填 \(c\) 最多使得两个点无法填 \(c\),于是可以把问题规模减小 \(\lceil \frac n 3 \rceil\),可以通过本题。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
void _main() {
    int n;
    scanf("%d",&n);
    vec<int> p(n),s(n);
    for(auto &i : p) scanf("%d",&i),--i;
    reverse(p.begin(),p.end());
    for(int i=1; !p.empty(); ++i) {
        set<int> st;
        vec<bool> v(n,0);
        for(auto u : p) st.insert(u);
        for(auto u : p) {
            if(!v[u]) {
                v[u]=true;
                s[u]=i;
                auto it=st.find(u);
                if(it!=st.begin()) {
                    auto j=it;
                    --j;
                    v[*j]=true;
                }
                {
                    auto j=it;
                    ++j;
                    if(j!=st.end()) v[*j]=true;
                }
            }
            st.erase(u);
        }
        vec<int> q;
        for(auto u : p) if(!s[u]) q.push_back(u);
        p=q;
    }
    for(auto i : s) printf("%d ",i);
    putchar('\n');
}

补充

这个做法十分的不直观,因为考虑 \(p = \{1,2,3,4\}\),若 \(f = \{1,2,1,2\}\),虽然符合任意时刻相邻两个都不相同,但却是非法的。

这是因为当我们把 \(\{1,3\}\) 两个位置的 \(f\) 标为 \(1\) 时,缩小了问题规模,删除了他们,这时 \(p = \{2,4\}\),因此 \(2\)\(4\) 相邻,我们就不能使 \(f = \{2,2\}\),而是 \(f = \{3,2\}\),合并得到 \(f = \{1,3,1,2\}\)

FFT & NTT 复习笔记

默认文中的形如 \([l,r)\) 的区间为其与整数集的交集。

快速变换

设原多项式为 \(F(x) = \sum_{i \in [0,n)} a_i x ^ i\),其中 \(n = 2 ^ k, k \in \mathbb Z ^ +\)

我们要求 \(\forall i \in [0,n),\hat a_i = F(t_i)\),其中 \(t\) 是一个长度为 \(n\) 且两两互不相同的序列。

显然 \(F\) 可以被一组 \(\hat a,t\) 唯一确定,即点值表示法。


另设两个多项式

\[ G_0(x)=a_0 + a_2 x + \dots + a_{n - 2} x^{\frac n 2 - 1} \\ G_1(x)=a_1 + a_3 x + \dots + a_{n - 1} x^{\frac n 2 - 1} \\ \]

则有

\[ \begin{aligned} F(x) & = \sum_{i \in [0,n)} a_i x ^ i \\ & = a_0 + a_1 x + a_2 x^2 + \dots + a_{n-1} x^{n-1} \\ & = (a_0 + a_2 x^2 + \dots + a_{n-2} x^{n-2}) + (a_1 x + a_3 x^3 + \dots + a_{n-1} x^{n-1}) \\ & = G_0 (x ^ 2) + x G_1 (x ^ 2) \\ \end{aligned} \]

考虑构造单位根 \(\omega _n ^k\) 满足 \(\omega _n ^{\frac n 2} = -1, \omega _{2n} ^ {2k} = \omega _n ^k\)

显然也有 \(\omega _n ^n = \omega _n ^0 = 1\)

\(\forall i \in [0,n), t_i = \omega _n ^i\)

\(x = \omega _n ^k, k \in [0,\frac n 2)\) 时显然有

\[ \begin{aligned} F(\omega _n ^k) & = G_0(\omega _n ^{2k}) + \omega _n ^k G_1(\omega _n ^{2k}) \\ & = G_0(\omega _ {\frac n 2} ^ k) + \omega _n ^k G_1(\omega _{\frac n 2} ^k) \\ \end{aligned} \]

\(x = \omega _n ^{k + \frac n 2}, k \in [0,\frac n 2)\) 时有

\[ \begin{aligned} F(\omega _n ^{k + \frac n 2}) & = G_0(\omega _n ^{2k + n}) + \omega _n ^{k + \frac n 2} G_1(\omega _n ^{2k + n}) \\ & = G_0(\omega _n ^{2k} \cdot \omega _n ^n) - \omega _n ^k G_1(\omega _n ^{2k} \cdot \omega _n ^n) \\ & = G_0(\omega _n ^{2k}) - \omega _n ^k G_1(\omega _n ^{2k}) \\ & = G_0(\omega _{\frac n 2} ^k) - \omega _n ^k G_1(\omega _{\frac n 2} ^k) \\ \end{aligned} \]

由于两者只有一个符号的差异,于是 \(F\) 的点值可以直接 \(\mathrm O(n)\)\(G_0, G_1\) 的点值得到。

递归解决,时间复杂度 \(\mathrm O(n \log n)\)

逆变换

设变换后的点值序列为 \(\hat a\),即

\[ \begin{aligned} \forall i \in [0,n), \hat a_i & = F(\omega _n ^i) \\ & = \sum _{j \in [0,n)} a_j (\omega _n ^i)^j \\ & = \sum _{j \in [0,n)} a_j \omega _n ^{ij} \\ \end{aligned} \]

设多项式 \(\hat F(x) = \sum _{i \in [0,n)} \hat a_i x^i\)

\(\hat F\) 进行点值变换(\(\forall i \in [0,n),t_i = \omega _n ^{-i}\)),设点值序列为 \(s\)

则有

\[ \begin{aligned} \forall i \in [0,n), s_i & = \hat F(\omega _n ^{-i}) \\ & = \sum _{j \in [0,n)} \hat a_j (\omega _n ^{-i}) ^j \\ & = \sum _{j \in [0,n)} \omega _n ^{-ij} \hat a_j \\ & = \sum _{j \in [0,n)} \omega _n ^{-ij} \sum _{k \in [0,n)} a_k \omega _n ^{jk} \\ & = \sum _{j \in [0,n), k \in [0,n)} \omega _n ^{-ij} a_k \omega _n ^{jk} \\ & = \sum _{j \in [0,n), k \in [0,n)} \omega _n ^{j(k-i)} a_k \\ & = \sum _{k \in [0,n)} a_k \sum _{j \in [0,n)} \omega _n ^{j(k-i)} \\ & = \sum _{k \in [0,n)} a_k \sum _{j \in [0,n)} (\omega _n ^{k-i}) ^j \\ \end{aligned} \]

显然第二个求和是一个等比数列,由等比数列求和公式 \(\sum _{i \in [m,n)} p^i = \frac {p^m - p^n} {1 - p}\) 得:

  • \(\omega _n ^{k-i} \not = 1 \iff i \not = k\)
\[ \begin{aligned} \sum _{j \in [0,n)} (\omega _n ^{k-i}) ^j & = \frac {1 - \omega _n ^{(k-i) n}} {1 - \omega _n ^{k-i}} \\ & = \frac {1 - (\omega _n ^{k-i}) ^n} {1 - \omega _n ^{k-i}} \\ & = \frac {1 - 1} {1 - \omega _n ^{k-i}} \\ & = 0 \end{aligned} \]
  • \(\omega _n ^{k-i} = 1 \iff i = k\)
\[ \sum _{j \in [0,n)} (\omega _n ^{k-i}) ^j = \sum _{j \in [0,n)} 1 = n \]

因此

\[ \begin{aligned} \forall i \in [0,n), s_i & = \sum _{k \in [0,n)} a_k \sum _{j \in [0,n)} (\omega _n ^{k-i}) ^j \\ & = n a_i \\ \end{aligned} \]

于是我们有

\[ \forall i \in [0,n), a_i = \frac {s_i} n \]

构造单位根

  • FFT

在复数域下,有 \(\omega _n = \cos \frac {2 \pi} n + \mathrm i \sin \frac {2 \pi} {n}\)

其中 \(\mathrm i = \sqrt {-1}\) 是 虚数单位,可以用 C++ 中的 complex 库中的 std::complex<double/long double> 存储复数。

  • NTT

对于模数 \(P \in \mathbb P, \exists n,k \in \mathbb Z^+, P=2^nk+1\),在模 \(P\) 意义下有 \(\omega _n \equiv g ^ {\frac {P-1} n}\),其中 \(g\) 是原根。

\(g\) 是模 \(P\) 意义下的原根当且仅当 \(g ^i \not \equiv 1 \pmod P,\forall i \in [1,\phi(P))\)\(g ^{\phi(P)} \equiv 1 \pmod P\)

specially,\(\forall P \in \mathbb P\),其原根 \(g\) 满足 \(\forall i \in [1,P-1), g ^i \not \equiv 1 \pmod P\)\(g^{P-1} \equiv 1 \pmod P\)

于是对 \(n = 2 ^m, m \in \mathbb Z ^+\),我们有 \(\omega _n ^n \equiv g ^{\frac {P - 1} {n} \cdot n} \equiv g ^{P - 1} \equiv 1, \pmod P\),且 \(\omega _n ^{\frac n 2} \equiv g ^{\frac {P - 1} 2} \equiv \pm \sqrt {g ^ {P - 1}} \equiv \pm 1 \pmod P\),又 \(g ^ {\frac {P-1} 2} \not \equiv 1 \pmod P\),所以 \(\omega _n ^{\frac n 2} \equiv -1 \pmod P\)

还有 \(\omega _{2n} ^{2k} \equiv g ^{\frac {2k(P-1)} {2n}} \equiv g ^{\frac{k(P-1)} n} \equiv \omega _n ^k \pmod P\)

由于原根的特殊性,模数 \(P \in \mathbb P\) 有特殊的限制,一般有 \(P = k 2 ^m + 1, k,m\in \mathbb Z ^+\)

常见的模数有

\[ \begin{aligned} 167772161 = 5 \times 2 ^{25} + 1, g = 3 \\ 469762049 = 7 \times 2 ^{26} + 1, g = 3 \\ 754974721 = 45 \times 2 ^{24} + 1, g = 11 \\ 998244353 = 119 \times 2 ^{23} + 1, g = 3 \\ 1004535809 = 479 \times 2 ^{21} + 1, g = 3 \\ \end{aligned} \]

平面直角坐标系的点绕原点旋转公式及证明

设点 \(A(x,y)\) 绕原点 \(O(0,0)\) 逆时针旋转 \(\beta\),则设在极坐标系下 \(A\) 的坐标为 \((r,\alpha)\)

这意味着 \(x=r \cos \alpha, y=r \sin \alpha\)

目标点 \(A'(x',y')\) 的极坐标即为 \((r,\alpha + \beta)\)

展开(其中 \(\sin\)\(\cos\) 的展开参考 here):


\[ \begin{aligned} x' & = r \cos (\alpha + \beta) \\ & = r(\cos \alpha \cos \beta - \sin \alpha \sin \beta) \\ & = r\cos \alpha \cos \beta - r\sin \alpha \sin \beta \\ & = x \cos \beta - y \sin \beta \\ \end{aligned} \]

\[ \begin{aligned} y' & = r \sin (\alpha + \beta) \\ & = r (\sin \alpha \cos \beta + \cos \alpha \sin \beta) \\ & = r \sin \alpha \cos \beta + r\cos \alpha \sin \beta \\ & = y \cos \beta + x \sin \beta\\ \end{aligned} \]

结论:点 \((x,y)\) 绕原点逆时针旋转 \(\theta\) 后坐标为 \((x \cos \theta - y \sin \theta, x \sin \theta + y \cos \theta)\)

\(\sin/\cos(\alpha+\beta)\) 的展开证明


\[ \begin{aligned} \cos(α+β) &= OB \\ & = OD - BD \\ & = OD - EC \\ & = OC \cos \beta - AC \sin \beta \\ & = OA \cos \alpha \cos \beta - OA \sin \alpha \sin \beta \\ & = \cos \alpha \cos \beta - \sin \alpha \sin \beta \end{aligned} \]

\[ \begin{aligned} \sin(α+β) &= AB \\ & = AE + BE \\ & = AE + CD \\ & = AC \cos \beta + OC \sin \beta \\ & = OA \sin \alpha \cos \beta + OA \cos \alpha \sin \beta \\ & = \sin \alpha \cos \beta + \cos \alpha \sin \beta \\ \end{aligned} \]

\(\rm FWT\) (Fast Walsh Transform) 学习笔记

本文中使用 \(\cap\) 表示按位与,用 \(\cup\) 表示按位或

与/或 卷积

问题引入

给定长度为 \(2^n\) 的数列 \(A,B\),求 \(C_i = \sum_{j \cup k = i} A_j \times B_k\)

显然有 \(O(4^n)\) 的暴力

正变换

这一部分可以参考 FMT 中的 Zeta 变换 ,即定义 \(\hat A_i\) 为序列 \(A\)\(i\) 的子集和

快速变换

考虑分治

对于长度为 \(2^n\) 的待求解的 \(\hat A\),我们把它分成独立的两部分,即 \([0,2^{n-1})\)\([2^{n-1},2^n)\) 两个区间分别求解

然后考虑这两个区间之间的贡献,发现只有 \([0,2^{n-1})\)\([2^{n-1},2^n)\) 有贡献,于是对于 \(i \in [0,2^{n-1})\),将 \(\hat A_{i}\) 加到 \(\hat A_{i+2^{n-1}}\) 中即可

时间复杂度 \(O(n 2^n)\)

逆变换

我们只需要对照正变换中的操作步骤,一步一步撤销变换即可

即对于 \(i \in [0,2^{n-1})\),将 \(\hat A_{i+2^{n-1}}\) 减去 \(\hat A_{i}\) 的贡献,然后再递归地求解


当然还有另外的做法,即我们要将长度为 \(2^n\)\(\hat A_i\) 在全集为 \(2^n-1\) 的情况下(即不考虑目前的长度外的子集)只包括自己的 \(A_i\),那么我们递归地求解完左右两个区间后,肯定有 \(\forall i\in [0,2^{n-1}),\hat A_{i+2^{n-1}}=\hat A_i+A_{i+2^{n-1}}\),因此减去贡献即可

in brief,可以先递归再减贡献,这样逆变换与正变换只有一个 +/- 的变化

推广

对于 \(\cap\) 的情况,与 \(\cup\) 十分相似,请读者尽量自己构造变换,推一下式子,可以加深理解

如果没有思路,here

异或 卷积

问题引入

给定长度为 \(2^n\) 的数列 \(A,B\),求 \(C_i = \sum_{j \oplus k = i} A_j \times B_k\)

原理

\(\operatorname{popcnt}(i)\) 表示 \(i\) 在二进制下 1 的个数,则有

\[ \operatorname{popcnt}(i \cap k) + \operatorname{popcnt}(j \cap k) \equiv \operatorname{popcnt}((i\oplus j)\cap k) \pmod 2 \]

证明:显然 \(k\)\(0\) 的位上,\(i,j\) 的取值不影响结果,那么设 \(i'=i\cap k,j'=j\cap k\),那么问题转化为 $$ \operatorname{popcnt}(i') + \operatorname{popcnt}(j') \equiv \operatorname{popcnt}(i'\oplus j') \pmod 2 $$ 对于 \(i',j'\) 每一位分开讨论,原命题易证

即异或不会改变 \(1\) 的总数的奇偶性

正变换

我们尝试构造一个变换 $$ \hat A_i = \sum_j g(i,j) A_j $$ 例如对于 \(\cup\) 卷积 ,\(g(i,j)=[i\cup j=i]\)

这个 \(g\) 函数满足以下性质:

\[ \begin{aligned} & \because \hat C_i = \hat A_i \times \hat B_i \\ & \therefore \sum_{p} g(i,p)C_p = \sum_{j,k} g(i,j)\times g(i,k)\times A_j \times B_k \\ & \therefore \sum_{p} g(i,p) \sum_{j \oplus k = p} A_j \times B_k = \sum_{j,k} g(i,j)\times g(i,k)\times A_j \times B_k \\ & \therefore \sum_{j,k} g(i,j \oplus k) A_j \times B_k = \sum_{j,k} g(i,j)\times g(i,k)\times A_j \times B_k \end{aligned} \]

即我们要使得 \(g(i,j)\times g(i,k) = g(i,j \oplus k)\)

因为 \(\operatorname{popcnt}(i \cap k) + \operatorname{popcnt}(j \cap k) \equiv \operatorname{popcnt}((i\oplus j)\cap k) \pmod 2\),我们发现 \(g(i,j) = (-1)^{\operatorname{popcnt}(i \cap j)}\) 满足这个性质

\[ \hat A_i = \sum_j (-1)^{\operatorname{popcnt}(i \cap j)} A_j \]

手动推一下式子:

\[ \begin{aligned} \hat A_i \times \hat B_i & = \sum_j g(i,j)A_j \times \sum_k g(i,k) B_k\\ & = \sum_j (-1)^{\operatorname{popcnt}(i \cap j)} A_j \times \sum_k (-1)^{\operatorname{popcnt}(i \cap k)} B_k\\ & = \sum_{j,k} (-1)^{\operatorname{popcnt}(i \cap j) + \operatorname{popcnt}(i \cap k)} A_j \times B_k \\ & = \sum_{j,k} (-1)^{\operatorname{popcnt}(i \cap (j \oplus k))} C_{j \oplus k} \\ & = \sum_{p} (-1)^{\operatorname{popcnt}(i \cap p)} C_p \\ & = \sum_p g(i,p) C_p \\ & = \hat C_i \end{aligned} \]

快速变换

有点抽象,画个 \(n=3\) 的图来模拟一下

快速变换图解

假设现在我们要考虑从 \(n=2\)\(n=3\) 的变换,我们发现左边是 0??,右边是 1??,分别多出一个最高位

设原变换为 \(F_i\) ,目标变换为 \(G_i\)

  • \(\to\)

我们发现左边内部之间的 \(\cap\) 结果不变,因此 \(\forall i<4, G_i \gets F_i\)

  • \(\to\) 右 / 右 \(\to\)

我们发现 0?? \(\cap\) 1?? = 0??,因此其 1 的个数不变,因此 \(\forall i<4,G_i \gets F_{i+4}, G_{i+4} \gets F_i\)

  • \(\to\)

我们发现之前是 ?? \(\cap\) ?? = ??,但是现在在前面加了一个 1,因此 \(-1\) 的指数加一,所以 \(\forall i<4,G_{i+4} \gets -F_{i+4}\)

综上,我们有 \(\forall i<4,G_i=F_i+F_{i+4}, G_{i+4}=F_i-F_{i+4}\)

generally,对于从 \(n-1\)\(n\) 的变换,我们有

\[ \forall i<2^{n-1}, G_i=F_i+F_{i+2^{n-1}}, G_{i+2^{n-1}}=F_i-F_{i+2^{n-1}} \]

时间复杂度 \(O(n2^n)\)

逆变换

对照上面式子易得 (留给读者思考)

模板题

核心代码如下

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
void fwtOr(ll *a,int n,int type) {
    for(int i=1; i<(1<<n); i<<=1)
        for(int j=0; j<(1<<n); j+=(i<<1))
            for(int k=0; k<i; ++k)
                (a[j+k+i]+=type*a[j+k])%=P;
}
void fwtAnd(ll *a,int n,int type) {
    for(int i=1; i<(1<<n); i<<=1)
        for(int j=0; j<(1<<n); j+=(i<<1))
            for(int k=0; k<i; ++k)
                (a[j+k]+=type*a[j+k+i])%=P;
}
void fwtXor(ll *a,int n,int type) {
    for(int i=1; i<(1<<n); i<<=1)
        for(int j=0; j<(1<<n); j+=(i<<1))
            for(int k=0; k<i; ++k) {
                ll u=a[j+k],v=a[j+k+i];
                a[j+k]=(u+v)%P;
                a[j+k+i]=(u-v)%P;
                if(type==-1) {
                    (a[j+k]*=Pi2)%=P;
                    (a[j+k+i]*=Pi2)%=P;
                }
            }
}

进阶理解

\(\rm FWT\) 的变换本质实际上是构造 \(g\) 的过程。

\(\hat A_i = \sum_j g(i,j) A_j\),那么这个 \(n \times n\)\(\lceil\)转移矩阵\(\rfloor\) 被称为位矩阵。

对于不同的卷积形式,我们需要构造不同的位矩阵。

特别注意由于我们需要逆变换,因此矩阵必须有逆。

形式化地,对于 \(\odot\) 卷积,即 \(C_i = \sum_{j \odot k = i} A_j B_k\),我们需要构造如下位矩阵:

\[ \begin{aligned} & \hat C_i = \hat A_i \hat B_i \\ \iff & \sum_j g(i,j) C_j = \sum_j g(i,j) A_j \sum_k g(i,k) B_k \\ \iff & \sum_j g(i,j) C_j = \sum_j \sum_k g(i,j) g(i,k) A_j B_k \\ \iff & g(i,j \odot k) = g(i,j) g(i,k) \end{aligned} \]

写的不是很严谨,请读者自己手推一下。

对于上文的三类卷积,其位矩阵的构造显然为

\[ g_{\cup}(i,j) = [i \cup j = i] \\ g_{\cap}(i,j) = [i \cap j = i] \\ g_{\oplus}(i,j) = (-1)^{\mathrm{popcnt} (i \cap j)} \]

\(\rm FWT\) 的线性性

\(\mathrm{FWT}(A)_i = \hat A_i = \sum_j g(i,j) A_j\),则有如下引理。

\(\mathrm{FWT}(A+B) = \mathrm{FWT}(A) + \mathrm{FWT}(B)\)

根据公式显然。

\(\mathrm{FWT}(cA) = c\mathrm{FWT}(A)\)

这个也是显然(

\(\mathrm{FWT}(A^{-1})_i = g(i,0) \mathrm{FWT}(A)_i^{-1}\)
\[ \begin{aligned} & \mathrm{FWT}(A^{-1})_i \mathrm{FWT}(A)_i \\ & = \sum_j g(i,j) A^{-1}_j \sum_k g(i,k) A_k \\ & = \sum_j \sum_k g(i,j \odot k) A^{-1}_j A_k \\ & = \sum_{p} g(i,p) \sum_{j \odot k = p} A^{-1}_j A_k \\ & = \sum_{p} g(i,p) [p=0] \\ & = g(i,0) \\ \end{aligned} \]

FMT(Fast Möbius Transform) 学习笔记

小 Tips:在计算机语言中 \(\cap\) = & / and\(\cup\) = | / or

定义

定义长度为 \(2^n\) 的序列的 and 卷积 \(A = B * C\)\(A_i=\sum_{j \cap k = i}{B_j \times C_k}\)

考虑快速计算

Zeta 变换

定义长度为 \(2^n\) 的序列的 Zeta 变换

\[ \hat A_i = \sum_{j \subseteq i}{A_j} \]

即子集和

它具有一下性质:

\[ \hat B_i \times \hat C_i = \sum_{j \subseteq i} B_j \times \sum_{k \subseteq i} C_k = \sum_{p \subseteq i} \sum_{j \cap k = p} B_j \times C_k = \sum_{p \subseteq i} A_p = \hat A_i \]

相当于 FFT 中的点值表示法,用以加速卷积过程

快速变换

暴力求解 \(\hat A\) 最优时间复杂度为 \(3^n\),考虑加速

考虑 子集DP ,有一篇很好的博客1,这里大概讲解一下。

其实 \(\hat A\) 本质上是一个高维(\(n\) 维)前缀和

如果我们要求二维前缀和,显然可以通过一下代码实现:

1
2
3
4
5
6
7
8
// 1st
for(int i=1; i<=n; ++i)
    for(int j=1; j<=m; ++j)
            a[i][j]+=a[i-1][j];
// 2nd
for(int i=1; i<=n; ++i)
    for(int j=1; j<=n; ++j)
            a[i][j]+=a[i][j-1];

我们发现 1st 完成后 a[i][j] 表示的是 \(\sum_{k \le i} a_{k,j}\),即 \((i,j)\) 上方的元素和

那么 2nd 便是将 \((i,j)\) 左边操作后的 \(\sum_{k \le j} a_{i,k}\) 累加到 a[i][j] 中,即完成二维前缀和


下面拓展到三维前缀和

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
// 1st
for(int i=1; i<=n; ++i)
    for(int j=1; j<=n; ++j)
        for(int k=1; k<=n; ++k)
            a[i][j][k]+=a[i-1][j][k];
// 2nd
for(int i=1; i<=n; ++i)
    for(int j=1; j<=n; ++j)
        for(int k=1; k<=n; ++k)
            a[i][j][k]+=a[i][j-1][k];
// 3rd
for(int i=1; i<=n; ++i)
    for(int j=1; j<=n; ++j)
        for(int k=1; k<=n; ++k)
            a[i][j][k]+=a[i][j][k-1];

类似地,1st,2nd 中处理除了第 k 平面中的二维前缀和,3rd 在立体空间中把他们加起来,成为三维前缀和


回到 Zeta 变换

\[ \hat A_i = \sum_{j \subseteq i}{A_j} \]

如果我们用二进制 \((x_{n-1}x_{n-2}x_{n-3}\ ...\ x_2x_1x_0)_2\) 表示集合 \(i\),二进制 \((y_{n-1}y_{n-2}y_{n-3}\ ...\ y_2y_1y_0)_2\) 表示集合 \(j\)

那么 \(\hat A_i\) 可以被视为 \(n\) 维前缀和,即 $$ \hat A_i = \sum_{y_{n-1}\le x_{n-1},y_{n-2}\le x_{n-2}, ..., y_1\le x_1,y_0\le x_0} A_j $$

当然,下标都为 0/1

因此,我们有如下 Code:

1
2
3
4
for(int i=0; i<n; ++i)
    for(int j=0; j<(1<<n); ++j)
        if(j&(1<<i))
            a[j]+=a[j^(1<<i)];

以上 Code 的 naive 形式为:

1
2
3
4
for(int i1=0; i1<2; ++i1) for(int i2=0; i2<2; ++i2) ... for(int in=0; in<2; ++in) if(i1>0) a[i1][i2][...][in]+=a[i1-1][i2][...][in];
for(int i1=0; i1<2; ++i1) for(int i2=0; i2<2; ++i2) ... for(int in=0; in<2; ++in) if(i2>0) a[i1][i2][...][in]+=a[i1][i2-2][...][in];
...
for(int i1=0; i1<2; ++i1) for(int i2=0; i2<2; ++i2) ... for(int in=0; in<2; ++in) if(in>0) a[i1][i2][...][in]+=a[i1][i2][...][in-1];

不要质疑我代码编译不通过

因此,我们用 \(O(n2^n)\) 的时间复杂度解决了 Zeta 变换

逆变换

warning:此处不是 莫比乌斯反演

众所周知,前缀和的逆运算即为差分

我们发现,只需要先将最后一维差分,即可将序列处理为 \(n-1\) 维前缀和

1
2
3
4
for(int i1=0; i1<2; ++i1) for(int i2=0; i2<2; ++i2) ... for(int in=0; in<2; ++in) if(in>0) a[i1][i2][...][in]-=a[i1][i2][...][in-1];
...
for(int i1=0; i1<2; ++i1) for(int i2=0; i2<2; ++i2) ... for(int in=0; in<2; ++in) if(i2>0) a[i1][i2][...][in]-=a[i1][i2-2][...][in];
for(int i1=0; i1<2; ++i1) for(int i2=0; i2<2; ++i2) ... for(int in=0; in<2; ++in) if(i1>0) a[i1][i2][...][in]-=a[i1-1][i2][...][in];

但实际上因为前缀和都是无序的,因此我们直接正着做就可以啦

只需要把正变换的 += 改为 -= 就可以了

扩展

有时候,题目给的不是 \(\cap\) ,而是 \(\cup\) 怎么办?

那我们重定义 Zeta 变换 为: $$ \hat A_i = \sum_{i \supseteq j} A_j $$ 再推一下式子: $$ \hat B_i \times \hat C_i =\sum_{i\supseteq j} B_j \times \sum_{i\supseteq k} C_k =\sum_{i\supseteq j,i\supseteq k} B_j \times C_k =\sum_{i\supseteq p} \sum_{j \cup k=p} B_j \times C_k =\sum_{i\supseteq p} A_p =\hat A_i $$ 变换即超集和,高维后缀和

求法也很简单,只用将源代码中的 if(j&(1<<i)) 该为 if(!(j&(1<<i))) 即可


\(\oplus\) 怎么办?

快去学 FWT

例题

再送你一个并/交集的小口诀:下并或,上交与