Solution -「校内题」矩阵求和
Description
共 \(T\) 组数据。对于每组数据，给定 \(a, b, n\)，求 \(\sum_{i = 1}^{n} \sum_{j = 1}^{n} \gcd(a^i - b^i, a^j - b^j)\)，并取一个让人迷惑的模数 \(10^9 + 9\)。
题目保证 \(\gcd(a, b) = 1\) 且 \(1 \leq T \leq 10, 1 \leq n \leq 10^7, 1 \leq a < b \leq 10^9\)。
Solution
因为笔者太蒻赛时并没有 A 掉此题，故先感谢 lihan 巨巨在赛后讨论中对于时间复杂度优化上（引入整除分块）给我的启发。 /fad
以下为正题，所有未知量均处于正整数域。
Trick 1
试证明：\(\gcd(a^i - b^i, a^j - b^j) = a^{\gcd(i, j)} - b^{\gcd(i, j)}\)，其中 \(\gcd(a, b) = 1, a > b\)。
首先明确一条显然的引理 \(1\)：若 \(a^i = b^i\)，则一定有 \((a^i)^k = (b^i)^k\)。
以及另一条由因式定理可证的显然的引理 \(2\)：\(a^i - b^i\) 一定有一个因式为 \(a - b\)。
引理 \(2\) 成立本质上是因为 \(a = b\) 是 方程 \(a^i - b^i = 0\) 的一个通解。故可结合引理 \(1\) 可推广为 \(a^p = b^p\) 是方程 \(a^q - b^q = 0\) 的一个通解，其中 \(p\) 是 \(q\) 的因数。
故可得 \(a^i - b^i\) 一定有一个因式为 \(a^{q_1} - b^{q_1}\)，其中 \(q_1\) 是 \(i\) 的因数，同理 \(a^j - b^j\) 一定有一个因式为 \(a^{q_2} - b^{q_2}\)，其中 \(q_2\) 是 \(j\) 的因数。
则 \(a^{\gcd(i, j)} - b^{\gcd(i, j)} \mid \gcd(a^i - b^i, a^j - b^j)\)。
接下来尝试证明：\(\gcd(a^i - b^i, a^j - b^j) \mid a^{\gcd(i, j)} - b^{\gcd(i, j)}\)。
引入引理 \(3\)：设 \(c = \gcd(a^i - b^i, a^j - b^j)\)，则 \(\gcd(c, b) = 1\)。
对于引理 \(3\)，若 \(\gcd(c, b) \neq 1\)，则有质数 \(p \mid \gcd(c, b)\)，故有 \(p \mid \gcd(c, b) \mid b, p \mid \gcd(c, b) \mid c \mid (a^i - b^i)\)。
可得 \(p \mid a^i\)，即 \(p \mid a\)，所以 \(p \mid \gcd(a, b)\)。此与 \(\gcd(a, b) = 1\) 矛盾，故引理 \(3\) 成立。
记 \(d = \gcd(i, j)\)，则由裴蜀定理可得，存在 \(x, y\)，使得 \(x \times i - y \times j = \gcd(i, j) = d\)。
又因为 \(c \mid (a^j - b^j)\)。故可由引理 \(2\) 得 \(c \mid (a^{j \times y} - b^{j \times y})\)。所以 \(c \mid a^d \times (a^{j \times y} - b^{j \times y}) = (a^{x \times i} - a^{d} \times b^{j \times y})\)。
又由 \(c \mid (a^i - b^i)\)，可知 \(c \mid (a^{x \times i} - b^{x \times i})\)，所以有 \(c \mid (a^{d} \times b^{j \times y} - b^{x \times i}) = b^{j \times y} \times (a^{d} - b^{d})\)。
而 \(\gcd(b, c) = 1\)，故有 \(c \mid (a^{d} - b^{d})\)，即 \(\gcd(a^i - b^i, a^j - b^j) \mid a^{\gcd(i, j)} - b^{\gcd(i, j)}\)。
又因为 \(a^{\gcd(i, j)} - b^{\gcd(i, j)} \mid \gcd(a^i - b^i, a^j - b^j)\)。
故有 \(\gcd(a^i - b^i, a^j - b^j) = a^{\gcd(i, j)} - b^{\gcd(i, j)}\)，得证。
Trick 2
试证明：若 \(f(x) = \sum_{i = 1}^{n} \sum_{j = 1}^{n} [\gcd(i, j) = x]\) ，则有 \(f(x) = 2 \times \left( \sum_{i = 1}^{\lfloor \frac {n}{x} \rfloor} \varphi(i) \right) - 1\)。
首先明确一条引理：对于任意 \(i, j\)，若 \(\gcd(i, j) = 1\)，则有 \(\gcd(i \times x, j \times x) = x\)。
于是我们考虑统计给定范围内有哪些组 \(i, j\) 满足要求。
对于 \(i\)，满足条件的 \((i, j)\) 有 \(\varphi(i)\) 组，其中 \(1 < i \leq \lfloor \frac {n} {x} \rfloor, i > j\)。
且不难发现 \((i, j)\) 与 \((j, i)\) 需要分开计算贡献，即 \(2 \times \left( \sum_{i = 2}^{\lfloor \frac {n}{x} \rfloor} \varphi(i) \right)\)。
最后特殊考虑 \(i = 1\) 的情况，仅有 \((1, 1)\) 符合条件，故 \(f(x) = 2 \times \left( \sum_{i = 2}^{\lfloor \frac {n}{x} \rfloor} \varphi(i) \right) + 1\)。
因为我们规定 \(\varphi(1) = 1\)，故也有 \(f(x) = 2 \times \left( \sum_{i = 1}^{\lfloor \frac {n}{x} \rfloor} \varphi(i) \right) - 1\)。得证。
Trick 3
根据前两个 Trick 可得：
但分析发现此题并不能使用 \(O(T \times n \log(n))\) 的做法 A 掉。
观察答案式子发现 \(\lfloor \frac {n}{x} \rfloor\) 的形式，故使用整除分块。
对于 \(l, r\)，有 \(\lfloor \frac {n}{l} \rfloor = \lfloor \frac {n}{r} \rfloor,\lfloor \frac {n}{l} \rfloor < \lfloor \frac {n}{r + 1} \rfloor\)。则可知其对答案的贡献为：
记 \(p = 2 \times \left( \sum_{i = 1}^{\lfloor \frac {n}{l} \rfloor} \varphi(i) \right) - 1\)，再套用等比数列求和公式，故上式可化为：\(\frac {a^{r + 1} - a^l} {a - 1} \times p - \frac {b^{r + 1} - b^l} {b - 1} \times p\)。
然后其他的就没有什么难点了，线性筛求一个 \(\varphi\) 的前缀和，快速幂求答案与逆元，然后跑整除分块即可。
Code
#include <cstdio>
typedef long long LL;
LL Max(LL x, LL y) { return x > y ? x : y; }
LL Min(LL x, LL y) { return x < y ? x : y; }
LL Abs(LL x) { return x < 0 ? -x : x; }
int read() {
int k = 1, x = 0;
char s = getchar();
while (s < '0' || s > '9') {
if (s == '-')
k = -1;
s = getchar();
}
while (s >= '0' && s <= '9') {
x = (x << 3) + (x << 1) + s - '0';
s = getchar();
}
return x * k;
}
LL read_LL() {
int k = 1;
LL x = 0;
char s = getchar();
while (s < '0' || s > '9') {
if (s == '-')
k = -1;
s = getchar();
}
while (s >= '0' && s <= '9') {
x = (x << 3) + (x << 1) + s - '0';
s = getchar();
}
return x * k;
}
void write(LL x) {
if (x < 0) {
putchar('-');
x = -x;
}
if (x > 9)
write(x / 10);
putchar(x % 10 + '0');
}
void print(LL x, char s) {
write(x);
putchar(s);
}
const int mod = 1e9 + 9;
const int MAXN = 1e7 + 5;
const int MAXM = 664579 + 5;
LL Quick_Pow(LL a, LL b) {
LL res = 1;
while (b) {
if (b & 1)
res = res * a % mod;
a = a * a % mod;
b >>= 1;
}
return res;
}
bool flag[MAXN];
LL phi[MAXN], sum[MAXN];
int num[MAXM], len = 0;
void Euler(int n) {
flag[1] = true;
phi[1] = 1;
for (int i = 2; i <= n; i++) {
if (!flag[i]) {
num[++len] = i;
phi[i] = i - 1;
}
for (int j = 1; j <= len; j++) {
if (i * num[j] > n)
break;
flag[i * num[j]] = true;
if (i % num[j] == 0) {
phi[i * num[j]] = phi[i] * num[j] % mod;
break;
} else
phi[i * num[j]] = phi[i] * phi[num[j]] % mod;
}
}
sum[1] = 1;
for (int i = 2; i <= n; i++) sum[i] = (sum[i - 1] + phi[i]) % mod;
return;
}
LL calc(int l, int r, int x) { return (Quick_Pow(x, r + 1) - Quick_Pow(x, l) + mod) % mod; }
int main() {
Euler(MAXN - 5);
int t = read();
while (t--) {
LL a = read_LL(), b = read_LL(), ans = 0;
LL inva = Quick_Pow(a - 1, mod - 2) % mod, invb = Quick_Pow(b - 1, mod - 2) % mod;
int n = read(), l, r;
for (l = 1; l <= n; l = r + 1) {
r = n / (n / l);
LL x = ((sum[n / l] << 1) % mod - 1 + mod) % mod;
ans = (ans + calc(l, r, a) * inva % mod * x % mod) % mod;
ans = (ans - calc(l, r, b) * invb % mod * x % mod + mod) % mod;
}
print(ans, '\n');
}
return 0;
}
