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;
}