[关闭]
@KirinBill 2017-10-12T04:44:29.000000Z 字数 11726 阅读 525

莫比乌斯反演学习笔记

学习笔记 数论

感觉这个东西很需要感觉,要多练习。。。

目录


基础知识

莫比乌斯函数

定义整数域上的数论函数:


这个函数就是莫比乌斯函数

  1. 性质

    • 显然,为积性函数,但不是完全积性的,不过我们依然可以用欧拉筛法求出;
    • ,所以的狄利克雷前缀和就是元函数。将看作容斥系数的话就很好理解了;
    • ,设反演一下就出来了,学完下面的就能证了。
  2. 练习:[BZOJ 2440] 完全平方数

    • 其实是用到作为容斥系数的应用,这个和后面的反演有些类似;
    • 考虑二分答案,则我们只需快速计算~内无平方因子的数的个数即可;
    • 显然算出有平方因子的数的个数更好算,用容斥原理则可知:
    • 观察发现,枚举的这个前面的系数刚好和值一样,那么,因为如果中有重复质因子,刚好不影响,反正就很神奇;
    • 复杂度为

    代码:

  1. #include <cstdio>
  2. #include <cmath>
  3. using std::sqrt;
  4. using std::floor;
  5. const int INF=2e9,MAXN=50005;
  6. int miu[MAXN];
  7. long long sqr[MAXN];
  8. inline void Euler_miu(){
  9. static int prm[MAXN];
  10. static bool notP[MAXN];
  11. notP[1]=true,miu[1]=1;
  12. for(int i=2;i<MAXN;++i){
  13. if(!notP[i]){
  14. prm[++prm[0]]=i;
  15. miu[i]=-1;
  16. }
  17. for(int j=1,tmp;j<=prm[0] && i*prm[j]<MAXN;++j){
  18. tmp=i*prm[j];
  19. notP[tmp]=true;
  20. if(i%prm[j]==0){
  21. miu[tmp]=0;
  22. break;
  23. }
  24. else miu[tmp]=-miu[i];
  25. }
  26. }
  27. }
  28. inline void pre_tab(){
  29. for(int i=1;i<MAXN;++i)
  30. sqr[i]=(long long)i*i;
  31. }
  32. inline int cal(int r){
  33. long long ret=0;
  34. for(int i=1,lim=floor(sqrt(r));i<=lim;++i)
  35. ret+=r/sqr[i]*miu[i];
  36. return ret;
  37. }
  38. inline int bin_chop(int k){
  39. int l=1,r=INF;
  40. long long mid;
  41. while(l<=r){
  42. mid=(long long)l+r>>1;
  43. if(cal(mid)>=k) r=mid-1;
  44. else l=mid+1;
  45. }
  46. return l;
  47. }
  48. int main(){
  49. Euler_miu();
  50. pre_tab();
  51. int T;
  52. scanf("%d",&T);
  53. int k;
  54. while(T--){
  55. scanf("%d",&k);
  56. printf("%d\n",bin_chop(k));
  57. }
  58. return 0;
  59. }

公式与证明

对于两个函数,若狄利克雷前缀和,即:,则有如下结论:

  1. 形式一:

    • 这个也就是说狄利克雷卷积
    • 证明:
  2. 形式二:

    • 这个形式有些“奇怪”,但却是最常用的一种;
    • 证明(令):
  3. 应用方法

    • 即使知道了公式,我们发现这玩意似乎并无卵用,其实不然;
    • 如果我们要求,但是不容易迅速得出,可以考虑设,若很好求,那么反演一下就出来了;
    • 其实这个反演过程就是一种容斥,加加减减就神奇的算出来了;
    • 在应用过程中通常是的形式,为了迅速反演求出,常常利用数论分段求和的技巧;
    • 以上这些见到下面的题就清楚了,这个东西多练才能有感觉。

练习

[HAOI 2011] Problem b

思路

代码

  1. #include <cstdio>
  2. #include <algorithm>
  3. using std::min;
  4. const int MAXN=50005;
  5. int miu[MAXN],sum[MAXN];
  6. inline void Euler_miu(){
  7. static int prm[MAXN];
  8. static bool notP[MAXN];
  9. notP[1]=true,miu[1]=1;
  10. for(int i=2;i<MAXN;++i){
  11. if(!notP[i]){
  12. prm[++prm[0]]=i;
  13. miu[i]=-1;
  14. }
  15. for(int j=1,tmp;j<=prm[0] && i*prm[j]<MAXN;++j){
  16. tmp=i*prm[j];
  17. notP[tmp]=true;
  18. if(i%prm[j]==0){
  19. miu[tmp]=0;
  20. break;
  21. }
  22. else miu[tmp]=-miu[i];
  23. }
  24. }
  25. }
  26. inline void pre_sum(){
  27. Euler_miu();
  28. for(int i=1;i<MAXN;++i)
  29. sum[i]=sum[i-1]+miu[i];
  30. }
  31. inline int cal(int n,int m){
  32. int ret=0;
  33. for(int i=1,lim=min(n,m),r;i<=lim;i=r+1){
  34. r=min(n/(n/i),m/(m/i));
  35. if(r>lim) r=lim;
  36. ret+=(sum[r]-sum[i-1])*(n/i)*(m/i);
  37. }
  38. return ret;
  39. }
  40. inline int solve(int a,int b,int c,int d,int k){
  41. --a,--c;
  42. a/=k,b/=k,c/=k,d/=k;
  43. return cal(b,d)-cal(b,c)-cal(a,d)+cal(a,c);
  44. }
  45. int main(){
  46. pre_sum();
  47. int n;
  48. scanf("%d",&n);
  49. int a,b,c,d,k;
  50. while(n--){
  51. scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
  52. printf("%d\n",solve(a,b,c,d,k));
  53. }
  54. return 0;
  55. }

[BZOJ 2818] Gcd

思路

代码

  1. #include <cstdio>
  2. const int MAXN=1e7+5;
  3. int prm[MAXN],miu[MAXN];
  4. inline void Euler_miu(int n){
  5. static bool notP[MAXN];
  6. notP[1]=true,miu[1]=1;
  7. for(int i=2;i<=n;++i){
  8. if(!notP[i]){
  9. prm[++prm[0]]=i;
  10. miu[i]=-1;
  11. }
  12. for(int j=1,tmp;j<=prm[0] && i*prm[j]<=n;++j){
  13. tmp=i*prm[j];
  14. notP[tmp]=true;
  15. if(i%prm[j]==0){
  16. miu[tmp]=0;
  17. break;
  18. }
  19. else miu[tmp]=-miu[i];
  20. }
  21. }
  22. }
  23. int main(){
  24. int n;
  25. scanf("%d",&n);
  26. Euler_miu(n);
  27. long long ans=0;
  28. for(int i=1;i<=prm[0];++i){
  29. for(int j=prm[i],k=1,tmp;j<=n;j+=prm[i],++k){
  30. tmp=n/j;
  31. ans+=(long long)miu[k]*tmp*tmp;
  32. }
  33. }
  34. printf("%lld",ans);
  35. return 0;
  36. }

[SDOI 2014] 数表

思路

代码

  1. #include <cstdio>
  2. #include <algorithm>
  3. #include <cstring>
  4. #include <climits>
  5. using std::min;
  6. using std::sort;
  7. const int MAXQ=2e4+5,MAXN=1e5+5;
  8. int miu[MAXN];
  9. struct data{
  10. int id;
  11. long long c;
  12. friend bool operator< (const data &a,const data &b){
  13. return a.c<b.c;
  14. }
  15. }sgm[MAXN];
  16. struct query{
  17. int id,n,m,a,ans;
  18. static bool cmp_a(const query &a,const query &b){
  19. return a.a<b.a;
  20. }
  21. static bool cmp_id(const query &a,const query &b){
  22. return a.id<b.id;
  23. }
  24. }qry[MAXQ];
  25. class BIT{
  26. private:
  27. unsigned c[MAXN];
  28. int lowbit(int x){return x&-x;}
  29. unsigned qry(int r){
  30. unsigned ret=0;
  31. for(;r;r-=lowbit(r))
  32. ret+=c[r];
  33. return ret;
  34. }
  35. public:
  36. void add(int l,unsigned x){
  37. for(;l<MAXN;l+=lowbit(l))
  38. c[l]+=x;
  39. }
  40. unsigned qry(int l,int r){
  41. return qry(r)-qry(l-1);
  42. }
  43. }ta;
  44. inline void Euler_sv(){
  45. static int prm[MAXN];
  46. static long long dProd[MAXN],prodS[MAXN];
  47. static bool notP[MAXN];
  48. notP[1]=true,miu[1]=sgm[1].c=1;
  49. for(int i=2;i<MAXN;++i){
  50. if(!notP[i]){
  51. prm[++prm[0]]=i;
  52. dProd[i]=i,prodS[i]=1+i;
  53. miu[i]=-1,sgm[i].c=1+i;
  54. }
  55. for(int j=1,tmp;j<=prm[0] && i*prm[j]<MAXN;++j){
  56. tmp=i*prm[j];
  57. notP[tmp]=true;
  58. if(i%prm[j]==0){
  59. dProd[tmp]=dProd[i]*prm[j];
  60. prodS[tmp]=prodS[i]+dProd[tmp];
  61. miu[tmp]=0,sgm[tmp].c=sgm[i].c/prodS[i]*prodS[tmp];
  62. break;
  63. }
  64. else{
  65. dProd[tmp]=prm[j];
  66. prodS[tmp]=1+prm[j];
  67. miu[tmp]=-miu[i],sgm[tmp].c=sgm[i].c*sgm[prm[j]].c;
  68. }
  69. }
  70. }
  71. }
  72. inline unsigned cal(int n,int m){
  73. unsigned ret=0;
  74. for(int l=1,r,lim=min(n,m);l<=lim;l=r+1){
  75. r=min(n/(n/l),m/(m/l));
  76. if(r>lim) r=lim;
  77. ret+=(unsigned)(n/l)*(m/l)*ta.qry(l,r);
  78. }
  79. return ret;
  80. }
  81. int main(){
  82. Euler_sv();
  83. for(int i=1;i<MAXN;++i)
  84. sgm[i].id=i;
  85. sort(sgm+1,sgm+MAXN);
  86. int q;
  87. scanf("%d",&q);
  88. for(int i=1;i<=q;++i){
  89. scanf("%d%d%d",&qry[i].n,&qry[i].m,&qry[i].a);
  90. qry[i].id=i;
  91. }
  92. sort(qry+1,qry+q+1,query::cmp_a);
  93. for(int i=1,j,k=1;i<=q;++i){
  94. for(;sgm[k].c<=qry[i].a;++k){
  95. for(j=1;j*sgm[k].id<MAXN;++j)
  96. ta.add(j*sgm[k].id,(unsigned)miu[j]*sgm[k].c);
  97. }
  98. qry[i].ans=cal(qry[i].n,qry[i].m)&INT_MAX;
  99. }
  100. sort(qry+1,qry+q+1,query::cmp_id);
  101. for(int i=1;i<=q;++i)
  102. printf("%d\n",qry[i].ans);
  103. return 0;
  104. }

[BZOJ 2154] Crash的数字表格

思路

代码

  1. #include <cstdio>
  2. #include <algorithm>
  3. using std::min;
  4. const int MAXN=1e7+5,P=20101009;
  5. int inv_6;
  6. int miu[MAXN],tab[MAXN];
  7. inline void Euler_miu(){
  8. static int prm[MAXN];
  9. static bool notP[MAXN];
  10. notP[1]=true,miu[1]=1;
  11. for(int i=2;i<MAXN;++i){
  12. if(!notP[i]){
  13. prm[++prm[0]]=i;
  14. miu[i]=-1;
  15. }
  16. for(int j=1,tmp;tmp=i*prm[j],j<=prm[0] && tmp<MAXN;++j){
  17. notP[tmp]=true;
  18. if(i%prm[j]==0){
  19. miu[tmp]=0;
  20. break;
  21. }
  22. else miu[tmp]=-miu[i];
  23. }
  24. }
  25. }
  26. inline void pre_tab(){
  27. Euler_miu();
  28. for(int i=1;i<MAXN;++i)
  29. tab[i]=(tab[i-1]+(long long)i*i*miu[i]%P+P)%P;
  30. for(int i=2;i<MAXN;++i)
  31. miu[i]+=miu[i-1];
  32. }
  33. inline int sum(int n,int m){
  34. return (long long)((long long)n*(n+1)>>1)%P*(((long long)m*(m+1)>>1)%P)%P;
  35. }
  36. inline int cal(int n,int m){
  37. int ret=0;
  38. for(int l=1,r,lim=min(n,m);l<=lim;l=r+1){
  39. r=min(n/(n/l),m/(m/l));
  40. if(r>lim) r=lim;
  41. ret=(ret+(long long)(tab[r]-tab[l-1]+P)%P*sum(n/l,m/l)%P)%P;
  42. }
  43. return ret;
  44. }
  45. int main(){
  46. pre_tab();
  47. int n,m;
  48. scanf("%d%d",&n,&m);
  49. int ans=0;
  50. for(int l=1,r,lim=min(n,m);l<=lim;l=r+1){
  51. r=min(n/(n/l),m/(m/l));
  52. if(r>lim) r=lim;
  53. ans=(ans+(long long)cal(n/l,m/l)*(((long long)(l+r)*(r-l+1)>>1)%P)%P)%P;
  54. }
  55. printf("%d",(ans+P)%P);
  56. return 0;
  57. }

[SDOI 2015] 约数个数和

思路

代码

  1. #include <cstdio>
  2. #include <algorithm>
  3. using std::min;
  4. const int MAXN=50005;
  5. int miu[MAXN];
  6. long long tao[MAXN];
  7. inline void Euler_sv(){
  8. static int prm[MAXN],dExp[MAXN];
  9. static bool notP[MAXN];
  10. notP[1]=true,dExp[1]=miu[1]=tao[1]=1;
  11. for(int i=2;i<MAXN;++i){
  12. if(!notP[i]){
  13. prm[++prm[0]]=i;
  14. dExp[i]=1;
  15. miu[i]=-1,tao[i]=2;
  16. }
  17. for(int j=1,tmp;j<=prm[0] && i*prm[j]<MAXN;++j){
  18. tmp=i*prm[j];
  19. notP[tmp]=true;
  20. if(i%prm[j]==0){
  21. dExp[tmp]=dExp[i]+1;
  22. miu[tmp]=0;
  23. tao[tmp]=tao[i]/(dExp[i]+1)*(dExp[tmp]+1);
  24. break;
  25. }
  26. else{
  27. dExp[tmp]=1;
  28. miu[tmp]=-miu[i];
  29. tao[tmp]=tao[i]*tao[prm[j]];
  30. }
  31. }
  32. }
  33. }
  34. inline void pre_tab(){
  35. Euler_sv();
  36. for(int i=2;i<MAXN;++i){
  37. miu[i]+=miu[i-1];
  38. tao[i]+=tao[i-1];
  39. }
  40. }
  41. inline long long cal(int n,int m){
  42. long long ret=0;
  43. for(int l=1,r,lim=min(n,m);l<=lim;l=r+1){
  44. r=min(n/(n/l),m/(m/l));
  45. if(r>lim) r=lim;
  46. ret+=(long long)tao[n/l]*tao[m/l]*(miu[r]-miu[l-1]);
  47. }
  48. return ret;
  49. }
  50. int main(){
  51. pre_tab();
  52. int T;
  53. scanf("%d",&T);
  54. int n,m;
  55. while(T--){
  56. scanf("%d%d",&n,&m);
  57. printf("%lld\n",cal(n,m));
  58. }
  59. return 0;
  60. }
添加新批注
在作者公开此批注前,只有你和作者可见。
回复批注