SPOJ 422 Transposing is Even More Fun（Polya计数）

TRANSP2 - Transposing is Even More Fun   

Problem

Suppose you are given a 2ax2b array. It is stored sequentially in memory in the usual way, first values in the first row, then values in the second one and so on. You would like to transpose it, but you don't have any additional memory. The only operation that you can perform is swapping contents of two memory cells. What is minimal number of such operations required for transpostion?

Input

First line of input contains number of test cases c (1<=c<=400000). Each test case consists of two integers a,b (0<=a+b<=1000000).

Output

For each test case output minimal number of swaps required to transpose an 2ax2b array. As it can be quite large, you have to output its remainder when divided by 1000003 (yes, it's a prime number :).

Example

Input:31 12 25 7Output:163744

        大致题意：给你一个2ax2b的矩阵，然后让你把它转置为2bx2a的矩阵，每次可以交换两个数字，问最少需要多少步交换。

        首先，我还是借鉴一下别人的图做一下解释。我们就以a=1，b=2为例子解释：

             

                    

        意思是，我把某一个位置的坐标转换成二进制表示。例如原本3的位置是010，转置之后就是100。通过观察，我们发现一个规律，对于一个一开始的位置的地址表示，把他进行旋转之后，对应地址的二进制形式相当于移动了b位。正如图中010->100,100->001……这样的话，就能够写出下面那张图的下半部分的循环，只要按照循环的移动，我们就可以将矩阵进行转置。

        考虑非最优情况下的交换次数。每次交换我都确定一个位置，那么总共就需要2ax2b次交换。再考虑一个循环中的元素，我们发现，对于一个长度为k的循环，只需要k-1次交换就可以把每个元素就位。那么，我们的最优解就是非最优的次数减去这样的循环个数。图中有4个循环，那么结果就是2^(1+2)-4=4次。然后仔细观察循环，可以发现一个循环内的地址二进制表示在移动b位的情况下是本质相同的，那么求循环的个数，相当于求长度a+b的01串中，有多少个移动b位后本质不同的，问题就转换成了一个Polya计数的问题了。

        但是，这里又不是像常规Polya一样的。普通的Polya，旋转置换有多个，这里却只能有移动b位这一个置换，故burnside引理在这里不能直接对移动b位使用。既然如此，我们就得转换一下了，考虑正常的时候，移动b位的循环节个数为gcd(b,a+b)，然后循环长度为(a+b)/gcd。再次说一下，这里不能直接用burnside引理计算不动点个数。我们要做的，是考虑把每gcd个数字绑在一起，原本移动b位，现在就能够移动任意位。那么问题就转换为了，有长度为(a+b)/gcd的一个项链，每个项链可以用2^gcd(每一捆里面有gcd个数字，每个数字取0、1)种颜色来染色，问在旋转置换下本质不同的染色方案数。此时再按照套路直接用burnside引理计算即可。

        想了想还是再解释一下吧，防止自己以后老年痴呆……解释一下，为什么要把每gcd个数字绑在一起。我们假设a=2、b=4。那么gcd(b,a+b)=2。每次移动4位本质相同。不能直接用burnside引理原因是，这里有a+b种置换，但是有效的却只有一种。然后观察移动b位的置换，在这种置换下第0、4、2位是在一个循环内，原本是0->4->2，但是我们完全可以看成0->2->4（每次都转两个b位），如此如果把每两个数字捆绑，捆绑后的堆记为a、b、c，那么循环变成了a->b->c。考虑另一个循环也是如此，说明捆绑之后对计算结果不影响。这样，我们就相当于没有了旋转b位的限制，对于a、b、c移动1、2、3位都是有效置换，于是就可以用burnside引理解决。

        最后要注意，本题卡常数，所以在用欧拉函数优化的同时，要提前预处理出欧拉函数和逆元，减少计算时间。还有，当a或b为0时要特殊处理，不然会RE……具体见代码：

#include<bits/stdc++.h>#define LL long long#define mod 1000003using namespace std;int a,b,prime[mod],phi[mod],tot;LL Pow[mod],inv[mod];bool v[mod];void init(){    tot=0;    memset(v,0,sizeof(v));    Pow[0]=inv[1]=phi[1]=1;    for(int i=1;i<mod;i++)        Pow[i]=Pow[i-1]*2LL%mod;for(int i=2;i<mod;i++)    {        inv[i]=(mod-mod/i)*inv[mod%i]%mod;//线性求逆元if(!v[i]) prime[tot++]=i,phi[i]=i-1;for(int j=0;j<tot&&i*prime[j]<mod;j++){v[i*prime[j]]=1;if(i%prime[j]==0){phi[i*prime[j]]=phi[i]*prime[j];//线性求phi                break;} else phi[i*prime[j]]=phi[i]*(prime[j]-1);}    }}int main(){    init();    int T_T;    cin>>T_T;    while(T_T--)    {        scanf("%d%d",&a,&b);        if (!a||!b){puts("0");continue;}//有0时特判，不然下面除法会RE        int gcd=__gcd(b,a+b);        int len=(a+b)/gcd; LL ans=0;        for(int i=1;i*i<=len;i++)        {            if (len%i) continue;            ans=(ans+Pow[gcd*i]*phi[len/i]%mod)%mod;            if (i*i!=len) ans=(ans+Pow[gcd*len/i]*phi[i]%mod)%mod;        }        ans=ans*inv[len]%mod;        printf("%lld\n",(Pow[a+b]-ans+mod)%mod);    }    return 0;}