利用辛达拉姆筛进行素数判定

1 常规判定方法

       素数判定问题就是对给定的正整数n判定是否为素数。所谓素数，是指恰好有2个约数的整数。因为n的约数都不超过n，所以只需要检查2~n-1的所有整数是否整除n就能判定是不是素数。不过，我们还能进一步优化。如果d是n的约数，那么n/d也是n的约数。由n=d*n/d可知min(d,n/d)，所以只需要检查2~的所有整数就足够了。此时，素数判定的复杂度为O()。代码实现如下：

int classic(long long a)
{
    long long half=(long long)sqrt(a);

    for(int i=2;i<=half;i++)
    {
        if(a%i==0)  return 0;
    }

    return 1;
}

常规的素数判定思想和实现都非常简单，但是效率较低。现在存在不少更高效的素数判定方法，比如费马测试、算法、数域筛法等，但是这些算法要么不精确要么需要的数学知识过于复杂，在此不做过多介绍。下面我们介绍一种只需要初等数学即可明白的素数判定方法，该方法在任何情况下都会获得比常规方法好一倍的性能。

2 埃氏筛法

       在介绍新方法之前，我们先回顾一下埃氏筛法，这是一个与辗转相除法一样古老的算法。

       大约在公元前3世纪，古希腊数学家埃拉托色尼提出了一种编造素数表的方法（如表一）。这种方法类似于筛东西，把不要的筛掉，把需要的留下。具体做法是：将从2到N的自然数，按顺序排列成2,3,4,5,…,N，然后留下第一个2，划去所有2的倍数；2之后没划去的第一个数是3，留下3，划去所有3的倍数；在3后面没划掉的第一个数是5，留下5，划去所有5的倍数；如此继续，直至上述一列数中再也没有可划的数为止，留下来的便是N以内的一切素数。

表一 埃氏筛法

       如果只对一个整数进行素性判定，常规判定方法或者后面介绍的方法就足够了。但如果要对许多整数进行素性判定，就需要采用埃氏筛法。以下枚举算法摘自《挑战程序设计竞赛》一书：

int prime[MAX_N];
bool is_prime[MAX_N];

//返回n以内的素数个数
int sieve(int n) {
	int p=0;
	for(int i=0;i<=n;i++) is_prime[i]=true;
	is_prime[0]=is_prime[1]=false;
	for(int i=2;i<=n;i++) {
		if(is_prime[i]) {
			prime[p++]=i;
			for(int j=2*i;j<=n;j+=i) is_prime[j]=false;
		}
	}
	return p;
}

3 辛达拉姆筛

       公元1934年，一名年轻的东印度数学生辛达拉姆，提出了一种与埃拉托色尼迥然不同的筛法。它首先列出了一张表，如表二。表的第一行和第一列都是首项为4，公差为3的等差数列。从第二行开始，以后各行也是等差数列，公差分别为5,7,9,11,13……。可以看出，该表其实是一个对称矩阵。

       辛达拉姆指出： 如果N出现在表中，则2N+1是合数；若N不在表中，则2N+1是素数。证明相当精彩！

表二 辛达拉姆筛构造的表

       先证明前半部分。首先，他写出来第m行的第一个数：

4+(m-1)×3=3m+1

注意到该行是公差为2m+1的等差数列，所以此行第n列的数是：

3m+1+(n-1)*(2m+1)=2mn+m+n

即第m行n列的数是2mn+m+n。于是2N+1=4mn+2m+2n+1=(2m+1)(2n+1)是合数。

       再证后半部分。要想正面证明2N+1是素数是相当困难的。如果换成等价的逆否命题，即证“若2N+1不是素数，则N必在表中”似乎要容易得多。事实上，若

       2N+1=x*y(x ,y)为整数

       则因2N+1为奇数，x、y也必为奇数。不妨设：

       x=2p+1;y=2q+1

       从而2N+1=(2p+1)(2q+1)=4pq+2p+2q+1。由此可以得到

       N=2pq+p+q。

也即N是表中第p行第q列的数。

       综上所述，我们证明了辛达拉姆筛法的正确性。例如18不在表中，则2*18+1=37为素数。相反，71在表中，则2*71+1=143是合数。后半部分的证明是反证法的绝佳实例。

4 基于辛达拉姆筛的素数判断   

       有了辛达拉姆筛法，我们该如何判定一个正整数n是否为素数呢？我们只需要判断(n-1)/2是否在表中即可。由于需要判定的数肯定是奇数（偶数肯定都是合数），所以(n-1)/2肯定可以整除。判断是否在表中，也即判断是否存在p和q使得2pq+p+q=(n-1)/2。按照这个思路实现的代码如下：

int xindalamu(long long a)
{
	if(a%2==0) return 0;

    long long aa=(a-1)>>1;
    long long half=(long long)sqrt(aa/2);

    for(int i=1;i<=half;i++)
    {
        if((aa-i)%(2*i+1)==0) return 0;
    }

    return 1;
}

下面详细解释一下上面的代码。首先判断用户输入是否是偶数，如果是直接返回false。然后将用户输入减一除以二赋值给aa，后面就需要判断aa是否在表中。判断的方法和常规判断的方法类似，也是需要验证某一行是否出现整除的情况，不过这里公式稍有不同。从上述代码可以看到，在对aa开方的时候，我们将其除以二。这是为什么呢？判断aa是否出现在表中，也即判断是否存在p和q使得2pq+p+q=aa。这和常规方法中的判断aa=p*q稍有不同。忽略掉两个小项，我们需要判断2pq是否等于aa，也即pq是否等于aa/2，这就是除以2的原因。for循环从第一行开始检查，判断该行是否存在列满足2pq+p+q=aa。当行号固定之后，判断的方法是看aa-p是否整除2p+1（稍作变化即可看出）。

       我们分析一下基于辛达拉姆筛的素数判定复杂度。由for循环可以看出，复杂度也是O( )，不过常系数比常规方法要小。这是因为，在求aa时有个除以二的操作，在开方时又有一个除以二的操作，所以相当于对a/4进行开方，因而最多只需要对 的行进行判断。而常规方法最多则需要验证 个不同的因子，所以可以看出基于辛达拉姆筛的素数判定方法比常规方法快一倍。

       但是，还有一个疑问需要解决，如果用户输入的是一个合数，基于辛达拉姆筛的素数判定方法是否还优于常规判定方法。可以证明，不管输入的是合数还是素数，基于辛达拉姆筛的素数判定方法始终比常规方法快一倍。证明如下：

       如果用户输入的是素数，在分析复杂度的时候已经证明。如果用户输入的是合数，则常规方法在遇到a的最小质因子b时返回，而辛达拉姆判定方法会在遍历到第(b-1)/2行时返回。假设a=b*y，因为a是合数，则肯定存在p行和q列使得2pq+p+q=(a-1)/2。变换得(2p+1)(2q+1)=a。由于辛达拉姆判定方法在遍历时是从第一行开始的，所以我们可以认为b=2p+1，从而辛达拉姆判定方法在第(b-1)/2行时返回，这也就证明了当输入是合数（奇数）时，基于辛达拉姆筛的素数判定方法也比常规方法快一倍。

5 基于辛达拉姆筛的素数筛选

       辛达拉姆筛最初的目的和埃拉托色尼筛一样，也是筛选所有的素数，上面的素数判定算是一个变形应用，筛选所有素数才是筛的最重要目的。不过和埃氏筛法相比，辛达拉姆筛却不那么容易获得所有的素数。从辛达拉姆得到的结论：如果N出现在表中，则2N+1是合数；若N不在表中，则2N+1是素数，我们可以很容易知道不在表中的所有数，其二倍加一包含了所有的奇素数。例如，表中出现的前10个数为4,7,10,12,13,16,17,19,22,24，不在表中的数即为1,2,3,5,6,8,9,11,14,15,18,20,21,23，这些数的二倍加一分别为：3,5,7,11,13,17,19,23,29,31,37,41,43,47，也即前14个奇素数。虽然我们只需要找到所有不在表中的所有数即可获得所有的奇素数，但是如何找到却是个麻烦问题。问题在于如何按序遍历在表中的所有数。本人暂时没有思路，还请高手帮忙解答，在此谢过。

参考

       对于辛达拉姆筛的了解是从高中一本数学习题集中获得的，当时读了练习后面的拓展阅读，印象非常深刻，以至于到现在我都留着那本习题集，感谢它。