莫比乌斯反演与容斥原理

退役前一直想把莫比乌斯反演与容斥原理统一在一起,奈何自己水平不足,只能作罢。
这次把《组合数学》、《具体数学》、《初等数论》的相关内容读了一遍,总算是完成了这个遗愿:
mobius_and_inclusion_exclusion_principle

Download:http://oi.cyo.ng/wp-content/uploads/2018/03/mobius_and_inclusion_exclusion_principle.pdf
拓展阅读1:http://blog.miskcoo.com/2015/12/inversion-magic-binomial-inversion
拓展阅读2:http://vfleaking.blog.uoj.ac/blog/87

【BZOJ 4635】数论小测验

相关链接

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=4635
数据生成器:http://paste.ubuntu.com/24360378/
神犇题解:http://www.cnblogs.com/clrs97/p/5625508.html

解题报告

对于第一问,显然可以用莫比乌斯反演做到$O(Tm)$

对于第二问,考虑枚举$k$,然后$10^3$以内的数最多有$4$种不同的质因数
于是搞一个状压$DP$,用矩阵快速幂优化
单词询问时间复杂度:$O(32m^2 + 32^3 log (n))$
看起来蛮有希望过的,但卡了一下午常也没有卡过 QwQ

正解的话,我们可以学习Claris用容斥
具体来讲枚举$k$,然后枚举$k$的每一个质因数是否满足条件
然后配合一点预处理之类的,就可以做到$O(m^2 + m \log m)$了

Code

这份代码在BZOJ被卡常了 QwQ

#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];
 
inline int read() {
    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() {
    int T = read(), type = read();
    if (type == 1) {
        GetMu(); 
        for (int t=1,vout;vout=0,t<=T;t++) {
            n = read(); m = read(); 
            int L = read(), R = read();
            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++) {
            n = read(); m = read(); l = read(), r = read();
            for (int w=l;w<=r;w++) vout = (vout + solve(w)) % MOD;
            printf("%d\n",vout);
        }
    }
    return 0;
}

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

相关链接

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=3994
数据生成器:http://paste.ubuntu.com/23992676/
神犇题解:https://blog.sengxian.com/solutions/bzoj-3994

解题报告

这题要用到一个结论: $d(xy) = \sum\limits_{i|x} {\sum\limits_{j|y} {[\gcd (i,j) = 1]} } $
这个结论似乎没有办法从意义上去推导,证明也只能是展开后发现刚好相等
但这个结论前不久集训的时候就考过一次,然而做这题的时候一点印象都没有
我要是下一次还记不住这个结论就直播吃键盘!

好了现在开始说正解
知道上面的结论以后,答案就成了 $\sum\limits_{k = 1}^n {\mu (k)\sum\limits_{i = 1}^{\frac{n}{k}} {d(i)\sum\limits_{j = 1}^{\frac{m}{k}} {d(j)} } } $
然后设$f(i)$为除数函数的前缀和,答案就变成了: $\sum\limits_{k = 1}^n {\mu (k)f(\frac{n}{k})f(\frac{m}{k})} $
于是直接下底函数分块就可以啦!

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;

inline int read() {
	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();
	for (int T=read();T;T--,vout=0) {
		n = read(); m = read();
		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] 选数

相关链接

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=3930
神犇题解:http://www.cnblogs.com/iwtwiioi/p/4986316.html

解题报告

这题,仔细想一想,岂不是跟这个题一样:http://oi.cyo.ng/?p=498
于是答案就是这个东西:$ \sum\limits_{i = 1}^m {\sum\limits_{j = 1}^m { \cdots \sum\limits_{k = 1}^m {[gcd(i,j, \cdots ,k) = K]} } } = \sum\limits_{i = 1}^{\left\lfloor {\frac{m}{K}} \right\rfloor } {\sum\limits_{j = 1}^{\left\lfloor {\frac{m}{K}} \right\rfloor } { \cdots \sum\limits_{k = 1}^{\left\lfloor {\frac{m}{K}} \right\rfloor } {[gcd(i,j, \cdots ,k) = 1]} } } = \sum\limits_{d = 1}^{\left\lfloor {\frac{m}{K}} \right\rfloor } {\mu (d) \cdot {{\left\lfloor {\frac{m}{{Kd}}} \right\rfloor }^n}}$
于是剩下的锅就是预处理莫比乌斯函数的前缀和了!
然而我只会 $ O(n)$ 的做法,于是就跪了

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)}$
换一句话说,莫比乌斯函数的前缀和可以做到 $ O(\sqrt n \log n)$
这样的话,感觉复杂度是 $ O(n \log n)$ 的
然而PoPoQQQ大爷说,预处理前 $ S$ 个值的话,复杂度就可以变成:$ O(\frac{n}{S}\sqrt n lo{g_2}\frac{n}{S})$
然而并不知道怎么证明,感觉很有道理的样子!

当然这题还有容斥的做法
首先有结论:$ gcd(i,j) \le (r – l),(i,j \in (l,r))$
证明的话,参见这里:https://blog.sengxian.com/solutions/bzoj-3930
于是考虑 $ f(i)$ 表示gcd为i的方案数,则 $ f(i) = {\left\lfloor {\frac{m}{{Ki}}} \right\rfloor ^n} – \sum\limits_{i|j} {f(j)}$
这样的话, $ f(1)$ 就是答案辣!

—————————— UPD 2017.2.20 ——————————
突然发现,莫比乌斯函数的前缀和不是杜教筛的模板题吗?
当然PoPoQQQ的做法也有可取之处:不用复杂度公式推导

【BZOJ 3025】[Balkan2003] Farey数列

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=3025
数据加强版:http://oi.cdshishi.net/contestoi/problem_show.php?pid=1063
数据生成器(要求支持假分数):http://paste.ubuntu.com/23391949/

好久没有做推式子的题目啦!
首先这题可以用Farey sequence直接做:https://en.wikipedia.org/wiki/Farey_sequence
具体实现上可能需要用Stern–Brocot tree:https://en.wikipedia.org/wiki/Stern%E2%80%93Brocot_tree
上述做法把此题变成裸题,然而不知道上面的东西怎么办?
下面说一说使用莫比乌斯函数的O(nlog(n^2))的做法
注:以下讨论均针对于数据加强版,BZOJ上的原题可能不需要这么复杂

最外层先搞一个二分,假设二分到某一个值使得小于等于他的分数刚好有k个
那么我们只需要枚举分母,找到其中最大的一个输出即可
另外,这里的二分只能定分母、二分分子。不能够直接二分小数,原因待会儿说。
考虑到精度问题,我们定分母EPS=1e10(实际上也只能定为1e10,因为小了精度不够,大了要炸long long)

假设二分的分子为k,则现在我们需要求解的便是:\(\sum\limits_{y = 1}^n {\sum\limits_{x = 1}^{\left\lfloor {\frac{{yk}}{{EPS}}} \right\rfloor } {\sum\limits_{d|x,y} {\mu (d)} } } \)
稍微变形一下可以得到:\(\sum\limits_{d = 1}^n {\mu (d)\sum\limits_{y = 1}^{\left\lfloor {\frac{n}{d}} \right\rfloor } {\sum\limits_{x = 1}^{\left\lfloor {\frac{{\left\lfloor {\frac{{y \cdot d \cdot k}}{{EPS}}} \right\rfloor }}{d}} \right\rfloor } 1 } } \)
注意到x的上限有两个取整函数,因为取整函数有这个性质:\(\left\lfloor {\frac{{\left\lfloor {\frac{{\rm{y}}}{a}} \right\rfloor }}{b}} \right\rfloor = \left\lfloor {\frac{y}{{ab}}} \right\rfloor\)
所以才可以去掉上面的那个取整函数,这也是为什么不能直接二分实数,因为实数去不掉那个取整函数

我们继续推导,再次化简以后可以得到:\(\sum\limits_{d = 1}^n {\mu (d)\sum\limits_{y = 1}^{\left\lfloor {\frac{n}{d}} \right\rfloor } {\sum\limits_{x = 1}^{\left\lfloor {\frac{{yk}}{{EPS}}} \right\rfloor } 1 } } \)
这样的话,就可以线性处理后两个Σ,然后O(n)枚举d来计算了
另外还有一个小细节:x不仅受到分数值的限制,还受到n的限制,所以那里还需要特判一下

这题还有一种类似的做法,但可避免炸long long的风险:
考虑仪仗队那个题,因为方队左右对称,于是总可以把问题划归到方阵的右下部分
然后二分右下部分的斜率(定分母为n,二分分子)
二分到分子差只有1的时候,就暴力取出那个区间的所有分数,排序输出

这题我写的代码是第一种方式的,精度巨坑,必须是1e10,且不能有任何实数运算

#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];

inline int read(){
	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

题目传送门:http://www.spoj.com/problems/VLATTICE/
离线版题目:http://paste.ubuntu.com/20300942/

仪仗队升级版!
然而还是不会QAQ 还是可耻地看了题解QAQ

#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

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=2820
离线版题目:http://paste.ubuntu.com/20287584/

这个东西嘛,如果没有多组询问的话,枚举质数就可以水过去
但有多组询问QAQ
前排膜拜hzwer:http://hzwer.com/6142.html
现在身心俱疲,什么都不想说了 累…..

#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的问题

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=2705
离线版题目:http://paste.ubuntu.com/20283268/

这题真的是蜜汁复杂度。
自己想的时候想到了,然而没敢写QAQ

首先是枚举gcd,然后搞phi或者mu
这个很显然,但只能过60%的数据

因为这个题是求gcd(i,n),所以gcd的取值只会有σ(n)个(有一个是定值),也就是n的因数的个数
然后用sqrt(n)来求每个因数的phi()

这样的话,时间复杂度上界是(sqrt(n))^2=O(n)的
这题n=10^9那不T到死QAQ
然而万能的vfk告诉我们,10^9内σ(n)最大为1000的样子,这样就不会T了QAQ
前排膜拜vfk:http://vfleaking.blog.163.com/blog/static/174807634201341913040467/
前排膜拜hht:http://techotaku.lofter.com/post/4856f0_634219b
ps:hht告诉我们,除数函数可以线筛,然而一脸懵逼QAQ

然而这题的还有一个trick,求phi()那里,还可以dfs算QAQ
复杂度更优!快4倍的样子
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 4802
也是求$\phi (n)$但$n \le 10^{18}$

【BZOJ 2190】[SDOI2008] 仪仗队

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=2190
离线版题目:http://paste.ubuntu.com/20271981/

这个是“能量采集”的弱化版,所以mu和phi都可以很简单做
线筛写错了,wa了好久QAQ

#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;
} 

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

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=2005
前传传送门:http://oi.cyo.ng/?p=477

之前%JCVB的没有成功。
今天又来%,基本上有思路了,真的是_(:з」∠)_

引理:\(\varphi (n) = \sum\limits_{d|n} {\frac{n}{d} \cdot \mu (d)}\)
证明:\(\varphi (n) = \sum\limits_{i = 1}^n {[\gcd (i,n) = = 1] = \sum\limits_{i = 1}^n {\sum\limits_{d|i} {\sum\limits_{d|n} {\mu (d) = \sum\limits_{d|n} {\mu (d) \cdot \sum\limits_{i = 1}^n {[\gcd (i,n) = = i] = } } } } } } \sum\limits_{d|n} {\mu (d) \cdot \frac{n}{d}}\)

我们之前的式子是:\(\sum\limits_{d = 1}^n {\sum\limits_{k = 1}^{\left\lfloor {\frac{n}{d}} \right\rfloor } {d \cdot \mu (k) \cdot \left\lfloor {\frac{n}{{d \cdot k}}} \right\rfloor } } \cdot \left\lfloor {\frac{m}{{d \cdot k}}} \right\rfloor\)
注意看好,我要变形啦!( •̀ ω •́ )y
不妨设:D=k*d,则有原式=\(\sum\limits_{D = 1}^n {\sum\limits_{d|D} {d \cdot \mu (\frac{D}{d}) \cdot } } \left\lfloor {\frac{n}{D}} \right\rfloor \cdot \left\lfloor {\frac{m}{D}} \right\rfloor \)
根据引理可以进一步简化为:\(\sum\limits_{D = 1}^n {\varphi ({\rm{D}})} \cdot \left\lfloor {\frac{n}{D}} \right\rfloor \cdot \left\lfloor {\frac{m}{D}} \right\rfloor \)
然后就可以线筛出欧拉函数后,\(\sqrt n \)分块就好,总时间复杂度:O(n)

#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)
为什么会有这个结论呢?
因为这个相当于把线段分成了完全相同的gcd(x,y)段
或者可以这样理解:
离原点最近的一个肯定是(x/gcd(x,y),y/gcd(x,y))
之后的每一整点都是他的整数倍,所以总共有gcd(x,y)个

—————————— UPD 2017.2.1 ——————————
今天又看了看这个玩意儿,为什么要$\sqrt n$分块QwQ
暴力不也这复杂度吗……

而且完全不知道$\sum\limits_{d|n} {\frac{n}{d} \cdot \mu (d)}$是$\varphi (n)$也完全可以做啊!
根据狄利克雷卷积,可以知道这货是积性函数,然后再推一推发现可以线筛,然后直接做就好啊!

—————————— UPD 2017.7.3 ——————————
我发现自己之前数学好强
现在看这个题已经不会了qwq

【BZOJ 2818】Gcd

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=2818
离线版题目:http://paste.ubuntu.com/20184097/

这个题目,很容易想到枚举素数,然后求gcd=1的对数
用莫比乌斯函数+分段来搞的话,好像很常规,时间复杂度O(ln(n)*sqrt(n)+n)
然后就是看到这篇博客,发现可以直接用phi函数搞QAQ,时间复杂度O(ln(n)+n)
因为phi前缀和那里必须用long long搞得实际运行时间比莫比乌斯函数还要慢QAQ
莫比乌斯函数版本:

#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] 能量采集

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=2005
离线版题目:http://paste.ubuntu.com/20171175/

这题,O(n^0.5+n)的做法我不会QAQ
%JCVB,看了半小时,还是不懂QAQ
只会O(n^1.5+n)的做法:
枚举gcd(i,j)的值,然后就和“Zap”一样了
居然没有T,好感动இ௰இ

#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

题目传送门:http://www.lydsy.com/JudgeOnline/problem.php?id=2301
离线版题目:http://paste.ubuntu.com/20155288/

这个题目和1101简直一毛一样QAQ
就是要减一减。貌似黄学长的做法是容斥一小下。
然而我比较笨,不会容斥。于是就只能暴力搞一搞。
结果一言不合就 Rank 3 QAQ
mobimus_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];

inline int read(){
	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(){
	Get_mu(); int T = read();
	while (T--) {
		int a=read()-1,b=read(),c=read()-1,d=read(),k=read(),vout;
		a /= k; b /= k; c /= k; d /= k; vout = 0;
		for (int i=1,lim=min(b,d),tmp;i<=lim;i=tmp+1){
			tmp = min(b/(b/i),d/(d/i));
			if (a/i) tmp = min(tmp, a/(a/i));
			if (c/i) tmp = min(tmp, c/(c/i));
			vout += (sum[tmp]-sum[i-1])*(b/i-a/i)*(d/i-c/i);
		}
		printf("%d\n",vout);
	}
	return 0;
}