1、无法用传统方法求逆元,需借助二分模拟大数乘法取模运算
2、对于旋转1个单位这类置换,在该条件下不变的方案为5种。
#include<cstdio>
#include<algorithm>
#pragma comment(linker, "/STACK:1024000000,1024000000")
using namespace std;
typedef __int64 ll;
ll mod;
const int MAXN = 32000;
bool isPrime[MAXN] = {false};
int prime[MAXN], np = 0;
int nn;
ll bigMult(ll a, ll b)
{
ll res = 0;
if (a < 0) a += mod;
if (b < 0) b += mod;
while (b)
{
if (b & 1)
{
res += a;
if (res >= mod) res -= mod;
}
a <<= 1;
if (a >= mod) a -= mod;
b >>= 1;
}
return res;
}
const ll bf[2][2] = {3,1,-1,0};
const ll bb[2][2] = {3,1,0,0};
struct matrix
{
ll da[2][2];
matrix(){}
void init()
{
memset(da, 0, sizeof da);
da[0][0] = da[1][1] = 1;
}
matrix operator * (const matrix & a)
{
matrix res;
for (int i = 0; i< 2; ++i)
for (int j = 0; j< 2; ++j)
{
res.da[i][j] = 0;
for (int k = 0; k< 2; ++k)
res.da[i][j] = (res.da[i][j] + bigMult(da[i][k], a.da[k][j]))%mod;
}
return res;
}
matrix operator ^ (int a)
{
matrix res, p = *this;
res.init();
while ( a)
{
if (a & 1) res = res * p;
p = p*p;
a >>= 1;
}
return res;
}
};
void initPrime()
{
for (int i = 2; i< MAXN; ++i)
{
if(!isPrime[i]) prime[np++] = i;
for (int j = 0; j<np && i*prime[j]<MAXN; ++j)
{
isPrime[i*prime[j]] = true;
if (i % prime[j] == 0) break;
}
}
}
int fenzi[100][2];
int nfen;
void divide(int a)
{
for (int i = 0; i< np && prime[i]*prime[i] <= a; ++i)
{
if (a % prime[i] == 0)
{
fenzi[nfen][0] = prime[i]; fenzi[nfen][1] = 0;
while (a % prime[i] == 0)
{
a /= prime[i];
fenzi[nfen][1]++;
}
nfen++;
}
}
if (a != 1)
{
fenzi[nfen][0] = a; fenzi[nfen++][1] = 1;
}
}
ll resAll;
ll getCir(int t)
{
ll res = 0;
if (t == 1) return 1;
if (t == 2) return 5;
matrix ff, bg;
memcpy(bg.da, bb, 32); memcpy(ff.da, bf, 32);
bg = bg*(ff^(t-2));
res = (bg.da[0][0] + (bg.da[0][0] - bg.da[0][1] - 1)*2%mod)%mod;
if (res < 0) res += mod;
return res;
}
void dfs(int i, int p1, int p2, int x)
{
if (i == nfen)
{
int t = nn/x;
resAll = (resAll + bigMult(getCir(t), x/p1*p2))%mod;
}else
{
int cnt = 0, p = fenzi[i][1], q = fenzi[i][0];
for ( ;cnt <= p; ++cnt)
{
if (!cnt) dfs(i+1, p1, p2, x);
else dfs(i+1, p1*q, p2*(q-1), x);
x *= q;
}
}
}
ll solve()
{
nfen = 0;
resAll = 0;
divide(nn);
dfs(0, 1, 1, 1);
return resAll/nn;
}
int main()
{
#ifndef ONLINE_JUDGE
freopen("in.txt", "r", stdin);
#endif
initPrime();
while (scanf("%d%I64d", &nn, &mod) != EOF)
{
mod *= nn;
printf("%I64d\n", solve()%(mod/nn));
}
return 0;
}