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) 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 可得：

[sum_{i = 1}^{n} sum_{j = 1}^{n} gcd(a^i - b^i, a^j - b^j) = sum_{i = 1}^{n} sum_{j = 1}^{n} a^{gcd(i, j)} - b^{gcd(i, j)} = sum_{x = 1}^{n} (a^{x} - b^{x}) imes f(x)= sum_{x = 1}^{n} (a^{x} - b^{x}) imes left(2 imes left( sum_{i = 1}^{lfloor frac {n}{x} floor} varphi(i) ight) - 1 ight) ]

但分析发现此题并不能使用 (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)。则可知其对答案的贡献为：

[(a^l + a^{l + 1} + dots + a^r) imes left(2 imes left( sum_{i = 1}^{lfloor frac {n}{l} floor} varphi(i) ight) - 1 ight) - (b^l + b^{l + 1} + dots + b^r) imes left(2 imes left( sum_{i = 1}^{lfloor frac {n}{l} floor} varphi(i) ight) - 1 ight) ]

记 (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;
}