SP422 TRANSP2 - Transposing is Even More Fun 题解

题目链接：TRANSP2 - Transposing is Even More Fun

题意：

假设你有一个 \(2^a\times 2^b\) 的数组。它在内存中以通常的方式按顺序存储，先是第一行的值，然后是第二行的值，依此类推。你希望对它进行转置，但没有任何额外的内存。你唯一能执行的操作是交换两个内存单元的内容。进行转置所需的最少操作数是多少？

输入格式：

输入的第一行包含测试用例的数量 \(c~(1 \le c \le 4\times 10^6)\) 。

每个测试用例由两个整数 \(a\) 和 \(b\) 组成 \((0 \le a + b \le 1\times 10^6)\) 。

输出格式：

对于每个测试用例，输出转置一个 \(2a\times 2b\) 数组所需的最少交换次数。

由于这个数可能非常大，你需要输出它除以 \(1000003\) 的余数（是的，\(1000003\) 是一个质数:) ）。

我们可以把答案看作将一个序列转化为排列 \(p=0,1,\cdots\) 的最少交换数。

意思对于每个数 \(i\) 我们连边 \(i \to p_i\) 会形成若干个循环置换，每个大小为 \(x\) 的循环置换贡献 \(x-1\) 的次数

现在的问题是如何找到这个序列中循环置换的数量。

注意到对于矩阵元素 \(a_{ij}\) ，它在序列上的位置为 \(2^bi+j\) ，而在转置矩阵 \(a^{\rm T}\) 上的位置为 \(2^aj + i\) 。

如果我们把这两个数看作一个 \(a+b\) 位的二进制数，并首尾相接，那么交换操作就等价于旋转 \(b\) 个单位。

如果 \(x\) 在旋转后变为 \(y\) ，我们就连边 \(x\to y\) ，显然最终形成的环的个数就是我们要求的循环置换的个数。

注意到旋转 \(b\) 个单位会把 \(a+b\) 个二进制位划分成 \(\gcd(a+b,\,b)\) 个大小为 \(\frac{a+b}{\gcd(a + b,\,b)}\) 的环

因为 \(\gcd(a + b,\, b)=\gcd(a, b)\) ，所以也就是 \(\gcd(a,b)\) 个大小为 \(\frac{a+b}{\gcd(a,b)}\) 的等价类。

不妨记 \(g = \gcd(a,b),~n = \frac{a+b}{\gcd(a,b)}\) ，根据 Burnside 引理 \[ |X / G|=\frac{1}{|G|} \sum_{\varphi_i \in G} \left|X^{\varphi_i}\right| \] 其中 \(X^{\varphi_i}\) 为 \(\varphi_i\) 作用下的不动点集合，本题中置换 \(\varphi_i\) 表示将 \(g\) 个大小为 \(n\) 的环旋转 \(i\) 个单位。

而 \(\varphi_i\) 作用下的不动点数量为 \(\left(2^{\operatorname{gcd}(n, i)}\right)^g=2^{g \operatorname{gcd}(n, i)}\) ，这里 \(2\) 是每个等价类可以染成 \(0\) 或 \(1\) 。

那么答案就是 \[ \begin{aligned} \mathrm{Ans} &= \sum_{i=1}^n 2^{g \operatorname{gcd}(n, i)} \\[3pt] &= \sum_{d \mid n} 2^{g d} \sum_{i=1}^{\frac{n}{d}}\left[\operatorname{gcd}\left(i, \frac{n}{d}\right)=1\right] \\ \\[3pt]&= \sum_{d \mid n} 2^{g d} \varphi\left(\frac{n}{d}\right) \end{aligned} \] 提示：这里 \(\varphi(n)\) 是欧拉函数，中括号是艾弗森括号。

那么我们只需要预处理 \(\varphi\) 以及所有数的约数即可。

这里学到了一种妙妙办法，可以枚举每个数的倍数并和它连边，这样每次只需要遍历邻接表就可以了。

时间复杂度 \[ \mathcal{O}\left(V \log V + T \max_{1\le n \le V}\{ d(n)\}\right) \] 其中 \(V = a + b\) 的值域，\(d(n)\) 为 \(n\) 的因子个数。根据表格可知 \(V \le 10^6\) 时 \(\max d(n) \le 240\) 。

代码：

#include <bits/stdc++.h>
using namespace std;
// #define int long long
// #define INF 0x3f3f3f3f3f3f3f3f
typedef long long ll;
const int mod = 1000003;
void up(int &x, int y) { x < y ? x = y : 0; }
void down(int &x, int y) { x > y ? x = y : 0; }
void add(int &x, int y) { (x += y) >= mod ? x -= mod : 0; }
int dec(int x, int y) { return (x - y < 0) ? x - y + mod : x - y; }
#define N ((int)(1e6 + 15))

char vis[N];
struct Edge { int val, next; } e[N * 14];
int pos = 1, pcnt, pw2[N], head[N], phi[N], p[N];
int gcd(int a, int b) { return !b ? a : gcd(b, a % b); }
int qpow(int a, int b)
{
    int r = 1;
    for(a %= mod; b; b >>= 1, a = 1ll * a * a % mod)
        if(b & 1) r = 1ll * r * a % mod;
    return r;
}
void addEdge(int u, int v) { e[++pos] = {v, head[u]}; head[u] = pos; }
void init(int n)
{
    pw2[0] = 1;
    for(int i = 1; i <= n; i++)
    {
        pw2[i] = 2ll * pw2[i - 1] % mod;
        for(int j = i; j <= n; j += i) addEdge(j, i);
    }
    phi[1] = 1;
    for(int i = 2; i <= n; i++)
    {
        if(!vis[i]) { p[++pcnt] = i; phi[i] = i - 1; }
        for(int j = 1; j <= pcnt && i * p[j] <= n; j++)
        {
            const int x = i * p[j]; vis[x] = true;
            if(i % p[j]) { phi[x] = phi[i] * phi[p[j]]; }
            else { phi[x] = phi[i] * p[j]; break; }
        }
    }
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    init(N - 5); int a, b, qwq; cin >> qwq;
    while(qwq--)
    {
        cin >> a >> b; if(!a || !b) { cout << "0\n"; continue; }
        int g = gcd(a, b), n = (a + b) / g, res = 0;
        for(int i = head[n], d; i; i = e[i].next)
        {
            d = e[i].val;
            add(res, 1ll * pw2[d * g] * phi[n / d] % mod);
        }
        cout << dec(pw2[a + b], 1ll * res * g % mod * qpow(a + b, mod - 2) % mod) << '\n';
    }
    return 0;
}