[国家集训队]Crash的数字表格

知识点: 莫比乌斯反演

原题面:Luogu


取模,很有精神!


题意简述

给定 \(n,m\) , 求:

\[\sum_{i=1}^n\sum_{j=1}^{m} \operatorname{lcm}(i,j)\bmod 20101009 \]

\(1\le n,m\le 10^7\)


分析题意

易知原式等价于:

\[\sum_{i=1}^{n}\sum_{j=1}^{m}\dfrac{ij}{\gcd (i,j)} \]

考虑枚举 \(\gcd(i,j)\),设枚举量为 \(d\)
\(d=\gcd(i,j)\) 的充要条件是满足 \(d|i, d|j\)\(\gcd(\dfrac{i}{d},\dfrac{j}{d}) = 1\),则原式等价于:

\[\sum_{i=1}^n\sum_{j=1}^m \sum_{d=1} \dfrac{ij}{d}[d|i][d|j][\gcd(\frac{i}{d}, \frac{j}{d})=1] \]

先枚举 \(d\),则原式等价于:

\[\sum_{d=1}\sum_{i=1}^{n}[d\mid i]\sum_{j=1}^m [d\mid j][\gcd(\dfrac{i}{d}, \dfrac{j}{d}=1)] \dfrac{ij}{d} \]

这个 \(d\) 很烦人,把 \(i,j\) 中的 \(d\) 提出来,变为枚举 \(\frac{i}{d}\)\(\frac{j}{d}\)
消去 \(d\mid i\)\(d\mid j\) 的限定条件,则原式等价于:

\[\begin{aligned} &\sum_{d=1}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[\gcd(i,j)=1]\dfrac{id\times jd}{d}\\ =& \sum_{d=1}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[\gcd(i,j)=1]ijd\\ =& \sum_{d=1}d \sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[\gcd(i,j)=1]ij \end{aligned}\]

单独考虑后半部分,设 \(f(x,y) = \sum\limits_{i=1}^{x} \sum\limits_{j=1}^{y}[\gcd(i,j)=1]ij\)
发现 \(f(x,y)\) 的形式与 [HAOI2011] Problem b 中的式子类似,代入 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\) 套路一波:

\[\begin{aligned} f(x,y) =& \sum\limits_{i=1}^{x} \sum\limits_{j=1}^{y}[\gcd(i,j)=1]ij\\ =& \sum_{i=1}^{x} \sum_{j=1}^{y}\sum_{d\mid \gcd(i,j)}\mu(d) ij\\ =& \sum_{d=1}\mu(d)\sum_{i=1}^{x}[d\mid i] \sum_{j=1}^{y}[d\mid j] ij\\ =& \sum_{d=1}\mu(d) d^2\sum_{i=1}^{\left\lfloor \frac{x}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{y}{d}\right\rfloor}ij \end{aligned}\]

前半部分 \(\sum\limits_{d=1}\mu(d) d^2\),可以考虑筛出 \(\mu(d)\) 后求前缀和。
后半部分是等差数列乘等差数列的形式,设 \(g(p,q) = \sum\limits_{i=1}^{p} \sum\limits_{j=1}^{q}ij\)\(g_{p,q}\) 可以通过下式 \(O(1)\) 计算:

\[g(p,q) = \sum_{i=1}^{p} i \sum_{j=1}^{q}j= \dfrac{(1 + p)\times p}{2}\times \dfrac{(1+q)\times q}{2} \]

则对于 \(f(x,y)\),有:

\[f(x,y) = \sum_{d=1}\mu(d) d^2\cdot g(\left\lfloor \frac{x}{d}\right\rfloor, \left\lfloor\frac{y}{d}\right\rfloor) \]

数论分块求解即可。

再看回原式,原式等价于:

\[\sum_{d=1}d\cdot f(\left\lfloor\frac{n}{d}\right\rfloor, \left\lfloor\frac{m}{d}\right\rfloor) \]

又是一个可以数论分块求解的形式。
线性筛预处理后 数论分块套数论分块,复杂度 \(O(n + m)\),瓶颈是线性筛。


爆零小技巧

处理时会出现求平方的运算,需要特别注意取模问题,ll 都会爆,被坑惨了。

在预处理前缀和的这个地方:

sum[i] = (sum[i - 1] + 1ll * i * i % kMod * (mu[i] + kMod)) % kMod; //仅令 mu + kMod

注意先对平方取模,在乘上 \(\mu\),否则就会爆掉。
以及可以仅令 \(mu + mod\)

以及这个地方:

int g(int n_, int m_) {
  return (1ll * n_ * (n_ + 1ll) / 2ll % kMod) * (1ll * m_ * (m_ + 1ll) / 2ll % kMod) % kMod;
}

平方计算,注意随时取模。


代码实现

//莫比乌斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstring>
#define ll long long
const ll kMod = 20101009;
const int kMaxn = 1e7 + 10;
//=============================================================
int pnum, p[kMaxn];
ll mu[kMaxn], sum[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] = 1ll;
  for (int i = 2; i <= lim_; ++ i) {
    if (! vis[i]) {
      p[++ pnum] = i;
      mu[i] = - 1;
    }
    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;
        break;
      }
      mu[i * p[j]] = - mu[i];
    }
  }
  sum[1] = 1ll;
  for (int i = 1; i <= lim_; ++ i) {
    sum[i] = (sum[i - 1] + 1ll * i * i % kMod * (mu[i] + kMod)) % kMod; //仅令 mu + kMod
  }
}
int g(int n_, int m_) {
  return (1ll * n_ * (n_ + 1ll) / 2ll % kMod) * (1ll * m_ * (m_ + 1ll) / 2ll % kMod) % kMod;
}
int f(int n_, int m_) {
  int lim = std :: min(n_, m_), ret = 0;
  for (int l = 1, r; l <= lim; l = r + 1) {
    r = std :: min(n_ / (n_ / l), m_ / (m_ / l));
    ret = (ret + 1ll * (sum[r] - sum[l - 1] + kMod) * g(n_ / l, m_ / l) % kMod) % kMod;
  }
  return ret;
}
int Sum(ll l, ll r) {
  return (1ll * (r - l + 1ll) * (l + r) / 2ll) % kMod;
}
//=============================================================
int main() { 
  int n = read(), m = read();
  int lim = std :: min(n, m), ans = 0;
  Euler(lim);
  for (int l = 1, r; l <= lim; l = r + 1) {
    r = std :: min(n / (n / l), m / (m / l));
    ans = (ans + 1ll * Sum(l, r) * f(n / l, m / l) % kMod) % kMod;
  }
  printf("%d", ans);
  return 0;
}
/*
7718820 8445880
*/
posted @ 2020-09-10 08:33  Luckyblock  阅读(101)  评论(0编辑  收藏  举报