时间戳

首页 > 技术文章 > 「SDOI2015」约数个数和

「SDOI2015」约数个数和

luckyblock 2020-09-10 09:35 原文

知识点: 莫比乌斯反演

原题面：Loj Luogu

题意简述

\(T\) 次询问，每次询问给定 \(n,m\)。
定义 >\(\operatorname{d}(i)\) 为 \(i\) 的约数个数，求：

\[\sum_{i=1}^{n}\sum_{j=1}^m\operatorname{d}(ij) \]

\(1<T,n\le 5\times 10^4\)。

分析题意

一个结论：

\[\operatorname{d}(ij) = \sum_{x|i}\sum_{y|j}[\gcd(x,y)=1] \]

证明

先考虑 \(i = p^a\)，\(j=p^b(p\in \mathrm{primes})\) 的情况，有：

\[\operatorname{d}(p^{a+b})=\sum_{x|p^a}\sum_{y|p^b}[\gcd(x,y)=1] \]

对于等式左侧，\(p^{a+b}\) 的约数个数为 \(a+b+1\)。
对于等式右侧，在保证 \(\gcd(x,y)=1\) 成立的情况下，有贡献的数对 \((x,y)\) 只能是下列三种形式：

  • \(x>0,y-0\)，\(x\) 有 \(a\) 种取值方法。
  • \(x=0,y>0\)，\(y\) 有 \(b\) 种取值方法。
  • \(x=0,y=0\)。

则等式右侧贡献次数为 \(a+b+1\) 次，等于 \(p^{a+b}\) 的约数个数。
则当 \(i = p^a\)，\(j=p^b(p\in \mathrm{primes})\) 时等式成立。

又不同质因子间相互独立，上述情况可拓展到一般的情况。

对 \(\operatorname{d}(i,j)\) 进行化简，代入 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\)，原式等价于：

\[\begin{aligned} \operatorname{d}(ij) =& \sum_{x|i}\sum_{y|j}[\gcd(x,y)=1]\\ =& \sum_{x|i}\sum_{y|j}\sum\limits_{d\mid \gcd(x,y)} {\mu (d)} \end{aligned}\]

调换枚举顺序，先枚举 \(d\)，原式等价于：

\[\sum_{d=1}[d|i][d|j]{\mu (d)}\sum_{x|i}[d|x]\sum_{y|j}[d|y] \]

把各项中的 \(d\) 提出来，消去整除的条件，原式等价于：

\[\begin{aligned} &\sum_{d=1}[d|i][d|j]{\mu (d)}\sum_{x|\frac{i}{d}}\sum_{y|\frac{j}{d}}1\\ =& \sum_{d=1}[d|i][d|j]{\mu (d)}\cdot \operatorname{d}(\dfrac{i}{d})\operatorname{d}(\dfrac{j}{d}) \end{aligned}\]

将 \(\operatorname{d}(ij)\) 代回原式，原式等价于：

\[\sum_{i=1}^{n}\sum_{j=1}^m \sum_{d=1}[d|i][d|j]{\mu (d)}\cdot \operatorname{d}(\dfrac{i}{d})\operatorname{d}(\dfrac{j}{d}) \]

调换枚举顺序，先枚举 \(d\)，原式等价于：

\[\sum_{d=1}{\mu (d)}\sum_{i=1}^{n}[d|i]\sum_{j=1}^m [d|j]\cdot \operatorname{d}(\dfrac{i}{d})\operatorname{d}(\dfrac{j}{d}) \]

把 \(i,j\) 中的 \(d\) 提出来，变为枚举 \(\frac{i}{d}, \frac{j}{d}\)，消去整除的条件，原式等价于：

\[\sum_{d=1}{\mu (d)}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\operatorname{d}(i)\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}\operatorname{d}(j) \]

考虑预处理 \(S(x) = \sum\limits_{i=1}^{x}\operatorname{d}(i)\)，则原式等价于：

\[\sum_{d=1}{\mu (d)}S\left({\left\lfloor\frac{n}{d}\right\rfloor}\right)S\left({\left\lfloor\frac{m}{d}\right\rfloor}\right) \]

线性筛预处理 \(\mu,\operatorname{d}\)，数论分块求解即可，复杂度 \(O(n+T\sqrt{n})\)。

爆零小技巧

注意筛法中出现平方因子的处理。

代码实现

//知识点：莫比乌斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstring>
#define ll long long
const int kMaxn = 5e4 + 10;
//=============================================================
int pnum, p[kMaxn];
ll mu[kMaxn], num[kMaxn], d[kMaxn]; //num 为最小质因子的次数
ll summu[kMaxn], sumd[kMaxn];
bool vis[kMaxn];
//=============================================================
inline int read() {
  int f = 1, w = 0;
  char ch = getchar();
  for (; !isdigit(ch); ch = getchar())
    if (ch == '-') f = -1;
  for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
  return f * w;
}
void Getmax(int &fir_, int sec_) {
  if (sec_ > fir_) fir_ = sec_;
}
void Getmin(int &fir_, int sec_) {
  if (sec_ < fir_) fir_ = sec_;
}
void Euler(int lim_) {
  vis[1] = true;
  mu[1] = d[1] = 1ll;
  for (int i = 2; i <= lim_; ++ i) {
    if (! vis[i]) {
      p[++ pnum] = i;
      mu[i] = - 1;
      num[i] = 1;
      d[i] = 2;
    }
    for (int j = 1; j <= pnum && i * p[j] <= lim_; ++ j) {
      vis[i * p[j]] = true;
      if (i % p[j] == 0) { //勿忘平方因子
        mu[i * p[j]] = 0;
        num[i * p[j]] = num[i] + 1;
        d[i * p[j]] = 1ll * d[i] / num[i * p[j]] * (num[i * p[j]] + 1ll);
        break;
      }
      mu[i * p[j]] = - mu[i];
      num[i * p[j]] = 1;
      d[i * p[j]] = 2ll * d[i]; //

    }
  }
  for (int i = 1; i <= lim_; ++ i) {
    summu[i] = summu[i - 1] + mu[i];
    sumd[i] = sumd[i - 1] + d[i];
  }
}

//=============================================================
int main() { 
  Euler(kMaxn - 10);
  int T = read();
  while (T --) {
    int n = read(), m = read(), lim = std :: min(n, m);
    ll ans = 0ll;
    for (int l = 1, r; l <= lim; l = r + 1) {
      r = std :: min(n / (n / l), m / (m / l));
      ans += 1ll * (summu[r] - summu[l - 1]) * sumd[n / l] * sumd[m / l]; //
    }
    printf("%lld\n", ans);
  }
  return 0;
}
/*
1
32 43
*/
/*
15420
*/

推荐阅读