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);
    }
}