POJ 2154 Color (polya 欧拉函数)

Description

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和p，代表n个珠子，n种颜色，要求不同的项链数，并对结果%p处理。

由Polya定理可得到最后结果等于1/N*∑N^gcd(n,i)
可是，N≤10^9，枚举i 明显会超时，但gcd（n,i）最后得到的结果很少，最多1000多个，于是反过来枚举gcd（i，n）的值L，这里L是n的某个约数，那么我们需要找到0~n-1中有多少个数与n的约数是L，由扩展欧几里得可以知道，必然存在x和y 使得 i*x+n*y=L,由于L是 i 和n最大公约数，所以可以变成
(i/L)*x+（n/L）*y=1，同时mod（n/L），（i/L）*x≡1(mod n/L),即，要找与n/L互质的i/L有多少个，就变成欧拉函数了！
于是，最后答案就变成了∑φ(n/L)*N^(L-1)%p
#include <iostream>
#include <string.h>
using namespace std;
int modular(int a, int b, int p)  //a^b mod p
{
    int ret = 1;
    for (; b; b >>= 1, a = (int) (a)*a%p)
        if (b & 1)
            ret = (int) (ret)*a%p;
    return ret;
}


int eular(int n)
{
    int ret = 1, i;
    for (i = 2; i*i <= n; i++)
        if (n%i == 0)
        {
            n /= i, ret *= i - 1;
            while (n%i == 0)
                n /= i, ret *= i;
        }
    if (n > 1)
        ret *= n - 1;
    return ret;
}
int main()
{
    int T;
    cin>>T;
    while(T--)
    {
        int n,p,ans=0;
        cin>>n>>p;
        for(int i=1;i*i<=n;i++) //枚举长度
        {
            if(n%i==0)//有长度为i的循环，就会有长度为n/i的循环；
            {
                ans=(ans+eular(n/i)%p*modular(n%p,i-1,p))%p;
                if(i!=n/i)//枚举循环长度i，找出相应的i的个数：gcd(n,i)=n/i;
                   ans=(ans+eular(i)%p*modular(n%p,n/i-1,p))%p;
            }
        }
        cout<<ans<<endl;
    }
}