1. 程式人生 > >[SDOI2015] 約數個數和

[SDOI2015] 約數個數和


Portal

考慮這樣一個式子:
\[ d(ij) = \sum_{x | i}\sum_{y | j} [x \bot y] \]
怎麼證明? 一開始我們一定會想到\[d(ij) = \sum_{x | i} \sum_{y | i} 1 \] 但這樣會計算重複. 於是我們考慮:
\[ d(ij) = \sum_{x | i}\sum_{\frac{m}{y} | j}1 \]
這樣每個因數就變成\[ \frac{xm}{y} \], 如果x和y不互質. 那麼就會有\[ \frac{(xp)m}{yp} == \frac{xm}{p} \]

如果xp,yp同時是i, j的因數那就會算重複.

所以一定要求此二者互質.

相應的:
\[ \sigma(ij) = \sum_{x | i}\sum_{y | j} [x \bot y]x \frac{}{} \frac{i}{y} \]


那麼我們可以開始化式子了!!
\[ Ans = \sum_{i}\sum_{j} \sum_{x | i}\sum_{y | j}[(x, y) == 1] \\ = \sum_{x} \sum_{y}[(x, y) == 1]\sum_{i} \sum_{j} [x | i][y | j] \\ = \sum_{x} \sum_{y}\sum_{d | (x, y)} \mu(d) \lfloor\frac{n}{x} \rfloor\lfloor\frac{m}{y} \rfloor\\ = \sum_{d} \mu(d) \sum_{x}\sum_{y}[d | x][d | y] \lfloor\frac{n}{x} \rfloor\lfloor\frac{m}{y} \rfloor \\ = \sum_{d} \mu(d) \sum_{x}^{n / d}\lfloor\frac{n}{xd} \rfloor\sum_{y}^{m / d}\lfloor\frac{m}{yd} \rfloor\\ \]


\(F(n) = \sum_{i = 1}^{n} \sigma_0(i)\)

那麼有
\[ Ans = \sum_{d} \mu(d) F(\frac{n}{d})F(\frac{m}{d}) \]
然後線性篩/杜教篩Min25篩洲閣篩篩以下約數個數的字首和就可以做了.

Codes

#include<bits/stdc++.h>
using namespace std;
#define rep(i, a, b) for(int i = (a), i##_end_ = (b); i <= i##_end_; ++i)
#define drep(i, a, b) for(int i = (a), i##_end_ = (b); i >= i##_end_; --i)
#define clar(a, b) memset((a), (b), sizeof(a))
#define debug(...) fprintf(stderr, __VA_ARGS__)
typedef long long LL;
typedef long double LD;
int read() {
    char ch = getchar();
    int x = 0, flag = 1;
    for (;!isdigit(ch); ch = getchar()) if (ch == '-') flag *= -1;
    for (;isdigit(ch); ch = getchar()) x = x * 10 + ch - 48;
    return x * flag;
}
void write(int x) {
    if (x < 0) putchar('-'), x = -x;
    if (x >= 10) write(x / 10);
    putchar(x % 10 + 48);
}

const int Maxn = 50009;
int prime[Maxn], isnprime[Maxn], mu[Maxn], prefixMu[Maxn], tot, low[Maxn];
LL sigma0[Maxn];
void linearSieve() {
    mu[1] = 1; sigma0[1] = 1;
    rep (i, 2, Maxn - 1) {
        if (!isnprime[i]) mu[i] = -1, prime[++tot] = i, low[i] = i, sigma0[i] = 2;
        for (int k, j = 1; j <= tot && (k = prime[j] * i) < Maxn; ++j) {
            isnprime[k] = 1;
            if (i % prime[j] == 0) {
                mu[k] = 0; low[k] = low[i] * prime[j];
                sigma0[k] = (low[i] == i) ? (sigma0[i] + 1) : (sigma0[k / low[k]] * sigma0[low[k]]);
                break;
            } else {
                mu[k] = -mu[i]; low[k] = prime[j];
                sigma0[k] = sigma0[i] * sigma0[prime[j]];
            }
        }
    }
    rep (i, 1, Maxn - 1) {
        prefixMu[i] = prefixMu[i - 1] + mu[i];
        sigma0[i] += sigma0[i - 1];
    }
}

void init() { linearSieve(); }

void solve() {
    int T = read();
    while (T--) {
        int n = read(), m = read();
        LL ans = 0; int Limit = min(n, m);
        for (int l = 1, r; l <= Limit; l = r + 1) {
            r = min(Limit, min(n / (n / l), m / (m / l)));
            ans += (prefixMu[r] - prefixMu[l - 1] * 1ll) * sigma0[n / l] * sigma0[m / l];
        }
        printf("%lld\n", ans);
    }
}

int main() {
    freopen("BZOJ3994.in", "r", stdin);
    freopen("BZOJ3994.out", "w", stdout);

    init();
    solve();

#ifdef Qrsikno
    debug("\nRunning time: %.3lf(s)\n", clock() * 1.0 / CLOCKS_PER_SEC);
#endif
    return 0;
}