学步园

Ploya定理+数论知识...

置换群是由旋转置换构成的，因此|G| = n

基本知识：循环的个数为gcd（i,n），循环的长度为l = n/gcd(i,n)      ans = sigma(n^(gcd(n,i)-1))  1<=i<=n

由于n很大所以需要利用欧拉函数进行优化。。。

枚举的方法有两种：

<1>将n进行质因子分解，然后根据质因子分解枚举它的约数ys，这样答案便为sigma（euler(n/ys)*n^(ys-1)）

<2>枚举循环的长度l,答案为 sigma(euler(l)*n^(n/l-1))

我采用的是第一种枚举的方法，比第二种方法稍微繁琐一些，不过都很简单啦。。。

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long ll;
const int MAXN = 60000;
bool is_prime[MAXN];
ll prime[MAXN];
ll facs[110], ccount[110];
ll ans, N, P;
int len, cnt, x;

void Init()
{
    int i, j;
    len = 0;
    memset(is_prime, true, sizeof(is_prime));
    is_prime[0] = is_prime[1] = false; prime[len++] = 2;
    for(i = 4; i < MAXN; i += 2)
        is_prime[i] = false;
    for(i = 3; i * i <= MAXN; i += 2)
    {
        if(is_prime[i])
        {
            prime[len++] = i;
            for(j = i * i; j < MAXN; j += i)
                is_prime[j] = false;
        }
    }
    for( ; i < MAXN; ++i)
        if(is_prime[i])
            prime[len++] = i;
}


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

ll mul_mod(ll a, ll b, ll c)
{
    ll ret = 0;
    while(b)
    {
        if(b & 1)
            ret = (ret + a) % c;
        a = a * 2 % c;
        b >>= 1;
    }
    return ret;
}

ll power_mod(ll a, ll b, ll c)
{
    ll ret = 1;
    while(b)
    {
        if(b & 1)
            ret = mul_mod(ret, a, c);
        a = mul_mod(a, a, c);
        b >>= 1;
    }
    return ret;
}

void cal(ll n)
{
    cnt = 0;
    memset(ccount, 0, sizeof(ccount));

    for(int i = 0; i < len && prime[i]*prime[i] <= n; ++i)
    {
        if(n % prime[i] == 0)
        {
            facs[cnt] = prime[i];
            while(n % prime[i] == 0)
            {
                ccount[cnt]++;
                n /= prime[i];
            }
            cnt++;
        }
    }
    if(n > 1)
        facs[cnt] = n, ccount[cnt++] = 1;
    return ;
}

ll euler(ll n)
{
    ll ret = n;
    for(int i = 0; i < len && prime[i]*prime[i] <= n; ++i)
    {
        if(n % prime[i] == 0)
        {
            ret = ret/prime[i]*(prime[i]-1);
            while(n % prime[i] == 0)
            {
                n /= prime[i];
            }
        }
    }
    if(n > 1)
        ret = ret/n*(n-1);
    return ret;
}

ll power(ll a, ll b)
{
    ll ret = 1;
    while(b)
    {
        if(b & 1)
            ret = ret * a;
        a = a * a;
        b >>= 1;
    }
    return ret;
}

void dfs(ll cur, int id)
{
    if(id == cnt)
    {
        ans = (ans + euler(N/cur)*power_mod(N, cur-1, P)) % P;
        return ;
    }
    for(int i = 0; i <= ccount[id]; ++i)
    {
        dfs(cur*power(facs[id], i), id + 1);
    }
    return ;
}

int main()
{
    //freopen("aa.in", "r", stdin);
    //freopen("bb.out", "w", stdout);

    Init();
    cin >> x;
    while(x--)
    {
        cin >> N >> P;
        if(N == 1)
        {
            cout << N % P << endl;
            continue;
        }
        ans = 0;
        cal(N);
        dfs(1, 0);
        cout << ans << endl;
    }
    return 0;
}