SDOI 2015 约束个数和

Description:

\(T \le 5 \times 10^4\)组询问, 每组询问给定\(n\)\(m\), 请你求出
\[ \sum_{i = 1}^n \sum_{j = 1}^m \sigma_0 (ij) \]
Solution:

先给出一个结论:
\[ \sigma_0(ij) = \sum_{a | i} \sum_{b | j} [\gcd(a, b) = 1] \]
证明: 我们令\(i = p_1^{a_1} p_2^{a_2} \cdots\), \(j = p_1^{b_1} p_2^{b_2} \cdots\), \(d | ij\)\(d = p_1^{c_1} p_2^{c_2} \cdots\), 则\(c_n \le a_n + b_n\).

考虑如何不重复地统计每一个\(d\): 令\(c_n = A_n + B_n\), 其中\(A_n\)\(B_n\)分别为\(i\)\(j\)\(c_n\)的贡献, 则我们要求
\[ \begin{cases} B_n = 0 & A_n < a_n \\ B_n \ge 0 & A_n = a_n \end{cases} \]
这样一来, \(c_n\)的表示形式就变成唯一的了, 因而不会被重复统计. 我们再考虑如何统计这样的\(A_n\)\(B_n\): 我们令\(A_n' = a_n - A_n\), 则约束条件变为
\[ \begin{cases} B_n = 0 & A_n' \ne 0 \\ B_n \ge 0 & A_n' = 0 \end{cases} \]
等价于\(\gcd(A_n', B_n) = 1\).

因此得证.

好吧, 假如看不懂上面的这一些证明, 就这么想吧: \(i\)表示\(a\)中不取多少, \(j\)表示\(b\)中取多少, 只要保证\(\gce(a, b) = 1\), 即不会重复统计.

因此我们考虑原题的式子
\[ \begin{aligned} \sum_{i = 1}^n \sum_{j = 1}^m \sigma_0(ij) &= \sum_{i = 1}^n \sum_{j = 1}^m \sum_{a | i} \sum_{b | j} [\gcd(a, b) = 1] \\ &= \sum_{i = 1}^n \sum_{j = 1}^m \sum_{a | i} \sum_{b | j} \sum_{d | \gcd(a, b)} \mu(d) \\ &= \sum_{i = 1}^n \sum_{j = 1}^m \sum_{d | \gcd(i, j)} \mu(d) \sigma_0(\frac i d) \sigma_0(\frac j d) \\ &= \sum_{d = 1}^n \sum_{i = 1}^{\lfloor \frac n d \rfloor} \sum_{j = 1}^{\lfloor \frac m d \rfloor} \mu(d) \sigma_0(i) \sigma_1(j) \\ &= \sum_{d = 1}^n \mu(d) \sum_{i = 1}^{\lfloor \frac n d \rfloor} \sigma_0(i) \sum_{j = 1}^{\lfloor \frac m d \rfloor} \sigma_0(j) \end{aligned} \]
分块处理后半部分即可.

时间复杂度: 预处理\(O(n)\), 单次询问\(O(n^\frac 1 2)\)

#include <cstdio>
#include <cctype>
#include <cstring>
#include <algorithm>

using namespace std;
namespace Zeonfai
{
    inline int getInt()
    {
        int a = 0, sgn = 1; char c;
        while (! isdigit(c = getchar())) if (c == '-') sgn *= -1;
        while (isdigit(c)) a = a * 10 + c - '0', c = getchar();
        return a * sgn;
    }
}
const int N = (int)5e4, MOD = (int)1e9;
typedef int arr[N + 7];
typedef long long Larr[N + 7];
int tot;
arr isNotPrime, prm, mu, minDivisor, minDivisorDegree, sgm;
Larr a, b, c;
inline void initialize()
{
    memset(isNotPrime, 0, sizeof(isNotPrime));
    tot = 0;
    sgm[1] = mu[1] = 1;
    for (int i = 2; i <= N; ++ i)
    {
        if (! isNotPrime[i])
        {
            prm[tot ++] = i;
            mu[i] = -1;
            minDivisor[i] = i;
            minDivisorDegree[i] = 1;
            sgm[i] = 2;
        }
        for (int j = 0; j < tot && i * prm[j] <= N; ++ j)
        {
            int x = i * prm[j]; isNotPrime[x] = 1;
            if (i % prm[j])
            {
                mu[x] = - mu[i];
                minDivisor[x] = prm[j];
                minDivisorDegree[x] = 1;
                sgm[x] = sgm[i] * 2;
            }
            else
            {
                mu[x] = 0;
                minDivisor[x] = minDivisor[i] * prm[j];
                minDivisorDegree[x] = minDivisorDegree[i] + 1;
                sgm[x] = sgm[i / minDivisor[i]] * (minDivisorDegree[x] + 1);
            }
        }
    }
    a[0] = b[0] = c[0] = 0;
    for (int i = 1; i <= N; ++ i) a[i] = a[i - 1] + sgm[i], b[i] = a[i] * a[i], c[i] = c[i - 1] + mu[i];
}
int main()
{
    using namespace Zeonfai;
    initialize();
    int T = getInt();
    for (int cs = 0; cs < T; ++ cs)
    {
        int n = getInt(), m = getInt();
        long long ans = 0;
        int L = 1;
        while (L <= min(n, m))
        {
            int R = min(n / (n / L), m / (m / L));
            ans = ans + a[n / L] * a[m / L] * (c[R] - c[L - 1]);
            L = R + 1;
        }
        printf("%lld\n", ans);
    }
}

转载于:https://www.cnblogs.com/ZeonfaiHo/p/7623831.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值