Transposing is Fun——Polya定理

题面

　　SPOJ419 ( 洛谷 )

解析

　　$SPOJ422$是同一道题，只是$422$需要一些优化，思路是一样的。

　　一个点$(i, j)$转置后变为$(j,i)$，其中$0\leqslant i < 2^a, 0 \leqslant j < 2^b$，可以视作$i*2^b+j$变为了$j*2^a+i$。那么答案就等于$2^{a+b}-$循环节个数。

　　现在的问题就变成如何求循环节个数。将$i*2^b+j$视为首位相接的二进制串，$j*2^a+i$是其右移$b$位所得，问题再次转化：一个含有$a+b$个点的环，每个点有$2$种颜色，若一种染色方案向右旋转$b$位与另一方案相同，则将两种方案视为同一种方案，问有多少种不同的染色方案。

　　长为$a+b$的环上一个点不断向右移$b$位，回到起点所需要的步数为$\frac{lcm(a+b, b)}{b}$，因此该环上有$\frac{a+b}{lcm(a+b,b)/b}=gcd(a+b,b)=gcd(a,b)$条不相交循环，令$g=gcd(a,b)$，将相邻$g$位看作一个点，于是问题又一次转化：一个含有$\frac{a+b}{g}$个点的环，每个点有$2^g$种颜色，若一种染色方案向右旋转$1$位与另一方案相同，则将两种方案视为同一种方案，问有多少种不同的染色方案。

　　这是$polya$计数的经典题了，但之前没有总结过，这里就稍微写一下，令$n=\frac{a+b}{g}$，$m=2^g$

　　则有：$$\begin{align*}Ans&=\frac{1}{n}\sum_{i=1}^n m^{gcd(i,n)}\\&=\frac{1}{n}\sum_{d|n}m^d\sum_{i=1}^n[gcd(i,n)==d]\\&=\frac{1}{n}\sum_{d|n}m^d\sum_{i=1}^{\frac{n}{d}}[gcd(\frac{n}{d},i)==1]\\&=\frac{1}{n}\sum_{d|n}m^d\varphi(\frac{n}{d})\end{align*}$$

　　还有一些实现上的优化，这个就看代码了。

　代码：

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<vector>
using namespace std;
typedef long long ll;
const int maxn = 1000005, mod = 1000003;

int add(int x, int y)
{
    return x + y < mod? x + y: x + y - mod;
}

int rdc(int x, int y)
{
    return x - y < 0? x - y + mod: x - y;
}

ll qpow(ll x, int y)
{
    ll ret = 1;
    while(y)
    {
        if(y&1)
            ret = ret * x % mod;
        x = x * x % mod;
        y >>= 1;
    }
    return ret;
}

int cnt, pri[maxn], phi[maxn<<1], fct[maxn<<1];
bool notp[maxn<<1];

void Euler()
{
    phi[1] = 1;
    for(int i = 2; i <= 2000000; ++i)
    {
        if(!notp[i])
        {
            pri[++cnt] = i;
            phi[i] = ((i - 1 < mod)? i - 1: i - 1 - mod);
            fct[i] = i;
        }
        for(int j = 1; j <= cnt; ++j)
        {
            if(i * pri[j] > 2000000)    break;
            notp[i*pri[j]] = 1;
            fct[i*pri[j]] = pri[j];
            if(i % pri[j] == 0)
            {
                phi[i*pri[j]] = phi[i] * pri[j] % mod;
                break;
            }
            phi[i*pri[j]] = phi[i] * (pri[j] - 1) % mod;
        }
    }
}

int T, n, m, g, d, ans;
int pw2[maxn<<1];

int gcd(int x, int y)
{
    return (y == 0)? x: gcd(y, x % y);
}

int a[50], num[50], tot;

void dfs(int x, int idx)
{
    if(x == tot + 1)
    {
        ans = add(ans, qpow(pw2[g], idx) * phi[d/idx] % mod);
        return ;
    }
    for(int i = 0; i <= num[x]; ++i)
    {
        dfs(x + 1, idx);
        idx *= a[x];
    }
}

int main()
{
    Euler();
    pw2[0] = 1;
    for(int i = 1; i <= 2000000; ++i)
        pw2[i] = add(pw2[i-1], pw2[i-1]);
    scanf("%d", &T);
    while(T --)
    {
        scanf("%d%d", &n, &m);
        if(!n || !m)
        {
            printf("0\n");
            continue;
        }
        g = gcd(n, m);
        d = (n + m) / g;
        ans = 0;
        tot = 0;
        for(int i = d; i != 1; i /= fct[i])
        {
            if(fct[i] == a[tot])
            {
                ++ num[tot];
            }
            else
            {
                a[++tot] = fct[i];
                num[tot] = 1;
            }
        }
        dfs(1, 1);
        printf("%d\n", rdc(pw2[n+m], ans * qpow(d, mod - 2) % mod));
    }
    return 0;
}