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) eq 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 imes i - y imes j = gcd(i, j) = d)。
又因为 (c mid (a^j - b^j))。故可由引理 (2) 得 (c mid (a^{j imes y} - b^{j imes y}))。所以 (c mid a^d imes (a^{j imes y} - b^{j imes y}) = (a^{x imes i} - a^{d} imes b^{j imes y}))。
又由 (c mid (a^i - b^i))，可知 (c mid (a^{x imes i} - b^{x imes i}))，所以有 (c mid (a^{d} imes b^{j imes y} - b^{x imes i}) = b^{j imes y} imes (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 imes left( sum_{i = 1}^{lfloor frac {n}{x} floor} varphi(i) ight) - 1)。
首先明确一条引理：对于任意 (i, j)，若 (gcd(i, j) = 1)，则有 (gcd(i imes x, j imes x) = x)。
于是我们考虑统计给定范围内有哪些组 (i, j) 满足要求。
对于 (i)，满足条件的 ((i, j)) 有 (varphi(i)) 组，其中 (1 < i leq lfloor frac {n} {x} floor, i > j)。
且不难发现 ((i, j)) 与 ((j, i)) 需要分开计算贡献，即 (2 imes left( sum_{i = 2}^{lfloor frac {n}{x} floor} varphi(i) ight))。
最后特殊考虑 (i = 1) 的情况，仅有 ((1, 1)) 符合条件，故 (f(x) = 2 imes left( sum_{i = 2}^{lfloor frac {n}{x} floor} varphi(i) ight) + 1)。
因为我们规定 (varphi(1) = 1)，故也有 (f(x) = 2 imes left( sum_{i = 1}^{lfloor frac {n}{x} floor} varphi(i) ight) - 1)。得证。
Trick 3
根据前两个 Trick 可得：
但分析发现此题并不能使用 (O(T imes n log(n))) 的做法 A 掉。
观察答案式子发现 (lfloor frac {n}{x} floor) 的形式，故使用整除分块。
对于 (l, r)，有 (lfloor frac {n}{l} floor = lfloor frac {n}{r} floor,lfloor frac {n}{l} floor < lfloor frac {n}{r + 1} floor)。则可知其对答案的贡献为：
记 (p = 2 imes left( sum_{i = 1}^{lfloor frac {n}{l} floor} varphi(i) ight) - 1)，再套用等比数列求和公式，故上式可化为：(frac {a^{r + 1} - a^l} {a - 1} imes p - frac {b^{r + 1} - b^l} {b - 1} imes 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, '
');
}
return 0;
}