【题解】BZOJ 3309 DZY Loves Math

传送门
令 \(f(n)\) 为 \(n\) 的质因子的最大幂指数，求 \(\sum_{i=1}^n\sum_{j=1}^m f(\gcd(i,j))\)，\(T\) 组询问。
\(T \leq 10^4 \space n, m \leq 10^7\)
题解：
这里设 \(n < m\)

\[\sum_{i=1}^n\sum_{j=1}^m f(\gcd(i,j)) \\ = \sum_{d=1}^n\sum_{i=1}^{\lfloor \frac nd \rfloor}\sum_{j=1}^{\lfloor \frac md \rfloor} f(d) [i \bot j] \\ = \sum_{d=1}^n f(d) \sum_{i=1}^{\lfloor \frac nd \rfloor}\sum_{j=1}^{\lfloor \frac md \rfloor} \sum_{k|i, k|j} \mu(k) \\ = \sum_{d=1}^n f(d) \sum_{k=1}^{\lfloor \frac nd \rfloor} \mu(k) \lfloor \frac n{kd} \rfloor \lfloor \frac m{kd} \rfloor \\ = \sum_{T=1}^n \lfloor \frac n{T} \rfloor \lfloor \frac m{T} \rfloor \sum_{d|T} f(T) \mu(\frac Td) \]

前面已经可以整数分块了，考虑后面怎么处理
考虑后面这个 \(\sum_{d|T} f(\frac Td) \mu(T)\) 的值可以是什么
先给出结论：
设 \(T\) 的标准分解形式为 \(p_1^{a_1}p_2^{a_2}\cdots p_k^{a_k}\)
如果 \(a_1=a_2=\cdots =a_k\) 那么值为 \((-1)^{(k+1)}\)
否则值为 \(0\)
考虑如何证明
这里只考虑 \(\mu(d)\) 不为 \(0\) 的情况，此时 \(d\) 的每个质因子的指数必然为 \(0\) 或 \(1\)（这里是指 \(d\) 写成 \(T\) 的标准分解形式）
显然有 \(2^{k-1}\) 个取值是 \(1\)，\(2^{k-1}\) 个取值是 \(-1\)
如果 \(T\) 的指数存在不同，那么考虑
无论 \(f(\frac Td)\) 的值为什么，都可以钦定一个不是最大值的 \(a_i\) 让这个 \(a_i-1\) 来使得这样 \(f(\frac Td')\mu(d')\) 正好为此时 \(f(\frac Td)\mu(d)\) 的相反数，这样两两抵消，最后的值就是 \(0\)
否则 \(T\) 的指数 \(a_i\) 的值都是一样的
当且仅当 \(\frac Td=p_1p_2\dots p_k\) 时 \(f(\frac Td)\) 会比别的情况少 \(1\)
那么如果此时 \(\mu(d)\) 为 \(-1\) 时（即 \(k\) 为奇数），那么贡献的负数的绝对值就会比正数少 \(1\)，最后这个式子的值为 \(1\)，反之亦然
得到了这个结论，考虑如何线筛。
我们可以用 \(b_i\) 记录 \(i\) 这个数最小质因子的指数
用 \(a_i\) 记录 \(i\) 的最小质因子的 \(b_i\) 次方的值
用 \(h_i\) 记录上式的值
这样就可以线筛上面的式子及其前缀和了

Code

#include <queue>
#include <cstdio>
#include <cctype>
#include <cstring>
#include <algorithm>
using namespace std;
#define FOR(i,j,k) for(int i=j; i<=k; ++i)
#define ROF(i,j,k) for(int i=j; i>=k; --i)
inline int read (void) {
  int x = 0, f = 1, ch = getchar();
  while(!isdigit(ch)) { if(ch == '-') f = -f; ch = getchar(); }
  while(isdigit(ch)) { x = x * 10 + ch - '0'; ch = getchar(); }
  return x * f;
}
const int maxn = 10000000;
const int maxl = 670000;
bool np[maxn+5];
int l, p[maxl];
int a[maxn+5], b[maxn+5], h[maxn+5];
int main (void) {
  np[0] = np[1] = 1;
  FOR(i,2,maxn) {
    if(!np[i]) p[++l] = a[i] = i, b[i] = h[i] = 1;
    for(int j=1; j<=l&&i*p[j]<=maxn; ++j) {
      int k = i * p[j]; np[k] = 1;
      if(i%p[j] == 0) {
        a[k] = a[i] * p[j]; b[k] = b[i] + 1;
        if(k == a[k]) h[k] = 1;
        else if(b[k] == b[k / a[k]]) h[k] = - h[k / a[k]];
        else h[k] = 0;
        break;
      } else {
        a[k] = p[j]; b[k] = 1;
        if(b[i] == 1) h[k] = -h[i];
        else h[k] = 0;
      }
    }
  }
  FOR(i,2,maxn) h[i] += h[i-1];
  int t = read();
  while(t--) {
    int n = read(), m = read();
    if(n > m) swap (n, m);
    long long ans = 0;
    FOR(i,1,n) {
      int r = min(n / (n / i), m / (m / i));
      ans += 1ll * (n / i) * (m / i) * (h[r] - h[i-1]);
      i = r;
    }
    printf("%lld\n", ans);
  }
  return 0;
}