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 可得:

\[\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}) \times f(x)= \sum_{x = 1}^{n} (a^{x} - b^{x}) \times \left(2 \times \left( \sum_{i = 1}^{\lfloor \frac {n}{x} \rfloor} \varphi(i) \right) - 1 \right) \]

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

\[(a^l + a^{l + 1} + \dots + a^r) \times \left(2 \times \left( \sum_{i = 1}^{\lfloor \frac {n}{l} \rfloor} \varphi(i) \right) - 1 \right) - (b^l + b^{l + 1} + \dots + b^r) \times \left(2 \times \left( \sum_{i = 1}^{\lfloor \frac {n}{l} \rfloor} \varphi(i) \right) - 1 \right) \]

\(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;
}
posted @ 2021-08-26 09:37  STrAduts  阅读(74)  评论(0编辑  收藏  举报