1. 程式人生 > 實用技巧 >「SDOI2015」約數個數和

「SDOI2015」約數個數和

知識點: 莫比烏斯反演

原題面: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
*/