POJ2409_Let it Bead_Polya定理

题意

一串项链上有 n 棵珠子，给每颗珠子图上 m 种颜色中的一种，问一共有多少种不同的涂法。特别的，如果一种涂法通过任意角度的旋转，或沿项链某条对称翻转后，能得到另一种涂法，则认为这两种是同一种。

思路

polya定理

直观地理解

只考虑旋转，不考虑对称的情况

首先计算，旋转 k 个位置后，和原来相同的染色方案的个数。
首先，按照顺时针顺序从 0 到 n - 1 给珠子编号。显然，i 珠子 与 (i + k) mod n 珠子颜色相同。传递下去可知，i 珠子和 (i + k * t) mod n 珠子都相同， 它们实际上就构成了“置换”中的一个循环，这个循环中珠子的集合又叫做 i 的“轨迹。
那么，这一条轨迹上有多少颗珠子呢？实际上要求 k * t = 0 (mod n) 的最小 t。显然 t = n / gcd(k, n)。
只要我们把同一条轨迹上的珠子都染成同一种颜色，这种方案在旋转 k 个位置后就与自身相同了。那么，怎么求轨迹的个数呢？
旋转 k 个位置下，每一个轨迹中石头的个数 t = n / gcd(k, n), 则轨迹的个数就是 n / t = gcd(k, n)。旋转 k 个位置与本身相同的染色方案数就是 (m 的 gcd(k, n) 次方)。
接下来比较暴力的做法就是，枚举 k 的值，把 (m 的 gcd(k, n) 次方) 加起来，最后除以 n。但是，如果 n 比较大的话，这样容易 T。
不过，gcd(k, n) 只有有限个不同的值，可以把相同的值合并在一起计算。
设 d 是 n 的约数，我们要统计满足 gcd(k, n) = d 的 k 的个数。令 k = d * t (0 <= t < n / d)。则， d = gcd(k, n) = gcd(d * t, n) = d * gcd(t, n / d)。即 gcd(t, n / d) = 1, 满足条件的个数就是 phi(n / d)。
到此，计算表达式变为 1 / n * sigma( (m 的 d 次方) * phi(n / d) )。
计算这个式子需要对 n 的所有约数 d 求欧拉函数。d 的质因数也是 n 的质因数，如果之前求出了 n 的质因数，求 d 的欧拉函数时就只需要 logn 时间。

最后再来考虑对称，对称变换需要对 n 进行奇偶分类。
若 n 为奇数，对称轴只能过一个顶点和中心，共有 n 条轴，n 种变换。在每种变换中，轴上的点循环节之和自己循环，循环节长度为 1，其余点和对称点循环，长度为 2。循环节个数为 (n + 1) / 2。
当 n 为偶数时，对称轴有两种，一种是过两个点，有 n / 2 ， n / 2 中变换。在每种变换中，轴上的点都是自己循环，其他点都是和对称点循环，长度为2。 循环节个数为 n / 2 + 1。另一种是过两条边的中点，n / 2 中变换。在这种变换下， 循环节的个数为 n / 2。

至此，本问题解决。需要注意的一点是，加上对称后，不论奇偶，总共有 n * 2 种变换， 即最后应除以 (n * 2)。

代码部分给出上述两种思路的实现。

链接

http://poj.org/problem?id=2409

代码

优化版本

#include<cstdio>
#include<iostream>

using namespace std;

typedef long long LL;

const int maxn = 1000;

LL n, m;
LL prime[9] = {2, 3, 5, 7, 11, 13, 17, 19, 23}, factor[maxn], it;
LL divisor[maxn], it0;

void get_factor(LL n){
    it = 0;
    for(int i = 0; i < 9; ++i){
        if(n % prime[i] == 0){
            factor[it++] = prime[i];
            while(n % prime[i] == 0) n /= prime[i];
        }
    }
    if(n > 1) factor[it++] = n;
}

void get_divisor(LL n){
    it0 = 0;
    for(int i = 1; i <= n; ++i){
        if(n % i == 0){
            divisor[it0++] = i;
        }
    }
}

LL pow(LL a, LL b){
    LL res = 1;
    while(b > 0){
        if(b & 1) res *= a;
        a *= a;
        b >>= 1;
    }
    return res;
}

int main(){
    while(scanf("%lld %lld", &m, &n) == 2){
        if(n == 0 && m == 0) break;

        if(n == 0){
            cout << 0 << endl;
            continue;
        }

        get_factor(n);
        get_divisor(n);

        LL res = 0;
        for(int i = 0; i < it0; ++i){
            LL euler = divisor[i];
            for(int j = 0; j < it; ++j){
                if(divisor[i] % factor[j] == 0) euler = euler / factor[j] * (factor[j] - 1);
            }

            res += euler * pow(m, n / divisor[i]);
        }

        if(n & 1){
            res += n * pow(m, (n + 1) / 2);
        }else{
            res += n / 2 * pow(m, n / 2 + 1) + n / 2 * pow(m, n / 2);
        }
        res = res / (n * 2);

        cout << res << endl;
    }

    return 0;
}

枚举版本

#include<iostream>
#include<cstdio>

using namespace std;

typedef long long LL;

LL n, m;

LL gcd(LL a, LL b){
    if(b == 0) return a;
    return gcd(b, a % b);
}

LL pow(LL a, LL b){
    LL res = 1;
    while(b > 0){
        if(b & 1) res *= a;
        a *= a;
        b >>= 1;
    }
    return res;
}

int main(){
//  freopen("in.txt", "r", stdin);

    while(scanf("%lld %lld", &m, &n) == 2){
        if(n == 0 && m == 0) break;

        if(n == 0){
            cout << 0 << endl;
            continue;
        }

        LL res = 0;
        for(int i = 0; i < n; ++i){
            res += pow(m, gcd(i, n));
        }
        if(n & 1){
            res += n * pow(m, (n + 1) / 2);
        }else{
            res += n / 2 * pow(m, n / 2 + 1) + n / 2 * pow(m, n / 2);
        }
        res = res / (n * 2);

        cout << res << endl;
    }

    return 0;
}