BZOJ3309: DZY Loves Math

Description

　　
​　　 对于正整数n，定义f(n)为n所含质因子的最大幂指数。例如f(1960)=f(2^3 * 5^1 * 7^2)=3, f(10007)=1, f(1)=0。
​　　 给定正整数a,b，求\(∑_{i=1}^a∑_{j=1}^bf(gcd(i,j))\)

。
　　

Input

　　
　　第一行一个数T，表示询问数。
　　接下来T行，每行两个数a,b，表示一个询问。
　　

Output

　　
​　　 对于每一个询问，输出一行一个非负整数作为回答。
　　

Sample Input

　　
​　　 4
​　　 7558588 9653114
​　　 6514903 4451211
​　　 7425644 1189442
　　​ 6335198 4957
　　

Sample Output

　　
​　　 35793453939901
​　　 14225956593420
​　　 4332838845846
　　​ 15400094813
　　

HINT

　　
【数据规模】
　　T<=10000
​　　1<=a,b<=10^7
　　
　　
　　

Solution

自我感觉很有难度的一道莫比乌斯反演题，是我做过的里面对筛法的分析涉及最深的了。
注：下文中，\((i,j)\)指\(gcd(i,j)\)，\(n,m\)指题目描述中的\(a,b\)，这里均设\(n<m\)

\[\begin{aligned} 原式&=\sum_{d=1}^nf(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\sum_{k|(i,j)}\mu(k)\\ &=\sum_{d=1}^nf(d)\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{i=1}^{\lfloor\frac{n}{kd}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{kd}\rfloor}\mu(k)\\ &=\sum_{d=1}^nf(d)\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}\lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor\mu(k)\\ &设T=kd\\ &=\sum_{T=1}^n\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\sum_{k|T}f(k)\mu(\frac{T}{k}) \end{aligned} \]

考虑筛出\(\sum_{k|T}f(k)\mu(\frac{T}{k})\)

• \(T=1\)
  \(\sum_{k|T}f(k)\mu(\frac{T}{k})=f(1)\mu(1)=0\)
• \(T=p(p \in prime)\)
  \(\sum_{k|T}f(k)\mu(\frac{T}{k})=f(T)\mu(1)+f(1)\mu(T)=1+0=1\)
• \(T=i*p(p \in prime)\)
  考虑根据唯一分解定理将\(T和\frac{T}{k}\)分解，\(T=p_1^{c_1}p_2^{c_2}...p_k^{c_k}\)，\(\frac{T}{k}=p_1^{a_1}p_2^{a_2}...p_k^{a_k}\)
  我们设\(g(T)=\sum_{k|T}f(k)\mu(\frac{T}{k})\)，记\(n=i*p(p \in prime)\)
  显然只有当\(a_i=1或0\)时\(f(d)\)对答案有贡献，所以只讨论\(a_i=1或0\)时的\(n\)

\[\begin{aligned} g(n)&=\sum_{k|n且f(k)=f(n)}f(n)\mu(\frac{n}{k})+\sum_{k|n且f(d)\not =f(n)}f(d)\mu(\frac{n}{k})\\ &=f(n)\sum_{k|n且f(k)=f(n)}\mu(\frac{n}{k})+(f(n)-1)\sum_{k|n且f(d)\not =f(n)}\mu(\frac{n}{k})\\ &=f(n)\sum_{k|n}\mu(\frac{n}{k})-\sum_{k|n且f(d)\not =f(n)}\mu(\frac{n}{k})\\ &=f(n)*[n=1]-\sum_{k|n且f(d)\not =f(n)}\mu(\frac{n}{k})\\ &=-\sum_{k|n且f(d)\not =f(n)}\mu(\frac{n}{k})\\ \end{aligned} \]

所以，如果\(\mu(\frac{n}{k})\)对答案有贡献，则对于每一项指数最大的\(p_i\)，\(k\)中\(p_i\)的指数一定需要是\(n\)中\(p_i\)的指数\(-1\).
分两种情况讨论：

1. \(c_1=c_2...=c_k\)，那么显然为了满足\(f(d)\not = f(n)\)的条件，所有的\(a_i\)都必须为\(1\)，所以有$$g(n)=-\mu(p_1p_2p_3...p_k)=(-1)^{k+1}$$
2. \(c_i\)不完全相等。则我们按\(a_i\)将\(p_i\)分入集合\(A,B\)
   \(A=\{i|a_i=1\},B=\{i|a_i=0\}\)
   则\(A\)中的\(p_i\)必须取，\(B\)中的\(p_i\)取不取都无所谓。所以有：

\[\begin{aligned} g(n)&=-\sum\mu((\prod_{i\in A}p_i)*(\prod_{i\in B}p_i或1))\\ &=-\sum\mu(\prod_{i\in A}p_i)\mu(\prod_{i\in B}p_i或1)\\ &=-\mu(\prod_{i\in A}p_i)\sum\mu(\prod_{i\in B}p_i或1)\\ &=(-1)^{|A|+1}\sum\mu(\prod_{i\in B}p_i或1) \end{aligned} \]

如何处理集合\(B\)呢？
实际上我们就是在\(|B|\)个数里面，取出若干个数的积，由于\(\mu=\{1,-1,0\}\)且这里不会取到\(0\)，所以我们可以把它抽象成\(n\)个数，我们选的所有方案的价值之和，选奇数个数的方案价值为\(-1\)，选偶数个数的方案价值为\(1\)
所以其实等价于\(\sum_{i=0}^nC_n^i(-1)^i=0\)。
这个玩意怎么证明呢？

\(\sum_{i=0}^nC_n^i(-1)^i=0\)的证明：
奇数时显然成立。\((C_n^x=C_n^{n-x})\)
偶数时，考虑杨辉三角的递推式\(C_m^n=C_m^{n-1}+C_{m-1}^{n-1}\)，会发现对于杨辉三角的偶数行，它奇数和和偶数和都正好相当于上一行（奇数行）\(\sum_{i=1}^{n-1}C^i_n\)，所以奇偶数的正负正好抵消。
证毕。

所以有

\[ \begin{aligned} g(n)&=(-1)^{|A|+1}\sum\mu(\prod_{i\in B}p_i或1)\\ &=(-1)^{|A|+1}*0\\ &=0 \end{aligned} \]

所以说，对答案有贡献的\(n\)，当且仅当\(c_1=c_2...=c_k\)
这个东西要筛出来就挺容易的了。

具体做法：

维护一个\(a_i\)表示数\(n\)的最小质因子\(p_{min}\)的指数，\(b_i=p_{min}^{a_i}\)。 记\(x=i*p,g_i表示上文的函数g(i)\)

• 当筛到\(i*p(p\in prime)\)且\(p⊥i\)时（\(p⊥i\)即\(p\)与\(i\)互质）
  \(a_x=1\)，\(b_x=p\)（因为线性筛保证每次筛掉这个数的都是它的最小质因子）
  当\(a_x=a_i\)时，\(g_x=-g_i\)，否则\(g_x=0\)
• 当筛到\(i*p(p\in prime)\)且\(p|i\)时，
  \(a_x=a_i+1,b_x=b_i*p\)
  这时\(b\)数组就派上用场，因为\(i\)的最小质因子也是\(p\)，所以并不能比较\(a_i\)和\(a_p\)，我们设\(d=\frac{i}{b_i}\)，则有\(d⊥i,d⊥x\)，那么这时当\(a_d=a_x\)时，\(g_x=-g_d\)，否则\(g_x=0\)

那么数论分块一下，就可以在\(O(n+T\sqrt{n})\)的复杂度内解决本题了。

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N = 1e7 + 5;
int p[N], vis[N], cnt;
ll a[N], b[N];
ll g[N];
void init() {
    g[1] = 0;
    for(int i = 2; i < N; ++i) {
        if(!vis[i]) {
            p[++cnt] = i;
            b[i] = i;
            a[i] = 1;
            g[i] = 1;
        }
        for(int j = 1; j <= cnt && i * p[j] < N; ++j) {
            vis[i * p[j]] = 1;
            if(i % p[j] == 0) {
                int d = i / b[i];
                a[i * p[j]] = a[i] + 1;
                b[i * p[j]] = b[i] * p[j];
                if(d == 1) g[i * p[j]] = 1;
                else g[i * p[j]] = (a[i * p[j]] == a[d]) ? -g[d] : 0;
                break;
            }
            a[i * p[j]] = 1;
            b[i * p[j]] = p[j];
            g[i * p[j]] = (a[i] == 1) ? -g[i] : 0;
        }
    }
    for(int i = 2; i < N; ++i) g[i] += g[i - 1];
}
int main() {
    init();
    int T; scanf("%d", &T);
    while(T--) {
        ll n, m, ans = 0; 
        scanf("%lld%lld", &n, &m);
        if(n > m) swap(n, m);
        for(ll l = 1, r; l <= n; l = r + 1) {
            r = min(n / (n / l), m / (m / l));
            ans += (ll)(n / l) * (m / l) * (g[r] - g[l - 1]);
        }
        printf("%lld\n", ans);
    }
}