## 莫比乌斯反演与容斥原理

mobius_and_inclusion_exclusion_principle

## 【BZOJ 4635】数论小测验

### Code

#include<bits/stdc++.h>
#define LL long long
using namespace std;

const int MOD = 1000000007;
const int N = 10000009;
const int M = 32;

int n,m,SZ,s1[N],hc[N],to[1001];
int tot,pri[N>>3],mu[N]; bool vis[N];

char c=getchar(); int f=1,ret=0;
while (c<'0'||c>'9') {if(c=='-')f=-1;c=getchar();}
while (c<='9'&&c>='0') {ret=ret*10+c-'0';c=getchar();}
return ret * f;
}

inline int Pow(int w, int t) {
int ret = 1;
for (;t;t>>=1,w=(LL)w*w%MOD)
if (t&1) ret=(LL)ret*w%MOD;
return ret;
}

inline void GetMu() {
mu[1] = 1;
for (int i=2;i<N;i++) {
if (!vis[i]) mu[i] = -1, pri[++tot] = i;
for (int j=1;j<=tot&&pri[j]*i<N;j++) {
vis[i*pri[j]] = 1;
if (i%pri[j]) mu[i*pri[j]] = -mu[i];
else {mu[i*pri[j]] = 0; break;}
}
mu[i] = mu[i-1] + mu[i];
}
}

inline int cal(int mx) {
int ret = 0;
for (int l=1,r,tmp;l<=mx;l=r+1) {
r = mx / (mx / l);
tmp = (LL)(mu[r] - mu[l-1]) * hc[mx/l] % MOD;
ret = (ret + tmp) % MOD;
}
return (ret + MOD) % MOD;
}

struct Matrix{
int a[M][M];
inline Matrix() {memset(a,0,sizeof(a));}
inline Matrix(const Matrix *A) {
for (int i=0;i<M;i++) for (int j=0;j<M;j++) a[i][j] = A->a[i][j];
}
inline Matrix(int x) {
memset(a,0,sizeof(a));
for (int i=0;i<SZ;i++) a[i][i] = x;
}
inline void operator *= (const Matrix &A) {
Matrix ret;
for (int i=0,*t1,*t3;i<SZ;i++) {
t1 = ret.a[i]; const int *t2 = A.a[i];
for (int j=0;j<SZ;j++) {
t3 = a[j];
for (int k=0;k<SZ;k+=4) {
t1[k] = (t1[k] + (LL)t3[k] * t2[j]) % MOD;
t1[k+1] = (t1[k+1] + (LL)t3[k+1] * t2[j]) % MOD;
t1[k+2] = (t1[k+2] + (LL)t3[k+2] * t2[j]) % MOD;
t1[k+3] = (t1[k+3] + (LL)t3[k+3] * t2[j]) % MOD;
}
}
}
for (int i=0;i<SZ;i++) for (int j=0;j<SZ;j++) a[i][j] = ret.a[i][j];
}
inline Matrix operator ^ (int t) {
Matrix ret(1),w(this);
for (;t;t>>=1,w*=w) if (t&1) ret*=w;
return ret;
}
}ans,trans;

inline int solve(int w) {
tot = 0;
for (int i=2,tmp=w;i<=tmp;i++) {
if (tmp % i == 0) {
pri[++tot] = i; tmp /= i;
while (tmp % i == 0) pri[tot] *= i, tmp /= i;
}
}
ans = Matrix(); trans = Matrix();
ans.a[0][0] = 1; SZ = 1 << tot;
memset(to,0,sizeof(to));
for (int i=1,t=1,ww;ww=pri[i],i<=tot;i++,t<<=1)
for (int j=ww;j<=m;j+=ww) to[j] |= t;
for (int i=1,tt;tt=to[i],i<=m;i++) {
for (int p=0,ww;ww=p,p<SZ;p++)
++trans.a[p|tt][p];
}
trans = trans ^ n;
ans *= trans;
return ans.a[SZ-1][0];
}

int main() {
if (type == 1) {
GetMu();
for (int t=1,vout;vout=0,t<=T;t++) {
for (int l=1,r;l<=m;l=r+1) {
r = m / (m / l);
hc[m / l] = Pow(m / l, n);
}
for (int l=L,r,tmp;l<=R;l=r+1) {
r = m / (m / l); tmp = cal(m / l);
vout = (vout + (min(r,R)-max(L,l)+1ll) * tmp) % MOD;
}
printf("%d\n",vout);
}
} else if (type == 2) {
for (int t=1,vout,l,r;vout=0,t<=T;t++) {
for (int w=l;w<=r;w++) vout = (vout + solve(w)) % MOD;
printf("%d\n",vout);
}
}
return 0;
}


## 【SPOJ DIVCNT2】Counting Divisors (square)

### 解题报告

spoj_divcnt2_solution

## 【BZOJ 2986】Non-Squarefree Numbers

### 解题报告

prefix_sum_of_mobius_function

## 【BZOJ 3994】[SDOI2015] 约数个数和

### Code

#include<bits/stdc++.h>
#define LL long long
using namespace std;

const int N = 50000+9;

int n,m,tot,pri[N],mu[N],d[N],sum[N],f[N];
bool vis[N]; LL vout;

char c=getchar(); int f=1,ret=0;
while (c<'0'||c>'9') {if(c=='-')f=-1;c=getchar();}
while (c<='9'&&c>='0') {ret=ret*10+c-'0';c=getchar();}
return ret * f;
}

inline void prework() {
mu[1] = 1;
for (int i=2;i<N;i++) {
if (!vis[i]) pri[++tot] = i, mu[i] = -1;
for (int j=1;j<=tot&&pri[j]*i<N;j++) {
vis[i*pri[j]] = 1;
if (i % pri[j]) mu[i*pri[j]] = -mu[i];
else {mu[i*pri[j]] = 0; break;}
}
}
for (int i=1;i<N;i++)
for (int j=i;j<N;j+=i)
d[j]++;
for (int i=1;i<N;i++) {
f[i] = f[i-1] + d[i];
sum[i] = sum[i-1] + mu[i];
}
}

int main() {
prework();
if (m > n) swap(n, m);
for (int l=1,r;l<=m;l=r+1) {
r = min(n / (n / l), m / (m / l));
vout += (LL)(sum[r] - sum[l-1]) * f[n/l] * f[m/l];
}
printf("%lld\n",vout);
}
return 0;
}


## 【BZOJ 3930】[CQOI2015] 选数

### 解题报告

PoPoQQQ大爷给了一种神奇的做法：http://blog.csdn.net/popoqqq/article/details/44917831

$\sum\limits_{i = 1}^n {\mu (i) = 1 – \sum\limits_{i = 2}^n {(0 – \mu (i)) = } 1 – \sum\limits_{i = 2}^n {(\sum\limits_{j|i} {\mu (j)} – \mu (i))} = 1 – \sum\limits_{i = 2}^n {\sum\limits_{j|i,i \ne j} {\mu (j)} } } = 1 – \sum\limits_{j = 1}^n {\mu (j) \cdot (\left\lfloor {\frac{n}{j}} \right\rfloor – 1)} = 1 – \sum\limits_{j = 1}^{\left\lfloor {\frac{n}{2}} \right\rfloor } {\mu (j) \cdot (\left\lfloor {\frac{n}{j}} \right\rfloor – 1)}$

—————————— UPD 2017.2.20 ——————————

## 【BZOJ 3025】[Balkan2003] Farey数列

#include<bits/stdc++.h>
#define LL long long
#define LD long double
using namespace std;

const int N = 100000+9;
const LL EPS = 1e10;

int pri[N],tot,n,cnt,mu[N];
LL f[N],K,vx,vy;
bool vis[N];

char c=getchar(); int ret=0,f=1;
while (c<'0'||c>'9') {if(c=='-')f=-1;c=getchar();}
while (c<='9'&&c>='0') {ret=ret*10+c-'0';c=getchar();}
return ret*f;
}

inline void Get_Mu() {
mu[1] = 1;
for (int i=2;i<=n;i++) {
if (!vis[i]) pri[++tot] = i, mu[i] = -1;
for (int j=1;j<=tot&&i*pri[j]<=n;j++) {
vis[i*pri[j]] = 1;
if (i%pri[j]) mu[i*pri[j]] = -mu[i];
else {mu[i*pri[j]] = 0; break;}
}
}
}

inline LL Get(LL w) {
for (int i=1;i<=n;i++) f[i] = i * w / EPS;
for (int i=2;i<=n;i++) f[i] += f[i-1];

LL ret = 0, sta;
for (int i=1;i<=n;i++) {
sta = (LL)EPS * n / (w * i);
ret += f[min(n/(LL)i, sta)] * mu[i];
if (sta < n/i) ret += (n/i) * (n/i - sta) * mu[i];
}
return ret;
}

inline void Get_Ans(LL w) {
for (int y=1,x;y<=n;y++) {
x = min((LL)n, w * y / EPS);
if (x && (vx * y < vy * x || !vx))
vx = x, vy = y;
}
printf("%lld %lld\n", vx, vy);
}

int main(){
cin>>n>>K;
Get_Mu();

LL l=1, r=(LL)n*EPS, mid, ret;
while (l <= r) {
mid = l + r >> 1;
if (Get(mid) <= K) ret = mid, l = mid + 1;
else  r = mid - 1;
}
Get_Ans(ret);
return 0;
}


## 【SPOJ 7001】VLATTICE

#include<iostream>
#include<cstdio>
#include<bitset>
#define LL long long
using namespace std;

const int MAXN = 1000000+9;
const int MX = 1000000;

int pri[MAXN],tot,mu[MAXN];
bitset<MAXN> tag;

inline void Get_mu(){
mu[1] = 1;
for (int i=2;i<=MX;i++) {
if (!tag[i]) pri[++tot] = i, mu[i] = -1;
for (int j=1;j<=tot && i*pri[j]<=MX;j++){
tag[i*pri[j]] = 1;
if (i % pri[j]) mu[i*pri[j]] = -mu[i];
else {mu[i*pri[j]]=0; break;}
}
}
for (int i=1;i<=MX;i++) mu[i] += mu[i-1];
}

int main(){
Get_mu(); int n,T,tmp; cin>>T;
while (T--) {
scanf("%d",&n); LL vout = 0;
for (int i=1;i<=n;i=tmp+1){
tmp = n/(n/i);
vout += (LL)(mu[tmp]-mu[i-1])*(n/i)*(n/i);
} vout = vout*3+3;
for (int i=1;i<=n;i=tmp+1){
tmp = n/(n/i);
vout += (LL)(mu[tmp]-mu[i-1])*(n/i)*(n/i)*(n/i);
}
printf("%lld\n",vout);
}
return 0;
}


—– UPD 2016.9.21 —–

$$\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {\sum\limits_{k = 1}^n {[\gcd (i,j,k) = 1] = } } } \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {\sum\limits_{k = 1}^n {\sum\limits_{d|\gcd (i,j,k)} {\mu (d)} = \sum\limits_{d = 1}^n {\mu (d) \cdot {{\left\lfloor {\frac{n}{d}} \right\rfloor }^3}} } } }$$

## 【BZOJ 2820】YY的GCD

#include<iostream>
#include<cstdio>
#include<bitset>
#define LL long long
using namespace std;

const int MAXN = 10000000+9;
const int MX = 10000000;

int pri[MAXN],tot,fuc[MAXN],mu[MAXN];
bitset<MAXN> tag;

inline void Get_Fuc(){
fuc[1] = 0;
for (int i=2;i<=MX;i++){
if (!tag[i]) pri[++tot] = i, fuc[i] = 1, mu[i] = -1;
for (int j=1;j<=tot && i*pri[j]<=MX;j++){
tag[i*pri[j]] = 1;
if (i % pri[j]) fuc[i*pri[j]] = -fuc[i] + mu[i], mu[i*pri[j]] = -mu[i];
else {fuc[i*pri[j]] = mu[i]; mu[i] = 0; break;}
}
}
for (int i=2;i<=MX;i++) fuc[i] += fuc[i-1];
}

int main(){
Get_Fuc(); int T; cin>>T;
while (T--) { LL vout = 0;
int n,m,tmp; scanf("%d%d",&n,&m);
if (n > m) swap(n,m);
for (int i=1;i<=n;i=tmp+1){
tmp = min(n/(n/i),m/(m/i));
vout += (LL)(fuc[tmp]-fuc[i-1])*(n/i)*(m/i);
}
printf("%lld\n",vout);
}
return 0;
}


## 【BZOJ 2705】[SDOI2012] Longge的问题

ps:hht告诉我们，除数函数可以线筛，然而一脸懵逼QAQ

DFS-version：https://oi.abcdabcd987.com/eight-gcd-problems/（我就偷懒不写啦！）
original-version:

#include<iostream>
#include<cstdio>
#include<cmath>
#define LL long long
using namespace std;

inline LL Get_Phi(LL n){
LL ret = n;
for (int i=2,lim=ceil(sqrt(n));i<=lim;i++) if (n % i == 0) {
ret = ret*(i-1)/i;
while (n % i == 0) n /= i;
}
if (n > 1) ret = ret*(n-1)/n;
return ret;
}

int main(){
LL n,vout=0; scanf("%lld",&n);
for (int i=1,lim=ceil(sqrt(n));i<=lim;i++) if (n%i == 0) {
vout += (LL)i*Get_Phi(n/i);
if (i*i < n) vout += (LL)(n/i)*Get_Phi(i);
}
printf("%lld\n",vout);
return 0;
}


—————————— UPD 2017.4.8 ——————————

## 【BZOJ 2190】[SDOI2008] 仪仗队

#include<iostream>
#include<cstdio>
#include<bitset>
using namespace std;

const int MAXN = 40000+9;

int n,pri[MAXN],tot,phi[MAXN];
bitset<MAXN> tag;

inline void Get_Phi(){
phi[1] = 1;
for (int i=2;i<=n;i++){
if (!tag[i]) pri[++tot] = i, phi[i] = i-1;
for (int j=1;j<=tot && pri[j]*i<=n;j++) {
tag[i*pri[j]] = 1;
if (i % pri[j]) phi[i*pri[j]] = phi[i]*(pri[j]-1);
else {phi[i*pri[j]] = phi[i]*pri[j]; break;}
}
}
for (int i=2;i<=n;i++) phi[i] += phi[i-1];
}

int main(){
scanf("%d",&n); n-=1;
Get_Phi();
printf("%d\n",phi[n]*2-1+2);
return 0;
}


#include<iostream>
#include<cstdio>
#include<bitset>
#define LL long long
using namespace std;

const int MAXN = 100000+9;

int n,m,pri[MAXN],tot,tmp;
LL phi[MAXN],vout;
bitset<MAXN> tag;

inline void Get_Phi(){
phi[1] = 1;
for (int i=2;i<=m;i++){
if (!tag[i]) pri[++tot] = i, phi[i] = i-1;
for (int j=1;j<=tot && pri[j]*i<=m;j++){
tag[i*pri[j]] = 1;
if (i % pri[j]) phi[i*pri[j]] = phi[i]*(pri[j]-1);
else {phi[i*pri[j]] = phi[i]*pri[j]; break;}
}
}
for (int i=1;i<=m;i++) phi[i] += phi[i-1];
}

int main(){
scanf("%d%d",&n,&m);
if (n > m) swap(n, m); Get_Phi();
for (int i=1;i<=n;i=tmp+1){
tmp = min(n/(n/i),m/(m/i));
vout += (LL)(phi[tmp]-phi[i-1])*(n/i)*(m/i);
}
printf("%lld\n",vout*2-(LL)n*m);
return 0;
}


（x,y）这个点到原点的连线上整点的个数为gcd(x,y)

—————————— UPD 2017.2.1 ——————————

—————————— UPD 2017.7.3 ——————————

## 【BZOJ 2818】Gcd

#include<iostream>
#include<cstdio>
#include<bitset>
#define LL long long
using namespace std;

const int MAXN = 10000000+9;

int pri[MAXN],tot,mu[MAXN],n;
bitset<MAXN> tag;

inline void Get_mu(){
mu[1] = 1;
for (int i=2;i<=n;i++){
if (!tag[i]) pri[++tot] = i, mu[i] = -1;
for (int j=1;j<=tot && pri[j]*i <= n;j++){
tag[i*pri[j]] = 1;
if (i % pri[j]) mu[i*pri[j]] = -mu[i];
else {mu[i*pri[j]] = 0; break;}
}
}
for (int i=1;i<=n;i++) mu[i] += mu[i-1];
}

int main(){
scanf("%d",&n); Get_mu(); LL vout = 0;
for (int j=1;j<=tot && pri[j] <= n;j++) {
int k = pri[j], a=n/k;
for (int i=1,tmp;i<=a;i=tmp+1)
tmp = a/(a/i),
vout += (LL)(mu[tmp]-mu[i-1])*(a/i)*(a/i);
}
printf("%lld\n",vout);
return 0;
}


#include<iostream>
#include<cstdio>
#include<bitset>
#define LL long long
using namespace std;

const int MAXN = 10000000+9;

int pri[MAXN],tot,n;
LL phi[MAXN],vout=0;
bitset<MAXN> tag;

inline void Get_Phi(){
phi[1] = 1;
for (int i=2;i<=n;i++){
if (!tag[i]) pri[++tot] = i, phi[i] = i-1;
for (int j=1;j<=tot && pri[j]*i<=n;j++){
tag[pri[j]*i] = 1;
if (i % pri[j]) phi[i*pri[j]] = phi[i]*(pri[j]-1);
else {phi[i*pri[j]] = phi[i]*pri[j]; break;}
}
}
for (int i=1;i<=n;i++) phi[i] += phi[i-1];
}

int main(){
scanf("%d",&n); Get_Phi();
for (int i=1;i<=tot && pri[i] <= n;i++)
vout += phi[n/pri[i]]*2 - 1;
printf("%lld\n",vout);
return 0;
}


## 【BZOJ 2005】[Noi2010] 能量采集

%JCVB，看了半小时，还是不懂QAQ

#include<iostream>
#include<cstdio>
#include<bitset>
#define LL long long
using namespace std;

const int MAXN = 100000+9;
const int MX = 100000;

int pri[MAXN],tot,mu[MAXN],sum[MAXN];
bitset<MAXN> tag;

inline void Get_mu(){
mu[1] = 1;
for (int i=2;i<=MX;i++){
if (!tag[i]) pri[++tot] = i, mu[i] = -1;
for (int j=1;j<=tot && pri[j]*i<=MX;j++) {
tag[pri[j]*i] = 1;
if (i % pri[j]) mu[i*pri[j]] = -mu[i];
else {mu[i*pri[j]] = 0; break;}
}
}
for (int i=1;i<=MX;i++) sum[i] = sum[i-1] + mu[i];
}

int main(){
Get_mu(); int n,m,a,b,tmp;LL vout=0;
scanf("%d%d",&n,&m); if (n > m) swap(n,m);
for (int k=1;k<=n;k++) {
a = n/k; b = m/k;
for (int i=1;i<=a;i=tmp+1){
tmp = min(a/(a/i),b/(b/i));
vout += (LL)k*(sum[tmp]-sum[i-1])*(a/i)*(b/i);
}
}
printf("%lld\n",vout*2-(LL)n*m);
return 0;
}


## 【BZOJ 2301】[HAOI2011] Problem b

#include<iostream>
#include<cstdio>
#include<bitset>
using namespace std;

const int MAXN = 50000+9;
const int MX = 50000;

bitset<MAXN> tag;
int pri[MAXN],tot,mu[MAXN],sum[MAXN];

char c=getchar(); int buf=0,f=1;
while (c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
while (c<='9'&&c>='0'){buf=buf*10+c-'0';c=getchar();}
return buf*f;
}

inline void Get_mu(){
mu[1] = 1;
for (int i=2;i<=MX;i++) {
if (!tag[i]) pri[++tot] = i, mu[i] = -1;
for (int j=1;j<=tot && pri[j]*i<=MX;j++) {
tag[pri[j]*i] = 1;
if (i % pri[j]) mu[i*pri[j]] = -mu[i];
else {mu[i*pri[j]] = 0; break;}
}
}
for (int i=1;i<=MX;i++) sum[i] = sum[i-1] + mu[i];
}

int main(){
while (T--) {