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)

[egin{aligned} 原式&=sum_{d=1}^nf(d)sum_{i=1}^{lfloorfrac{n}{d} floor}sum_{j=1}^{lfloorfrac{m}{d} floor}sum_{k|(i,j)}mu(k)\ &=sum_{d=1}^nf(d)sum_{k=1}^{lfloorfrac{n}{d} floor}sum_{i=1}^{lfloorfrac{n}{kd} floor}sum_{j=1}^{lfloorfrac{m}{kd} floor}mu(k)\ &=sum_{d=1}^nf(d)sum_{k=1}^{lfloorfrac{n}{d} floor}lfloorfrac{n}{kd} floorlfloorfrac{m}{kd} floormu(k)\ &设T=kd\ &=sum_{T=1}^nlfloorfrac{n}{T} floorlfloorfrac{m}{T} floorsum_{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)

[egin{aligned} g(n)&=sum_{k|n且f(k)=f(n)}f(n)mu(frac{n}{k})+sum_{k|n且f(d) ot =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) ot =f(n)}mu(frac{n}{k})\ &=f(n)sum_{k|n}mu(frac{n}{k})-sum_{k|n且f(d) ot =f(n)}mu(frac{n}{k})\ &=f(n)*[n=1]-sum_{k|n且f(d) ot =f(n)}mu(frac{n}{k})\ &=-sum_{k|n且f(d) ot =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) ot = 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)取不取都无所谓。所以有：

[egin{aligned} g(n)&=-summu((prod_{iin A}p_i)*(prod_{iin B}p_i或1))\ &=-summu(prod_{iin A}p_i)mu(prod_{iin B}p_i或1)\ &=-mu(prod_{iin A}p_i)summu(prod_{iin B}p_i或1)\ &=(-1)^{|A|+1}summu(prod_{iin 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)，所以奇偶数的正负正好抵消。
证毕。

所以有

[ egin{aligned} g(n)&=(-1)^{|A|+1}summu(prod_{iin 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(pin 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(pin 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+Tsqrt{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
", ans);
    }
}