Polya + Euler优化。由于n太大,所以需要用Euler优化。设L = n / gcd(n, i),t = i / gcd(i, n),则gcd(L, t) = 1.又有i < n所以,t < L,且t和L互质。于是满足gcd(n, i) = n / L的i的个数为Euler(L)。对于题目原本的求和式应该是sum(pow(n, gcd(i, n) - 1))。所以我们枚举n的所有因子L,然后求和式就变成了
sum(Euler(L) * pow(n, n / L - 1))。
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<vector>
#include<cmath>
#include<set>
#include<map>
#define LL long long
#define CLR(a, b) memset(a, b, sizeof(a))
#define REP(i, n) for(int i = 0; i < n; i ++)
using namespace std;
const int N = 100100;
bool isp[N];
vector<int> p;
void get_P()
{
CLR(isp, true);p.clear();
for(int i = 2; i < N; i ++)
{
if(isp[i])
{
p.push_back(i);
if(i < 1111) for(int j = i * i; j < N; j += i)
{
isp[j] = false;
}
}
}
}
int Euler_phi(int n)
{
int ret = n;
for(int i = 0; p[i] * p[i] <= n; i ++) if(n % p[i] == 0)
{
ret = ret / p[i] * (p[i] - 1);
while(n % p[i] == 0) n /= p[i];
}
if(n > 1) ret = ret / n * (n - 1);
return ret;
}
int Pow(int a, int b, int mod)
{
int ret = 1;a %= mod;
while(b)
{
if(b & 1) ret = ret * a % mod;
a = a * a % mod;
b >>= 1;
}
return ret;
}
int main()
{
int n, ans, mod, x;
get_P();
scanf("%d", &x);
while(x --)
{
scanf("%d%d", &n, &mod);ans = 0;
for(int i = 1; i * i <= n; i ++)
{
if(n % i != 0) continue;
ans += Euler_phi(i) % mod * Pow(n, n / i - 1, mod) % mod;
if(n / i != i)
ans += Euler_phi(n / i) % mod * Pow(n, i - 1, mod) % mod;
ans %= mod;
}
printf("%d\n", ans);
}
}