BZOJ3309 DZY Loves Math

题目大意

对于正整数 \(n\)，定义 \(f(n)\) 为 \(n\) 所包含质因子的最大幂指数。例如 \(f(1960)=f(2^3 \times 5^1 \times 7^2)=3\)，\(f(10007)=1\)，\(f(1)=0\)。

给定正整数 \(a,b\)，求下式的值：

\[\sum^{a}_{i=1}\sum^{b}_{j=1}f(\gcd(i,j)) \]

推导

首先记 \(n=\min(a,b),m=\max(a,b)\)。

\[\begin{aligned} & \ \ \ \ \ \ \sum^{a}_{i=1}\sum^{b}_{j=1}f(\gcd(i,j)) \\ &= \sum^{n}_{d=1}f(d)\sum^{n}_{i=1}\sum^{m}_{j=1}\left[ \gcd(i,j)=d \right] \\ &= \sum^{n}_{d=1}f(d)\sum^{\left\lfloor \frac{n}{d} \right\rfloor}_{i=1}\sum^{\left\lfloor \frac{m}{d} \right\rfloor}_{j=1}\left[ \gcd(i,j)=1 \right] \\ &= \sum^{n}_{d=1}f(d)\sum^{\left\lfloor \frac{n}{d} \right\rfloor}_{i=1}\sum^{\left\lfloor \frac{m}{d} \right\rfloor}_{j=1}\sum_{t \mid \gcd(i,j)}\mu(t) \\ &= \sum^{n}_{d=1}f(d)\sum^{\left\lfloor \frac{n}{d} \right\rfloor}_{i=1}\sum^{\left\lfloor \frac{m}{d} \right\rfloor}_{j=1}\sum_{t \mid i \land t \mid j}\mu(t) \\ &= \sum^{n}_{d=1}f(d) \sum^{\left\lfloor \frac{n}{d} \right\rfloor}_{t =1}\mu(t) \sum^{\left\lfloor \frac{n}{d} \right\rfloor}_{i=1} \left[ t \mid i \right] \sum^{\left\lfloor \frac{m}{d} \right\rfloor}_{j=1} \left[ t \mid j \right] \\ &= \sum^{n}_{d=1}f(d)\sum^{\left\lfloor \frac{n}{d} \right\rfloor}_{t=1}\mu(t) \left\lfloor \frac{n}{dt} \right\rfloor \left\lfloor \frac{m}{dt} \right\rfloor \\ \end{aligned} \]

设 \(T=dt\)（十分常用的技巧），那么有

\[= \sum^{n}_{T=1} \left\lfloor \frac{n}{T} \right\rfloor \left\lfloor \frac{m}{T} \right\rfloor \sum_{d \mid T}\mu(\frac{T}{d})f(d) \]

记 \(h=f*\mu\)，那么有

\[= \sum^{n}_{T=1} \left\lfloor \frac{n}{T} \right\rfloor \left\lfloor \frac{m}{T} \right\rfloor h(T) \]

---

那么，现在的问题是如何求 \(h\)？考虑求 \(h(n)\)。

我们可以把 \(n\) 进行质因数分解：

\[n=\prod^{m}_{i=1}p_i^{c_i} \]

考虑 \(h(T)\) 的原本写法：

\[h(T)=\sum_{d \mid T}\mu(\frac{T}{d})f(d) \]

记 \(T=\prod^{m}_{i=1}p_i^{a_i}\)，\(d=\prod^{m}_{i=1}p_i^{b_i}\)。

考虑 \(\mu(n)\) 的定义，当 \(n\) 中含有平方质因子时 \(\mu(n)=0\)，即不会在 \(h\) 中产生贡献，可以不考虑。

所以 \(\dfrac{T}{d}\) 的质因数要满足最大幂指数小于 \(2\)，即 \(a_i-b_i \in \left\{ 0,1 \right\}\)。

所以将 \(h(T)\) 改写为如下形式时

\[\displaystyle h(T) = \sum_{ab = T} f(a) \mu(b) \]

有

\[b=\prod^{m}_{i=1}p_i^{d_i}(d_i \in \left\{ 0,1 \right\}) \]

记 \(l=\max^m\limits_{i=1}c_i,k=\sum^{m}\limits_{i=1} \left[c_i=l \right]\)。即 \(l\) 为 \(n\) 所包含质因子的最大幂指数，\(k\) 为 \(n\) 所包含质因子幂指数中最大幂指数的个数。

我们发现 \(f(a)\) 的取值只有 \(l\) 和 \(l-1\) 两种可能（\(b\) 可能把最大幂指数的都取走，导致 \(f(a)\) 的取值少了 \(1\)）。

我们先按 \(k=m\) 和 \(k \neq m\) 两种情况分类讨论，在每一项讨论中，我们再分 \(f(a)=l\) 和 \(f(a)=l-1\) 两种子情况。

---

当 \(k=m\) 时，贡献为（加号前为 \(l\) 的情况，加号后为 \(l-1\) 的情况）

\[\sum^{m-1}_{s=0}(l) \times (-1)^s \times {m \choose s} + (l-1) \times (-1)^m \times {m \choose m} \]

\(f(a)=l\) 的情况不会把最大幂指数的质因数都取走。

所以我们枚举到 \(m-1\)，中间的 \((-1)^s\) 和 \((-1)^m\) 实质上就是莫比乌斯函数，最右边的二项式系数是枚举选取的方案。

将上面的式子加号右边的部分拆开，把带 \(l\) 的部分合并到左面，得到

\[\sum^{m}_{s=0}(l) \times (-1)^s \times {m \choose s} + 1 \times (-1)^m \]

二项式定理，得

\[-1 \times (-1)^m=(-1)^{m+1} \]

---

当 \(k \neq m\) 时，

当 \(f(a)=l\) 时，设在 \(k\) 个 \(c_i=l\) 的质数中选了 \(t\) 个，另外 \(m-k\) 个质数中选了 \(s\) 个，贡献为：

\[\sum^{m-k}_{s=0}\sum^{k-1}_{t=0} {k \choose t} \times l \times (-1)^{s+t} \times {m-k \choose s} \]

中间的 \((-1)^{s+t}\) 仍是莫比乌斯函数。因为 \(f(a)=l\)，所以 \(k\) 个质数不能被选完，因此枚举到 \(k-1\)。

左右拆开，分别二项式定理，得

\[\begin{aligned} \sum^{k-1}_{t=0} {k \choose t} \times l \times (-1)^t \times \sum^{m-k}_{s=0}(-1)^s \times 1 ^{m-k-s} \times {m - k \choose s} =0 \end{aligned} \]

当 \(f(a)=l-1\) 时，\(k\) 个 \(c_i=l\) 的质数一定全部被选择，设在另外 \(m-k\) 个质数中选择 \(s\) 个，贡献为：

\[\begin{aligned} & \ \ \ \ \ \sum^{m-k}_{s=0}(l-1) \times (-1)^k \times (-1)^s \times {m-k \choose s} \\ &= \sum^{m-k}_{s=0}(l-1) \times (-1)^k \times (-1)^s \times 1^{m-k-s} \times {m-k \choose s} \\ &= (l-1) \times (-1)^k \times \sum^{m-k}_{s=0} (-1)^s \times 1^{m-k-s} \times {m-k \choose s} \\ &= 0 \end{aligned} \]

---

所以最终 \(h(n)\) 的表达式为

\[h(n)= \begin{cases} (-1)^{m + 1} & k = m \\ 0 & k \neq m \end{cases} \]

Code

#include <bits/stdc++.h>

using namespace std;

#define int long long

const int N = 1e7;
const int M = 10500;

int T,n,m;

bool p[N + M];
int pri[N >> 2],cnt,low[N + M];
int h[N + M],a[N + M];

void Seive() {
    p[1] = 1;

    for(int i = 2;i <= N; i++) {
        if(!p[i]) {
            cnt ++;
            low[i] = pri[cnt] = i;
            a[i] = h[i] = 1;
        }

        for(int j = 1;j <= cnt && i * pri[j] <= N; j++) {
            p[i * pri[j]] = 1;

            if(i % pri[j] == 0) {
                a[i * pri[j]] = a[i] + 1;
                low[i * pri[j]] = low[i] * pri[j];

                if(i == low[i]) 
                    h[i * pri[j]] = 1;
                else if(a[i / low[i]] == a[i * pri[j]])
                    h[i * pri[j]] = -h[i / low[i]];
                else
                    h[i * pri[j]] = 0;

                break;
            }

            a[i * pri[j]] = 1;
            low[i * pri[j]] = pri[j];
            if(a[i] == 1)
                h[i * pri[j]] = -h[i];
            else
                h[i * pri[j]] = 0;
        }
    }

    for(int i = 1;i <= N; i++)
        h[i] += h[i - 1];
    
    return ;
}

signed main() {
    ios::sync_with_stdio(false);
    Seive();
    
    cin >> T;

    while(T --> 0) {
        cin >> n >> m;

        if(n > m)
            swap(n,m);

        int l = 1,r = 0,ans = 0;

        while(l <= n) {
            r = min(n / (n / l),m / (m / l));
            ans += (n / l) * (m / l) * (h[r] - h[l - 1]);
            l = r + 1;
        }

        cout << ans << "\n";
    }
    return 0;
}