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