Color POJ - 2154 polya定理 + 欧拉优化

[题目链接](http://poj.org/problem?id=2154)

---

Beads of N colors are connected together into a circular necklace of N beads (N<=1000000000). Your job is to calculate how many different kinds of the necklace can be produced. You should know that the necklace might not use up all the N colors, and the repetitions that are produced by rotation around the center of the circular necklace are all neglected.

You only need to output the answer module a given number P.

Input

The first line of the input is an integer X (X <= 3500) representing the number of test cases. The following X lines each contains two numbers N and P (1 <= N <= 1000000000, 1 <= P <= 30000), representing a test case. 

Output

For each test case, output one line containing the answer. 

Sample Input

5
1 30000
2 30000
3 30000
4 30000
5 30000

Sample Output

1
3
11
70
629

---

题意:

n种不同颜色的珠子各有无数个,问能组成多少种不同的具有 n 个珠子的手链.

思路:

polya 定理

代码:

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;

typedef long long LL;
const int maxn = 33333;
bool isprim[maxn];
int pri[maxn];
int num;

void prime(){
    memset(isprim, false, sizeof(isprim));
    num = 0;
    for(int i = 2; i < maxn; ++i){
        if(!isprim[i]) pri[++num] = i;
        for(int j = 1; j <= num && i * pri[j] < maxn; ++j){
            isprim[i * pri[j]] = true;
            if (i % pri[j] == 0) break;
        }
    }
}

LL euler(LL n){
    LL ans = n;
    for(int i = 1; i <= num && pri[i] * pri[i] <= n; ++i){
        if(n % pri[i] == 0){
            ans -= ans / pri[i];
        while(n % pri[i] == 0) n /= pri[i];
        }
    }if (n > 1) ans -= ans / n;
    return ans;
}

LL qpow(LL x, LL n, LL mod){
    x %= mod;
    LL ans = 1;
    while(n){
        if (n&1) ans = (ans * x) % mod;
        x = (x * x) % mod;
        n >>= 1;
    }return ans;
}

int main(){
    ios::sync_with_stdio(false);
    int t;
    prime();
    cin >> t;
    while(t--){
        int n, p;
        cin >> n >> p;
        LL ans = 0;
        for(int i = 1; i * i <= n; ++i)
            if(n % i == 0){
                if(i * i == n)
                    ans += euler(i) * qpow(n, i-1, p);
                else ans += euler(i) * qpow(n, n/i - 1, p) + euler(n/i) * qpow(n, i - 1, p);
                ans %= p;
            }cout<<ans<<endl;
    }return 0;
}