Statement

给定整数 n,an,a。

定义『好的珠子』由一个无序三元组 (x1,x2,x3)(x_1,x_2,x_3) 表示，且满足 1≤xi≤a1\le x_i\le a，gcd⁡{x1,x2,x3}=1\gcd\{x_1,x_2,x_3\}=1。

定义『好的项链』是一串由 nn 个珠子构成的环，满足相邻两个珠子不同。

试统计本质不同（即循环不同构）的好的项链的个数，答案对 109+710^9+7 取模。

TT 组数据。数据范围：n≤1014n\le10^{14}，a≤107a\le10^7，T≤10T\le10。时限 3 s3~\text{s}。

Solution

二合一题。

统计好的珠子个数

考虑容斥。设 f(x)f(x) 表示 gcd=x\text{gcd}=x 的无序三元组个数，g(x)g(x) 表示 gcd⁡\gcd 为 xx 的倍数的无序三元组个数。

显然有 g(d)=∑d∣nf(x)g(d)=\sum_{d|n}f(x)，由莫比乌斯反演得 f(d)=∑d∣nμ(xd)⋅g(x)f(d)=\sum_{d|n}\mu(\frac xd)\cdot g(x)。所以 f(1)=∑d=1aμ(d)⋅g(d)f(1)=\sum_{d=1}^a\mu(d)\cdot g(d)。

求 g(d)g(d) 考虑 Polya 引理，显然有 11 种置换有 33 个环，33 种置换有 22 个环，22 种置换有 11 个环。所以设 m=⌊ad⌋m=\left\lfloor\frac ad\right\rfloor，则 g(d)=16(m3+3m2+2m)g(d)=\frac16(m^3+3m^2+2m)。整除分块即可。

统计好的项链个数

考虑 Burnside 引理，设 f(x)f(x) 表示长为 xx 的相邻、首尾颜色不同的序列个数，则：

ans=1n∑i=1nf(gcd⁡{n,i})=1n∑d=1n∑i=1n[gcd⁡{n,i}=d]⋅f(d)=1n∑d∣nφ(nd)⋅f(d)\begin{aligned} \text{ans}&=\frac1n\sum_{i=1}^nf(\gcd\{n,i\})\\ &=\frac1n\sum_{d=1}^n\sum_{i=1}^n[\gcd\{n,i\}=d]\cdot f(d)\\ &=\frac1n\sum_{d|n}\varphi\left(\frac nd\right)\cdot f(d)\\ \end{aligned}

考虑 f(d)f(d) 的递推式。记 cc 表示颜色个数，分类讨论：

  • 如果 11 和 d−1d-1 颜色相同，方案数为 f(d−2)⋅(c−1)f(d-2)\cdot(c-1)。
  • 否则颜色不同，方案数为 f(d−1)⋅(c−2)f(d-1)\cdot(c-2)。

所以得到：

f(x)={0,x=1c⋅(c−1),x=2f(x−2)⋅(c−1)+f(x−1)⋅(c−2),x≥3f(x)=\begin{cases} 0,&x=1\\ c\cdot(c-1),&x=2\\ f(x-2)\cdot(c-1)+f(x-1)\cdot(c-2),&x\ge3 \end{cases}

但 nn 太大了，递推显然不行。考虑求出通项：

f(x)=f(x−2)⋅(c−1)+f(x−1)⋅(c−2)f(x)+f(x−1)=f(x−2)⋅(c−1)+f(x−1)⋅(c−1)\begin{aligned} f(x)&=f(x-2)\cdot(c-1)+f(x-1)\cdot(c-2)\\ f(x)+f(x-1)&=f(x-2)\cdot(c-1)+f(x-1)\cdot(c-1)\\ \end{aligned}

设 g(x)=f(x)+f(x−1)g(x)=f(x)+f(x-1)，则有：

g(x)=(c−1)⋅g(x−1)=(c−1)x−2g(2)=(c−1)x−2f(2)\begin{aligned} g(x)&=(c-1)\cdot g(x-1)\\ &=(c-1)^{x-2}g(2)\\ &=(c-1)^{x-2}f(2) \end{aligned}

然后 f(x)=g(x)−g(x−1)+g(x−2)−g(x−3)+⋯f(x)=g(x)-g(x-1)+g(x-2)-g(x-3)+\cdots，等比数列求和即可。

得到通项是 f(x)=(c−1)x+(−1)x⋅(c−1)f(x)=(c-1)^x+(-1)^x\cdot(c-1)。代回上面 ans\text{ans} 的式子即可。

实现时注意两点：

  • 因为 nn 太大了，DFS 按唯一分解定理枚举质因数的幂次，来枚举约数 nd\frac nd 并维护 φ(d)\varphi(d)。
  • 有一种烦人的情况是，可能 (109+7)∣n(10^9+7)|n。所以把所有模数改成 (109+7)2(10^9+7)^2 做，这样能保证余数不变。因为此时 ans\text{ans} 一定是 nn 的倍数，最后再分别除掉 109+710^9+7 和 n109+7\frac n{10^9+7} 即可。注意 66 对 (109+7)2(10^9+7)^2 的逆元要用欧拉定理求。

Code

实现的时候，偷懒可以用 __int128。不然只能龟速乘了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
#include <bits/stdc++.h>
namespace IO {
char gc() {
    static char buf[100003], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100003, stdin), p1 == p2) ? EOF : *p1++;
}
template <typename Tp>
void input(Tp& x) {
    static bool f;
    static char c;
    for (f = 0; !isdigit(c = gc());)
        if (c == '-') f = 1;
    for (x = 0; isdigit(c); c = gc()) x = (x << 1) + (x << 3) + (c ^ '0');
    if (f) x = -x;
}
const char emptyc = '\n';
template <typename Tp>
void output(Tp x, const char c = emptyc) {
    static char f[97];
    static int pf;
    pf = 0;
    if (!x) return putchar('0'), putchar(c), void();
    if (x < 0) x = -x, putchar('-');
    while (x) f[pf++] = x % 10 + '0', x /= 10;
    while (pf) putchar(f[--pf]);
    putchar(c);
}
}  // namespace IO
using namespace IO;
typedef __int128_t LL;
const int N = 1e7 + 3;
const LL P = 1e9 + 7, mod = P * P;
LL qpow(LL a, LL b) {
    LL s = 1;
    while (b) {
        if (b & 1) s = s * a % mod;
        a = a * a % mod, b >>= 1;
    }
    return s;
}
int A;
LL n, c, inv6;
namespace Sieve {
bool flag[N];
int mu[N], smu[N];
std::vector<int> prime;
void init(int n) {
    mu[1] = smu[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!flag[i]) prime.emplace_back(i), mu[i] = -1;
        smu[i] = smu[i - 1] + mu[i];
        for (int j : prime) {
            if (j > n / i) break;
            flag[i * j] = 1;
            if (i % j == 0) break;
            mu[i * j] = -mu[i];
        }
    }
}
}  // namespace Sieve
namespace CalcColor {
void main() {
    using Sieve::smu;
    c = 0;
    for (int l = 1, r; l <= A; l = r + 1) {
        LL m = A / l, g = (m * m * m % mod + 3 * m * m % mod + 2 * m % mod) % mod;
        r = A / m, c = (c + g * (smu[r] - smu[l - 1] + mod) % mod) % mod;
    }
    c = c * inv6 % mod;
}
}  // namespace CalcColor
namespace CalcChain {
LL ans;
std::vector<std::pair<LL, int>> factor;
LL f(LL x) { return (qpow(c - 1, x) + (x & 1 ? mod + 1 - c : c - 1)) % mod; }
void dfs(int x, LL cur, LL phi) {
    ans = (ans + phi % mod * f(n / cur) % mod) % mod;
    for (int i = x; i < (int)factor.size(); i++) {
        LL cp = 1, pp = 1, p = factor[i].first;
        for (int j = 0; j < factor[i].second; j++) pp *= p - !j, cp *= p, dfs(i + 1, cur * cp, phi * pp);
    }
}
void main() {
    ans = 0, factor.clear();
    LL now = n;
    for (LL i : Sieve::prime)
        if (now % i == 0) {
            factor.emplace_back(i, 0);
            while (now % i == 0) factor.back().second++, now /= i;
        }
    if (now != 1) factor.emplace_back(now, 1);
    dfs(0, 1, 1);
    if (n % P)
        ans = qpow(n % P, P - 2) % P * ans % P;
    else
        ans = qpow(n / P, P - 2) % P * (ans / P) % P;
    output(ans);
}
}  // namespace CalcChain
int main() {
    int T;
    input(T), Sieve::init(1e7), inv6 = qpow(6, P * (P - 1) - 1);
    while (T--) input(n), input(A), CalcColor::main(), CalcChain::main();
    return 0;
}