HDU 2481 Toy

是上一题FJOI轮状病毒的加强版，
前置知识：
$Burnside$引理：所有等价类的数目为其所有置换的不动点数目的平均值。
具体见《算法竞赛入门经典——训练指南》
无旋情况下的递推式是$$F[n]=3\ast F[n-1]-F[n-2]+2$$
然后套$Burnside$引理
环置换下$$ans=\frac{1}{n}\times \sum _{d|n}F(d)\varphi (\frac{n}{d})$$

$F$用矩阵快速幂计算

一个细节：$n$和$m$不一定互质 逆元可能不存在
这里要用$$\frac{a}{b}modp=\frac{amod(b\times p)}{b}$$
证明：设$ans=\frac{a}{b}+k\times p$
$ans\times b=a+b\times k\times p=a(\mathrm{m}\mathrm{o}\mathrm{d}b\times p)$
$ans=\frac{a}{b}\mathrm{m}\mathrm{o}\mathrm{d}b\times p$
所以计算过程中对n*m取模 最后再除n就行
但是$n\ast m$会爆ll
此时需要快速乘（龟速乘

#include<bits/stdc++.h>
using namespace std;
const int maxn = 1e6 + 7;
#define ll long long
ll ans, m, n, tot, isnp[maxn], md, pri[maxn];
ll ksc(ll a, ll b) {
    ll res = 0;
    while (b) {
        if (b & 1) {
            res += a;
            if (res >= md) res -= md;
        }
        a += a;
        if (a >= md) a -= md;
        b >>= 1;
    }
    return res;
}
struct _ {
    ll a[4][4];
    _ operator * (const _ &rhs) const {
        _ res;
        for (int i = 1; i <= 3; i++) {
            for (int j = 1; j <= 3; j++) {
                res.a[i][j] = 0;
                for (int k = 1; k <= 3; k++)
                    res.a[i][j] = (res.a[i][j] + ksc(a[i][k],rhs.a[k][j])) % md;
            }
        }
        return res;
    }
};
_ ksm(_ a, ll b) {
    _ res;
    for (int i = 1; i <= 3; i++)
        for (int j = 1; j <= 3; j++)
            res.a[i][j] = (i == j);
    while (b) {
        if (b & 1) res = res * a;
        a = a * a;
        b >>= 1;
    }
    return res;
}
void getpri() {
    for (int i = 2; i <= 1e6; i++) {
        if (!isnp[i]) pri[++tot] = i;
        for (int j = 1; (ll)i * pri[j] <= 1e6; i++) {
            isnp[i * pri[j]] = 1;
            if (i % pri[j] == 0) break;
        }
    }
}
ll phi(ll x) {
    ll res = x;
    for (int i = 1; i <= tot && pri[i] * pri[i] <= x; i++) {
        if (x % pri[i] == 0) {
            res = res / pri[i] * (pri[i]-1);
            while (x % pri[i] == 0) 
                x /= pri[i];
        }
    }
    if (x > 1) res = res / x * (x-1);
    return res % md;
}
_ a, b, c;
int main() {
    getpri();
    a.a[1][1] = b.a[1][2] = b.a[3][1] = b.a[3][3] = 1;
    a.a[1][3] = 2; b.a[1][1] = 3;
    while (cin >> n >> m) {
        md = n * m;
        ans = 0;
        b.a[2][1] = md - 1;
        for (ll i = 1; i * i <= n; i++) {
            if (n % i == 0) {
                c = a * ksm(b, i-1);
                ans = (ans + ksc(phi(n/i), c.a[1][1])) % md;
                if (i * i < n) {
                    c = a * ksm(b, n/i-1);
                    ans = (ans + ksc(phi(i), c.a[1][1])) % md;
                }
            }
        }
        cout << ans / n << '\n';
    }
}