自制暴力区间素数筛法的由来

一、探索的缘由

2017年3月31日,教我们C语言的老师谈到如何使用C语言判断一个数是否是质数。

他最开始的思路是判断2~n-1范围的数是否可以整除n。随后老师说由于数学家证明了如果用2到根号n之间的任意整数均无法整除n则n为质数。于是仅需判断2~根号n范围的数是否可以整除n即可。接着老师说:“代码的优化就是这样的,优化代码也是一种乐趣“。这成功激起了我的斗志,以至于接下来的课都没好好听老师讲的内容,也就断断续续听到大致的内容……

下课结束后,我回到宿舍,打开我的笔记本电脑(CPU是Intel Celeron M 420,开发环境Dev-Cpp 5.11),按照老师的思路写了最初的质数判断代码,用于获取0~0xFFFFFF范围内的素数个数。以下是我写的代码:

#include <stdio.h>
#include <math.h>

int IsPrime(int n)
{
	if(n < 2) return 0;
	if(n == 2) return 1;
	int m = sqrt(n) + 1, i;
	for(i = 2; i <= m; ++i)
	{
		if(n % i == 0)
			return 0;
	}
	return 1;
} 

int main()
{
	int i, sum = 0;
	
	for (i = 1; i <= 0xFFFFFF; ++i)
	{
		if(IsPrime(i)) ++sum;
	}

	printf("%d\n", sum);

	return 0;
}

编译运行发现控制台的光标闪烁了很久还没出结果,我还以为自己写的代码也许是因为bug而造成了死循环而又检查了一遍……最后在代码没有发现错误,于是又开始执行……同时我打开手机看会儿知乎以消磨时间……过了许久,结果终于出来了,0~0xFFFFFF范围内的素数个数为1077871。

二、对算法的简单修改

虽然上述的代码可以正确的得出结果,但是效率实在是太低了,在我那台笔记本电脑上花费了1分钟左右才得出了结果……于是势必要探索效率更高的实现。

先开始我通过度娘查了些资料,很多人都推荐使用“欧拉筛法”和“埃拉托斯特尼筛法”……但是并没读懂那些人写的代码(一堆命名为a,b,c的变量;谁看的懂啊……而且我高数学的并不好);于是我接着查资料;发现有些人提出“除2以外的2的倍数都是合数,可以事先筛掉”;由于比较容易理解,所以我立马修改了代码,执行时间缩短到了原来的一半;代码如下

#include <stdio.h>
#include <math.h>

int IsPrime(int n)
{
	if (n < 2) return 0;
	if (n % 2 == 0) return (n == 2);
	int m = sqrt(n) + 1, i;
	for (i = 3; i <= m; i += 2)
	{
		if (n % i == 0)
			return 0;
	}
	return 1;
}
int main()
{
	int i, sum = 0;

	for (i = 1; i <= 0xFFFFFF; ++i)
	{
		if(IsPrime(i)) ++sum;
	}

	printf("%d\n", sum);

	return 0;
}

但即使这样,还是需要35秒才能出结果,还是令我很不爽。于是又开始了新的一波修改之旅。

三、通过概率所思考到的

由于那时已经到中午了,走到食堂点了饭坐下来开始思考:”假设有100个数,是2的倍数概率为50%,是3的倍数的概率为33%;那么出现2的倍数或3的倍数的概率是多少呢?那么如果把这些数剔除掉的话,判断的概率应该会大幅度减小“

把饭吃完后,坐到笔记本前写了个小程序判断了下,大致是65%的概率……觉得这个概率还是太低于是增加了更多的质数;发现在为2,3,5,7的倍数的数大致占了80%;接着我把代码稍微改了下:

#include <stdio.h>
#include <math.h>

int IsPrime(int n)
{
	if(n < 2) return 0;
	if(n % 2 == 0) return (n == 2);
	if(n % 3 == 0) return (n == 3);
	if(n % 5 == 0) return (n == 5);
	if(n % 7 == 0) return (n == 7);

	int m = sqrt(n) + 1, i;
	for(i = 11; i <= m; i += 2)
	{
		if(n % i == 0)
			return 0;
	}
	return 1;
} 

int main()
{
	int i, sum = 0;

	for (i = 1; i <= 0xFFFFFF; ++i)
	{
		if(IsPrime(i)) ++sum;
	}

	printf("%d\n", sum);

	return 0;
}

运行后,可喜的发现,只需要8秒就可以出结果了;于是我在思考,如果排除掉更多的质数的倍数的话会如何;先开始在原来基础上修改;发现由于逻辑的增加效率反而变低了。遗憾的是接下来要上高数课了,于是只好中断思路……晚上回家打开我的台式机上的Visual Studio 2017(CPU是Intel Xeon E3 1230 v2)继续探索。首先测试了下之前写的代码所需要的时间。

  • 老师的思路(判断2~根号n范围的数是否可以整除n即可)需要10秒
  • 基于老师的思路,剔掉2的倍数的合数需要5秒

接着我开始使用查表法。经过我的反复实验,得出了获取0~0xFFFFFF范围内的素数个数效率最高的质数表;以下是通过素数表获取素数个数的代码:

#include <stdio.h>
#include <math.h>

const short PrimeList[] =
{
	2,3,5,7,11,13,17,19,
	23,29,31,37,41,43,47,53,
	59,61,67,71,73,79,83,89,
	97,101,103,107,109,113,127,131,
	137,139,149,151,157,163,167,173,
	179,181,191,193,197,199,211,223,
	227,229,233,239,241,251,257,263,
	269,271,277,281,283,293,307,311,
	313,317,331,337,347,349,353,359,
	367,373,379,383,389,397,401,409,
	419,421,431,433,439,443,449,457,
	461,463,467,479,487,491,499,503,
	509,521,523,541,547,557,563,569,
	571,577,587,593,599,601,607,613,
	617,619,631,641,643,647,653,659,
	661,673,677,683,691,701,709,719,
	727,733,739,743,751,757,761,769,
	773,787,797,809,811,821,823,827,
	829,839,853,857,859,863,877,881,
	883,887,907,911,919,929,937,941,
	947,953,967,971,977,983,991,997,
	1009,1013,1019,1021,1031,1033,1039,1049,
	1051,1061,1063,1069,1087,1091,1093,1097,
	1103,1109,1117,1123,1129,1151,1153,1163,
	1171,1181,1187,1193,1201,1213,1217,1223,
	1229,1231,1237,1249,1259,1277,1279,1283,
	1289,1291,1297,1301,1303,1307,1319,1321,
	1327,1361,1367,1373,1381,1399,1409,1423,
	1427,1429,1433,1439,1447,1451,1453,1459,
	1471,1481,1483,1487,1489,1493,1499,1511,
	1523,1531,1543,1549,1553,1559,1567,1571,
	1579,1583,1597,1601,1607,1609,1613,1619,
	1621,1627,1637,1657,1663,1667,1669,1693,
	1697,1699,1709,1721,1723,1733,1741,1747,
	1753,1759,1777,1783,1787,1789,1801,1811,
	1823,1831,1847,1861,1867,1871,1873,1877,
	1879,1889,1901,1907,1913,1931,1933,1949,
	1951,1973,1979,1987,1993,1997,1999,2003,
	2011,2017,2027,2029,2039,2053,2063,2069,
	2081,2083,2087,2089,2099,2111,2113,2129,
	2131,2137,2141,2143,2153,2161,2179,2203,
	2207,2213,2221,2237,2239,2243,2251,2267,
	2269,2273,2281,2287,2293,2297,2309,2311,
	2333,2339,2341,2347,2351,2357,2371,2377,
	2381,2383,2389,2393,2399,2411,2417,2423,
	2437,2441,2447,2459,2467,2473,2477,2503,
	2521,2531,2539,2543,2549,2551,2557,2579,
	2591,2593,2609,2617,2621,2633,2647,2657,
	2659,2663,2671,2677,2683,2687,2689,2693,
	2699,2707,2711,2713,2719,2729,2731,2741,
	2749,2753,2767,2777,2789,2791,2797,2801,
	2803,2819,2833,2837,2843,2851,2857,2861,
	2879,2887,2897,2903,2909,2917,2927,2939,
	2953,2957,2963,2969,2971,2999,3001,3011,
	3019,3023,3037,3041,3049,3061,3067,3079,
	3083,3089,3109,3119,3121,3137,3163,3167,
	3169,3181,3187,3191,3203,3209,3217,3221,
	3229,3251,3253,3257,3259,3271,3299,3301,
	3307,3313,3319,3323,3329,3331,3343,3347,
	3359,3361,3371,3373,3389,3391,3407,3413,
	3433,3449,3457,3461,3463,3467,3469,3491
}; // 3499

int IsPrime(int n)
{
	if(n < 2) return 0;

	int i;

	for(i = 0; i < sizeof(PrimeList) / sizeof(PrimeList[0]); ++i)
		if(n % PrimeList[i] == 0)
			return (n == PrimeList[i]);

	int m = sqrt(n) + 1;
	for(i = 3499; i <= m; i += 2)
	{
		if(n % i == 0)
			return 0;
	}
	return 1;
}

int main()
{
	int i, sum = 0;

	for (i = 1; i <= 0xFFFFFF; ++i)
	{
		if(IsPrime(i)) ++sum;
	}

	printf("%d\n", sum);

	return 0;
}

通过2.22秒就可以获取结果,我于是打算扩增了获取规模,变成获取0~0x3FFFFFFF之间的素数个数,规模整整扩大了64倍……编译运行后好长时间(2分钟)还没出结果;于是我等不及关掉了。把获取范围又变成原来的0~0xFFFFFF。

重新打开百度搜索查看资料以搜寻效率更高的质数判断算法,获取了以下经验:

  1. 发现我上午的思路原来和“埃拉托斯特尼筛法”是一致的
  2. 欧拉筛法(线性筛法)效率最高

但是我查了很多有关欧拉筛法的资料,直到第二天的凌晨依旧没看懂欧拉筛法的思路……于是恋恋不舍的上了床。依旧还在思考如何改进……

四、暴力区间素数筛法

早晨起床,打开电脑,突然有了灵感;对原来的代码进行了重写……编译执行惊异的发现获取0~0xFFFFFF范围内的素数个数只需要0.3秒就可以出结果。于是把规模扩大64倍后也只需要10秒(获取0~0x3FFFFFFF之间的素数个数)

于是我兴奋的把成果告诉TooYoungTooSimp(4071E95D-A09B-4AA3-8008-);他对我的算法表示很感兴趣;但是我虽然知道但没法说清楚思路……于是向他提供代码且希望他帮我解释下原理。他看了下代码,做了以下解释

  • 初看下来你是分了两个步骤:第一个是2~sqrt(n)是一般的暴力,然后第二部分是埃拉托斯特尼筛法
  • 一个数不会被他的平方根以内的数整除是这个数为质数的充要条件,这是第一部分的正确性;第二部分就是埃拉托斯特尼筛法了

他不仅做了解释,于是提出”如果单纯用埃拉托斯特尼筛法“会更快的提议;我以完全使用埃拉托斯特尼筛法获取0~0xFFFFFF范围内的质数个数需要10秒就可以出结果的实例反驳了他;同时也终于和他说出了我的思路:

一个数不会被它的平方根以内的数整除是这个数为质数的充要条件,那么可以根据它的平方根以内的数筛选掉该范围内额所有不是质数的数;这样规模就变成了根号n;然后2到根号n的所有数都能由其范围内的所有素数组合得到。于是可得出:通过2到根号n的所有素数可以得出2~n范围内的所有不是素数的数,那么久筛选掉他们

他对我的思路做了以下评价:

通过时间复杂度更高O(n\sqrt{n})(是这么多吧)
但是常数更小的朴素方式求出了小规模内的质数
然后再用时间复杂度较低的埃氏筛法求出之后的部分
看来暴力也是可以出奇迹的(不过都不如欧拉筛)

我也顺便向他请教了欧拉筛法的原理,他和我讲了很久;但我仍似懂非懂……于是他发过了他写的代码……读了一会儿;看懂了;欧拉筛法之所以效率高,是因为质数单独一个存放到一个数组中,这样可以使每个合数都能被其最小质数筛掉。

他希望能让我帮他测试下欧拉筛法在我的机器上速度如何;发现欧拉筛法只要7秒;比我的筛法快2秒但偶然一瞥VS调试器显示的内存占用;把我着实吓着了(获取0~0x3FFFFFFF之间的素数个数,我的方法需要1GB内存,而他写的高效率欧拉筛法需要1.3GB内存)。于是对他写的欧拉筛法进行了魔改(在使用位操作的情况下,又去掉了欧拉筛法存放质数结果的数组(毕竟几千万个int太浪费内存了和我思考出的筛法通过位操作缩减内存占用。

当我把魔改后欧拉筛法代码给他看后,他说你这样毁掉了欧拉筛法(当然事实的确如此);我那时继续查资料研究欧拉筛法(当时我还不知道为什么不能去掉存放质数结果的数组);结果偶然发现我的素数筛法应称为区间素数筛法,其定义是:

因为b以内合数的最小质因数一定不超过sqrt(b),如果有sqrt(b)以内的素数表的话,就可
以把筛选法用在[a,b)上了,先分别做好[2,sqrt(b))的表和[a,b)的表,然后从[2,sqrt(b))
的表中筛得素数的同时,也将其倍数从[a,b)的表中划去,最后剩下的就是区间[a,b)内的
素数了。

但由于和真正的区间素数筛法相比起来有些暴力;于是我把我的筛法命名为暴力区间素数筛法;经过和TooYoungTooSimp(4071E95D-A09B-4AA3-8008-)的大量讨论,最终代码是这样的:

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define INT_BITS (sizeof(unsigned long) * 8)

#define GET_BIT_VALUE(Array, Index) \
	(bool)((Array[Index / INT_BITS] >> (Index % INT_BITS)) & 1)

#define SET_BIT_TRUE(Array, Index) \
	(Array[Index / INT_BITS] |= (1 << (Index % INT_BITS)))

#define SET_BIT_FALSE(Array, Index) \
	(Array[Index / INT_BITS] &= ~(1 << (Index % INT_BITS)))

/*
该函数获取小于等于MaxNum之间自然数的所有素数
返回值是一个用malloc分配的bool数组(里面包含了该范围内任意数是否为素数的判断)

假设我要获取0到n范围自然数内的所有素数

首先搞定特殊情况,0和1不是质数

因为可以证明:
1.对正整数n,如果用2到根号n之间的任意整数均无法整除n,则n为质数
2.任何一个合数都可以分解为几个质数的积

所以可得到
1.2到n之间的任意所有合数都能由2到根号n内的任意整数组合得到
2.2到根号n之间的所有合数都能由2到根号n内的任意素数组合得到

综上所述:
2到n之间的所有合数都能由2到根号n内的任意素数组合得到

本算法步骤:
1.先以2到根号n内的任意素数组成一个筛子
2.通过刚刚组成的筛子筛掉2到n之间的所有合数,剩下来的就是素数

总之,该算法相当于埃拉托斯特尼筛法的改进:

由于我不知道我改进的地方,
问了下4071E95D-A09B-4AA3-8008-,他的回答是

通过时间复杂度更高O(n\sqrt{n})(是这么多吧)
但是常数更小的朴素方式求出了小规模内的质数
然后再用时间复杂度较低的埃氏筛法求出之后的部分

在百度查了下,这个算法应称为区间素数筛法,其定义是:
因为b以内合数的最小质因数一定不超过sqrt(b),如果有sqrt(b)以内的素数表的话,就可
以把筛选法用在[a,b)上了,先分别做好[2,sqrt(b))的表和[a,b)的表,然后从[2,sqrt(b))
的表中筛得素数的同时,也将其倍数从[a,b)的表中划去,最后剩下的就是区间[a,b)内的
素数了。
*/


/*
使用改进了内存占用的暴力区间素数筛法获取0到MaxNum范围内的素数情况
返回的是用malloc分配的记录素数情况的位图
*/
unsigned long* GetPrimeList0(unsigned long MaxNum)
{
    // 计算内存块大小
    size_t size = ((MaxNum + 1) >> 3) + 1;

    // 分配一块内存
    unsigned long *BitArray = (unsigned long*)malloc(size);
    if (BitArray)
    {
        // 内存块初始化,假定该范围都是素数
        memset(BitArray, -1, size);

        // 0和1不是质数
        SET_BIT_FALSE(BitArray, 0);
        SET_BIT_FALSE(BitArray, 1);

        // 判断[2,sqrt(MaxNum))范围内的素数情况
        unsigned long MaxRange = sqrt(MaxNum);
        for (unsigned long CurNum = 2; CurNum <= MaxRange; ++CurNum)
        {
            // 判断CurNum是否为质数
            unsigned long CurMaxRange = sqrt(CurNum);
            for (unsigned long i = 2; i < CurMaxRange; ++i)
            {
                if (CurNum % i == 0)
                {
                    SET_BIT_FALSE(BitArray, CurNum);
                    break;
                }
            }

            // 如果CurNum不为质数,则跳过
            if (!GET_BIT_VALUE(BitArray, CurNum)) continue;

            /*
            如果CurNum为质数,则用其筛掉[0,MaxNum}范围内可以整除该质数的合数
	
            初始化为CurNum^2的原因是:
            n能把n^2以上的合数筛掉,[n,n^2]之间的会被比n小的质数筛掉
            */
            for (unsigned long i = CurNum * CurNum; i <= MaxNum; i += CurNum)
            {
                if (GET_BIT_VALUE(BitArray, i)) SET_BIT_FALSE(BitArray, i);
            }
        }
    }
	return BitArray;
}

int main()
{	
	unsigned long sum = 0;
	//0x3FFFFFFF
	unsigned long* BitMap = GetPrimeList0(0x3FFFFFFF);

	if (BitMap)
	{
		for (unsigned long i = 0; i < 0x3FFFFFFF + 1; ++i)
		{
			if (GET_BIT_VALUE(BitMap, i)) ++sum;
		}

		printf("%ld\n", sum);

		free(BitMap);
	}

	return 0;
}

五、效率评测和评价

测试环境

  • Intel Xeon E3 1230 V2
  • 8GB DDR3 1600 x2
  • Windows 10 Creator Update Build 15063
  • Visual Studio 2017 Community

在以上环境获取小于等于0x3FFFFFFF的自然数的的素数个数

如果不开启SSE和SSE2的指令集编译优化且都使用位操作节省内存,结果如下

  • 我的魔改版欧拉筛法 内存占用142mb 耗时12031ms
  • 暴力区间素数筛法(也就是本说说分享的方法) 内存占用142mb 耗时9878ms
  • TooYoungTooSimp写的原版欧拉筛法 内存占用399mb 耗时7984ms

相对来说,虽然还是欧拉筛法最快,但是比较耗内存……

如果对内存占用比较在乎而不介意稍微比欧拉筛法慢一些的话,我觉得可以试试暴力区间素数筛法(就当我王婆卖瓜,自卖自夸)

六、读者引用代码须知

这段代码,其实是我和TooYoungTooSimp一起探讨完成的

本代码你们随便用,没有任何协议限制

但如果要提及作者的话,请务必把两个人都提到

七、后续

我在我的那台CPU为AMD闪龙3000+且只有512MB内存的安装Windows XP x64 Edition台式机,测试了暴力区间素数筛法,260秒后顺利出结果(但由于机器比较老,不久就过热重启了)

清明节过后的C语言课课后我和老师讨论自己对以上内容的探究,老师对我的探究表示了认可

几天后参加蓝桥杯省赛,刚好碰到和素数有关的题目;我使用暴力区间素数筛法得出了正确结果(由于比赛用机是和E3 1230 v2同架构的i5-3470;于是也是秒出结果)

PS1:由于我并不是计算机学院的学生,于是本文写的并不是很专业;请见谅。

PS2:如果本文有不妥之处,希望帮忙指出


毛利,写于2017年4月22日

编辑于 2017-04-23

文章被以下专栏收录