本文移植自个人的原博客（原文）

题目链接

COGS 1608

题目描述

给你一个$2^a∗2^b$的矩阵，在内存中的存放方式是先存第一行的，再存第二行的…每行也是从左到右存放。现在你想把它变成它的转置矩阵（也是一样的储存方式），但是只能用交换操作（即交换两个储存单元的内容），至少需要交换多少步？

输入输出格式

输入格式

第$1$行数数据组数$c(1<=c<=400000)$

接下来有$c$行，每行一组测试数据：两个整数$a,b(0<=a+b<=1000000)$

输出格式

对每组数据输出一个整数，即转置矩阵需要的最少交换次数。因为这个次数可能很大，你只需要输出它模$1000003$的值（没错，这是个素数）。

样例

输入样例

输出样例

1
6
3744

题解

该题所求的是最小步数不是计数，但仍可以通过polya定理求解。将矩阵表示成在内存中的存储方式，并将转置考虑为一种置换。

以$2^1∗2^2$为例，则置换为

到

0 2 4 6
1 3 5 7

表示为

/ 0 1 2 3 4 5 6 7 \
\ 0 2 4 6 1 3 5 7 /

亦或

(0)(1 2 4)(3 6 5)(7)

考虑在任意一个循环节中，需要交换的次数即为（循环节的长度 - 1 ）次，且已知所有循环节长度为矩阵中元素的总个数（即为$2^{a+b}$个），则有

$$Ans=2^{a+b}−循环节个数$$

现在的问题就转换成了求解循环节的个数，已知$(1<=c<=400000)$且$(0<=a+b<=1000000)$，则循环节的个数显然不能通过暴力求解。考虑这种置换对应矩阵中的转置，有一定的规律可循，不妨将元素转换为二进制进行观察。

则置换为

不难发现该转置操作中所有元素的二进制数皆循环右移了$a$位（也可以说循环左移了$b$位），我们也可以稍微推算一下，对于矩阵元素$(i,j)$来说，转置使其变为了$(j,i)$，将其表示为队列，则从第 $(i-1) * 2^a+j$位转置为第$(j-1) * 2^b+i$位，即对应二进制的循环右移。那么可以将问题转化为一个polya定理的经典问题，即项链染色问题。将二进制的表示看为一种染色方案，长度为$a+b$，颜色为$0,1$ 两种，循环右移左移即对应项链的旋转。

以a=2，b=4为例。

则项链为

      a1
     /  \
    a6   a2
    |    |
    a5   a3
     \  /
      a4

但该题的置换群$G= \lbrace 不动，循环右移a位 \rbrace $，显然与原问题不同，那么不妨将项链压缩一下，根据前面的推理，可知循环右移$a$位与左移$b$位是一种操作，那么不妨将$gcd(a,b)$个点压为一个点。则例图转化为

     a1,2
     /  \
 a5,6 -- a3,4
此时对于任意一个点有2^gcd(a,b)中染色方案。

则

$$Ans=2^{a+b}−(\sum_{d|n}k^d \varphi(n/d))/n$$

原式中

$$n=(a+b)/gcd(a,b),k=2^{gcd(a,b)}$$

由于原式需要取模，则需要求$n$的逆元，已知$a+b<=1e6$，模数为$1e6+3$（为质数），易知两者互质，则可通过费马小定理求逆。

（当时没系统地接触过群论，还不太理解Polya定理的含义）

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
const int N = 1e6+1;
const int mod = 1000003;
typedef long long ll;
int a,b;
bool vis[N];
int p[N],phi[N],tot,pow2[N];
int gcd(int a,int b){return b?gcd(b,a%b):a;}
void get(int &x)
{
    x=0;char ch=getchar();
    while(ch<'0'||ch>'9')ch=getchar();
    while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
}
void put(int x)
{  
    int num = 0; char c[15];
    while(x) c[++num] = (x%10)+48, x /= 10;
    while(num) putchar(c[num--]);
    putchar('\n'); 
}
int qpow(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=((ll)ans*a)%mod;
        b>>=1;
        a=((ll)a*a)%mod;
    }
    return ans;
}
int getprime()
{
    phi[1]=1;
    for(int i=2;i<=N-1;++i)
    {
        if(!vis[i])p[++tot]=i,phi[i]=(i-1)%mod;
        for(int t=1;t<=tot&&i*p[t]<=N-1;++t)
        {
            vis[i*p[t]]=1;
            if(i%p[t]==0){phi[i*p[t]]=((ll)phi[i]*p[t])%mod;break;}
            phi[i*p[t]]=((ll)phi[i]*(p[t]-1))%mod;
        }
    }
    pow2[0]=1;
    for(int i=1;i<=N-1;++i)pow2[i]=(pow2[i-1]<<1)%mod;
}

int main()
{
//    freopen("transp2.in","r",stdin);
//    freopen("transp2.out","w",stdout);
    int T;
    getprime();
    scanf("%d",&T);
    while(T--)
    {
        get(a);get(b);
        int lin=0;
        int gcdab=gcd(a,b);
        int n=(a+b)/gcdab;
        int inv;
        inv=qpow(n,mod-2);
        for(int d=1;d<=n;++d)
        {
            if(n%d==0)
            {
                (lin+=(pow2[gcdab*d]*phi[n/d]))%=mod;
            }
        }
        put(((ll)pow2[a+b]-lin*inv%mod+mod)%mod);
    }
}