BZOJ2671 Calc [莫比乌斯反演]

原创

qq62c30ac77b2a7 2022-07-05 12:16:54 博主文章分类：莫比乌斯反演 ©著作权

文章标签 i++ 整除 #define 文章分类 虚拟化云计算

©著作权归作者所有：来自51CTO博客作者qq62c30ac77b2a7的原创作品，请联系作者获取转载授权，否则将追究法律责任

传送门

设 d = gcd(a, b), a = ud, b = vd

$a+b|ab\Leftrightarrow d(u + v)|d^2uv\rightarrow u + v|duv$

因为 gcd(u, v) = 1, 所以 gcd(u+v, v) = gcd(u+v, u) = 1

BZOJ2671 Calc [莫比乌斯反演]_i++_02

$ans=\sum_d\sum_u\sum_{v<u} [u + v|d,ud<=n][(u,v)=1]$

$=\sum_u\sum_{v<u} \left \lfloor \frac{n}{u(u+v)}\right \rfloor[(u,v)=1]$

$=\sum_{d=1}^n \mu(d) \sum_{u}\sum_{v<u}\left \lfloor \frac{n}{u(u+v)d^2}\right \rfloor$

然后卡一下上界

$=\sum_{d=1}^{\sqrt n} \mu(d) \sum_{u}^{\sqrt{n/d^2}}\sum_{k= u+1}^{ukd^2<=n}\left \lfloor \frac{n}{ukd^2}\right \rfloor$

枚举 k 的时候再整除分块就可以直接怼过去

#include<bits/stdc++.h>
#define N 100050
using namespace std;
typedef long long ll;
ll n; int prim[N], isp[N], mu[N], tot;
void prework(int m){
  mu[1] = 1;
  for(int i = 2; i <= m; i++){
    if(!isp[i]) prim[++tot] = i, mu[i] = -1;
    for(int j = 1; j <= tot; j++){
      if(i * prim[j] > m) break;
      isp[i * prim[j]] = 1;
      if(i % prim[j] == 0){ mu[i * prim[j]] = 0; break;}
      mu[i * prim[j]] = mu[i] * mu[prim[j]];
    }
  }
}
ll calc(ll n, ll m){
  ll sum = 0;
  for(int u = 1; u <= m; u++){
    int t = n / u, up = (u << 1) - 1;
    for(int l = u + 1, r = 0; l <= t, l <= up; l = r + 1){
      int v = t / l; if(!v) break;
      r = min(t / v, up); sum += 1ll * (r - l + 1) * v; 
    }
  } return sum;
}
int main(){
  scanf("%d", &n); int m = sqrt(n); prework(m);
  ll ans = 0; 
  for(ll i = 1; i <= m; i++) if(mu[i]) ans += mu[i] * calc(n / i / i, m / i);
  cout << ans; return 0;
}