先拜大牛。感谢贾志鹏严谨的思维。以及简单清晰的论文描述。

一定要结合论文看。我只是提出我觉得关键的部分。论文在网上随处可见。贾志鹏线性筛。

开头两种线性筛的比较。

一种是传统的线性筛。时间复杂度为N*log(log(N))。

另外一种是优化了合数的筛法。文中称作Euler线性筛。

其优化的地方。

举个例子:合数6。 是2的倍数也是3的倍数。 当你用传统的筛法的时候在遍历2的倍数的时候会遍历到6。遍历3的倍数的时候同样也会遍历到6。

而另外一种只会筛出6为2的倍数。3就不会筛6了。

另外个人认为筛法二有一个很重要的思想。当i为合数的时候。其实脑海里不认为是合数。而是素数的乘积。这样就能比较直观地确定这个算法的正确性了。

积性函数。

分为完全积性和条件积性。

我们最喜欢的积性。大概就是互素积性了。因为满足互素积性的话。根据算术基本定理。就能够简单做到推广到任意实数。

f(1) = 1 。 这个在我们高中数学题。抽象函数。就已经能简单知道了。

欧拉函数。就不再谈了。包括其线性筛的那一步至关重要的证明。也在我其它博文提到过了。

其 欧拉定理和费马小定理的作用。我得再多补充一点。

以及互质数和。 n的互质数和为 n*φ(n)/2.

莫比乌斯函数和容斥定理的关系。

可以发现莫比乌斯函数其实就是容斥定理的映射一般。

莫比乌斯函数 是我们再熟悉不过的了。不熟悉可以看这里。

首先看 (-1)^r m = p1p2p3p4p5pr 其实就是在模拟容斥定理。

假如一但不是素数。那就为0.

两个函数的线性筛。这其实是我们处理问题的基本。这里需要讲的是。不一定只有积性函数才可以用这种筛法。

只要你能找到f(kn) n整除k 和不整除的两个时刻所对应的递推式。这个在扩展问题中会出现。

问题一:求1~N对质数P的乘法逆元。

关于f(n)为完全积性函数。根据同余定理可以简单获得。要证明的话。减法证同余即可。

P = nt + k

n'  ≡ n*(t^2)*f(k)^2 (mod P)

这个证明过程很漂亮(很佩服这么顺畅,思维这么清晰)。也是根据同余定理。还有逆元的性质。就能推理的。

这个问题的意义。可以求N!的 mod P 的逆元了。逆元还是很有用的。因为毕竟除法并没有特别好的同余式。(依稀还记得那两个。)

问题二:给T组N,M.依次求出的值.(N,M<=10^6,T<=10^3)

求解gcd(a,b).把gcd(a,b)当做n.再通过欧拉函数和式。推导过程如下。

第二个等式是用d来看待式子的方法来化简和式的。

之后再穷举d即可。

#include<stdio.h>
#include<string.h>
#define    N 100bool mark[N+5];
int prime[N+5];
int num;
int euler[N+5];
int Min(int a,int b){return a>b?a:b;}
void Euler()
{int i,j;num = 0;memset(euler,0,sizeof(euler));memset(prime,0,sizeof(prime));memset(mark,0,sizeof(mark));euler[1] = 1; // multiply functionfor(i=2;i<=N;i++){if(!mark[i]){    prime[num++] = i;euler[i] = i-1;}for(j=0;j<num;j++){if(i*prime[j]>N){break;}mark[i*prime[j]] = 1;if(i%prime[j]==0){euler[i*prime[j]] = euler[i]*prime[j];break;}else{euler[i*prime[j]] = euler[prime[j]]*euler[i];}}}
}int main()
{int i;int M1,M2;Euler();for(i=0;i<num;i++)printf("%d ",prime[i]);printf("\n");for(i=1;i<=N;i++)printf("%d ",euler[i]);printf("\n");M1 = 2;M2 = 3;int sum = 0;int min = Min(M1,M2);for(i=1;i<=min;i++){sum += euler[i]*(M1/i)*(M2/i);}printf("%d\n",sum);
}

the second problem test

问题三:给T组N,M.依次求出的值.(N,M<=10^6,T<=10^3)

在证明之前,先证明以下式子。

问题的解决推导。

第一个等式。lcm(a,b) = ab/gcd(a,b).

第二个等式。令d=gcd(a,b)。

第三个等式。转化为d的视角。(这个手法经常有)。

第四个等式。转化为莫比乌斯函数。

第五个等式。利用上述的等式来转化。注意d和d'

第六个等式。论文中提到的有趣的化简性质。

第七个等式。其实是d = dd'换元。不过有点老奸巨猾啊。干嘛不设个T = dd'。这个我纠结了半天。

之后就是如论文中介绍的。g(d) 为积性函数。线性筛之。

总体上算法还是N的。

问题四:给T组N,M.依次求出的值.(N,M<=10^6,T<=10^3)

实质上就是求 其中x和y互素。的对数。

我们是时候需要有和式化成的思想了。[gcd(a,b)=1]真是漂亮的莫比乌斯函数的和式的结果。

第一个等式:莫比乌斯函数扩写

第二个等式:gcd(a,b)=p -> gcd(a/p,b/p)=1问题转换。

第三个等式:一个和式的处理手段。

第四个等式:很常见的。

#include<stdio.h>
#include<string.h>
#define    N 100bool mark[N+5];
int prime[N+5];
int num;
int mobi[N+5];
int Min(int a,int b){return a>b?a:b;}
void Mobi()
{int i,j;num = 0;memset(mobi,0,sizeof(mobi));memset(prime,0,sizeof(prime));memset(mark,0,sizeof(mark));mobi[1] = 1; // multiply functionfor(i=2;i<=N;i++){if(!mark[i]){    prime[num++] = i;mobi[i] = -1;}for(j=0;j<num;j++){if(i*prime[j]>N){break;}mark[i*prime[j]] = 1;if(i%prime[j]==0){mobi[i*prime[j]] = 0;break;}else{mobi[i*prime[j]] = mobi[prime[j]]*mobi[i];}}}
}int main()
{int i;int M1,M2;Mobi();for(i=0;i<num;i++)printf("%d ",prime[i]);printf("\n");for(i=1;i<=N;i++)printf("%d ",mobi[i]);printf("\n");M1 = 2;M2 = 3;int sum = 0;int min = Min(M1,M2);for(i=1;i<=min;i++){sum += mobi[i]*(M1/i)*(M2/i);}printf("%d\n",sum);
}

Test problem forth

问题五:给T组N.依次求出的值.(N<=10^6,T<=10^3)

 

其实根据问题三.可以直接获得该化简出来的式子的。

然后解法和问题三一样。

但是论文上寻找积性f(n)直接筛出答案。

首先佩服一下其思维的紧密。一个变量啊。就寻找积性函数。这个转化真是清晰而又巧。

画个图就能知道 -n 是用来去重复的统计的。

.

f(n)是积性的。具体证明如论文上解释。

第一个等式:d = gcd(n,i)

第二个等式:k = i/d.且全部都除以d.gcd(a,b)=d转化成求互素(gcd(a,b)=1)的问题。

第三个等式:令d=n/d。是对应的。  其实在第二个等式就能看出是欧拉函数求约数和问题了。

第四个等式:不解释了吧。

第五个等式:手算一下容易得。

欧拉函数求互质数和的函数是积性函数

有一道题。就是利用这个。后会介绍。

见到积性函数我们现在应该是very happy的。

扩展问题1: T组N.依次求出的值.(N<=10^6,T<=10^3)

借鉴了贾志鹏上面所有问题的证明。这个是我自己写的扩展证明。难免有错误。见谅。还望留言提醒。

我觉得这样的证明是非常轻快明了的。然后网上还有流行一种。用莫比乌斯反演的另外一种表示式的。也是非常神奇的。

不过。那个反演我还没有证明过。不过还是借鉴了其下半部分的设T。(也是这个设T点醒了我。贾志鹏第3个问题的证明的最后一步。)

这里并不能因为p是素数。而n/p不一定是素数。所以并不是对称的。(如果看过具体数学就能很快明白了。)

处理

分类对于 g(kx) .有

 g(kx)=μ(x)                 k|p

 g(kx)=μ(x)-g(x)          k!|p

结合莫比乌斯函数。可以知道分类成立:

我们可以借这个 并且借用之前两个积性函数的筛法 来筛 g(n)。

这是明显可行的。也就是说。我们不需要函数必须是积性的才能去筛。

我们只需要找到g(kx)是由g(x)获得的。或者是在g(x)之前就筛掉的值获得的。eg:g(x-1) (筛法总是从小到大。)

更甚的是。我们只需要获得大值和小值的关系!就可以筛法。但是该筛法。是建立在素数表达式之上的。

这段阐述也许很混乱。但是我也只能描述个大概的个人体会。理解不理解没关系。

给你个例子。筛一个数的素因子之和。

对于上述的筛法.

让F[n] 为n的素因子之和。

F[i*prime[j]] = F[i]+prime[j]                            i\prime[j]时。

F[i*prime[j]] = F[i]+prime[j]                             i!\prime[j]时。 

两种情况是一样的。原因显而易见。不过我们还是得判断。因为i\prime[j]的时候我们更新后可以直接break;

再考虑一个问题:筛一个数的所拥有的素因子之和。 比如12 为2*2*3 我们只计算2+3.

那么有:

F[i*prime[j]] = F[i]                                             i\prime[j]时。

F[i*prime[j]] = F[i]+prime[j]                             i!\prime[j]时。 

对于这个问题的code.

#include<stdio.h>
#include<string.h>
#define N 100int num,prime[N+5],f[N+5];
bool mark[N+5];
void Init()
{int i,j;num = 0;f[1] = 0;for(i=2;i<=N;i++){if(!mark[i]){prime[num++] = i;f[i] = i;}for(j=0;j<num;j++){if(i*prime[j]>N){break;}mark[i*prime[j]] = 1;if(i%prime[j]==0){f[i*prime[j]] = f[i];break;}else{f[i*prime[j]] = f[i]+prime[j];    }}}
}int main()
{int i;Init();for(i=1;i<=N;i++){printf("%d = %d \n",i,f[i]);}
}

Problem test

再考虑一个问题:筛一个数的所拥有的不重复的素因子之和。比如12 为2*2*3 我们只计算3

那么有:

 i\prime[j]时。

  情况1: (i/prime[j])\prime[j]时.

  F[i*prime[j]] = F[i]

  情况2: (i/prime[j])!\prime[j]时.

  FF[i*prime[j]] = F[i]-prime[j].

 i!\prime[j]时。 

F[i*prime[j]] = F[i]+prime[j]                           

#include<stdio.h>
#include<string.h>
#define N 100int num,prime[N+5],f[N+5];
bool mark[N+5];
int Max(int a,int b)
{return a>b?a:b;
}
void Init()
{int i,j;num = 0;f[1] = 0;for(i=2;i<=N;i++){if(!mark[i]){prime[num++] = i;f[i] = i;}for(j=0;j<num;j++){if(i*prime[j]>N){break;}mark[i*prime[j]] = 1;if(i%prime[j]==0){if((i/prime[j])%prime[j]==0){f[i*prime[j]] = f[i];}else{f[i*prime[j]] = f[i] - prime[j];}break;}else{f[i*prime[j]] = f[i]+prime[j];    }}}
}int main()
{int i;Init();for(i=1;i<=N;i++){printf("%d = %d \n",i,f[i]);}
}

Problem test

扩展问题1: T组N.求1~N范围上与N互素的数的和。

值得一提的是推导到最后的。按照以往的手段似乎没有继续下去的可能了。(但是如果你仔细观察的话。可以发现 n/k 不需要取底符号。那么就能提取出一个n的因子)

Code:

#include<stdio.h>
#include<string.h>
#define N 100
int num;
int prime[N+5];
int mobius[N+5];
bool mark[N+5];void Mobius()
{int i,j;num = 0;mobius[1] = 1;for(i=2;i<=N;i++){if(!mark[i]){prime[num++] = i;mobius[i] = -1;}for(j=0;j<num;j++){if(i*prime[j]>N){break;}mark[i*prime[j]] = 1;if(i%prime[j]==0){mobius[i*prime[j]] = 0;}else{mobius[i*prime[j]] = -mobius[i];}}}
}
int solve(int n)
{int i,r;r = 0;for(i=1;i<=n;i++){if(n%i==0){r += mobius[i]*i*(n/i)*(n/i+1);}}r /= 2;return r;
}
int main()
{int i,n;Mobius();while(scanf("%d",&n)!=EOF){printf("%d = %d\n",n,solve(n));}
}

Problem two

实际上求互质数和有 n*φ(n)/2 。

用莫比乌斯函数表示

上面公式得证。十分感谢yzq986的留言。告诉了我后续的解法!!!~~~

如果我们直接用n*φ(n)/2。该函数我们是可以直接筛出来的。

对于互质数我们探讨得较多了。个数(欧拉函数)。互素数和。就是以上的。

那么对于约数呢?另外开一篇随笔去探讨这个问题。

论文上的一个优化:

论文上sqrt的优化具体原理论文已经给得很清楚了。

即存在 a/x = a/(x+k)  这个是取整除法

稍微讲述一下代码的构造。

我们预处理出目标函数之后。再预处理出其前缀和用sum数组保存.通过以下代码进行结果的处理。即可。

 for(int i=1,last;i<=n;i=last+1)  {  last = min(n/(n/i),m/(m/i));  ans += (n/i)*(m/i)*(sum[last]-sum[i-1]);  }  

优化

欢迎点个赞~

转载于:https://www.cnblogs.com/Milkor/p/4474835.html

读贾志鹏线性筛有感 (莫比乌斯函数的应用)相关推荐

  1. 线性筛求莫比乌斯函数前缀和

    大家知道有一类问题,让你把1~N之间每一个数的莫比乌斯函数都输出来,或者把它们全部加起来再输出.这种问题应该属于求前缀和一类的啦 用正规的方式表达: 就论求1~N间每个数的莫比乌斯函数吧.我们一般使用 ...

  2. 为什么线性筛欧拉函数i%prime[j]==0的时候phi[i*prime[j]]=phi[i]*prime[j]

    看贾志鹏线性筛的时候想起来的. 我有一个繁琐的证明- -. 证明ϕ(pm)=p×ϕ(m),p为素数,m∈Z\phi(pm)=p\times\phi(m),p为素数,m\in \Bbb Z. 设 m=p ...

  3. 【数学知识】三种方法求 [1,n] 中所有数欧拉函数(线性筛欧拉函数优化至 O(n) )

    整理的算法模板合集: ACM模板 ①直接求小于或等于n,且与n互质的数个数(求[1,n]中所有数的欧拉函数时间复杂度:O(nn)O(n\sqrt{n})O(nn​)) ②求[1,n]之间每个数的质因数 ...

  4. Educational Codeforces Round 37-F.SUM and REPLACE (线段树,线性筛,收敛函数)

    F. SUM and REPLACE time limit per test2 seconds memory limit per test256 megabytes inputstandard inp ...

  5. bzoj2186,P2155-[SDOI2008]沙拉公主的困惑【线性筛,欧拉函数,逆元】

    正题 题目链接: https://www.lydsy.com/JudgeOnline/problem.php?id=2186 https://www.luogu.org/problem/P2155 题 ...

  6. facsum (线性筛 积性函数)

    facsum 10.20 思路: 对于式子进行分析,前面的Phi就是一个很常规的转化,变成i^m.后边一个约数个数,一个Mobius,一个单位函数都是积性函数,所以卷起来也就积性啦(不是积性怎么做嘛. ...

  7. 线性筛求积性函数的模板

    ACM常用模板合集 void sieve() {tot = 1;memset(vis, 0, sizeof(vis));low[1] = 1;G[1] = 函数G(n) n=1时的定义for (int ...

  8. 证明莫比乌斯函数为积性

    线性筛求莫比乌斯函数前缀和_daisysunchaser的博客​​​​​​ 请访问该链接! 今天做一道题,题目详情 - 题目背景能别人帮我写么( - HydroOJ 真tm好烦,劳资一直往大的头上加, ...

  9. 线性筛与欧拉函数、莫比乌斯函数

    网上关于素数筛的资料很多,这里只是给出弱鸟整理的几个线性筛和应用. 最朴素的素数筛--埃拉托斯特尼筛法(Sieve of Eratosthenes) 复杂度 Olognlogn int primes[ ...

最新文章

  1. 求解最大公约数算法(包含四种方法)
  2. spring mvc 异常统一处理方式
  3. socket 与 vfs
  4. Android中selector的使用
  5. Python手册 3.7
  6. 小程序button去除边框
  7. 如何设计出高端大气、有黑科技感的可视化大屏?
  8. leetcode [59] 螺旋矩阵 II / Spiral Matrix II,面试常考题 模拟螺旋过程
  9. 计量经济学 pdf_计量经济学笔记(十六)
  10. 【POJ3335】Rotating Scoreboard(多边形的内核-----半平面交+特殊情况)
  11. iec104规约遥测遥信解析笔记
  12. 删除服务列表中的Tomcat服务?(或删除服务列表中的任意服务)
  13. 国内外五款免费的协作类管理软件推荐
  14. supercharge快充_IQOO、小米快充技术对比,27W快充和44W快充,到底哪家强
  15. java软件工程师工作业绩_JAVA软件工程师简历自我评价
  16. 不可思议:99%的人不了解的真实中国历史
  17. 数字图像处理:4.色彩空间转换
  18. 代码: 0x80131500,Windows应用商店出现错误
  19. 【组队学习】【36期】组队学习内容详情
  20. android HorizontalScrollView讲解

热门文章

  1. php自动关机代码,vb下的恶搞关机程序
  2. linux链接器脚本,vmlinux-lds 连接器脚本
  3. 如何从文件名字符串中获取文件扩展名_Linux操作系统:文件系统的功能和命名...
  4. android ble readcharacteristic,android 4.3蓝牙ble不要调用onCharacteristicRead()
  5. java商品名称_Java统计商品信息
  6. linux操作系统的训练目的,《西南交大-Linux操作系统应用训练》.doc
  7. CO2 convex sets
  8. ResNeX论文概述
  9. mysql 导出所有函数_mysql 导入导出 包括函数或者存储过程
  10. Java面试题-2021