[原创]浅析求素数算法

chai2010 · 发表于 2006-10-27 15:54:46

浅析求素数算法

时间: 2006-10-27

注意: 如果没有特殊说明, 以下讨论的都是针对n为素数时的时间复杂度

1. 根据概念判断:

如果一个正整数只有两个因子, 1和p，则称p为素数.

bool isPrime(int n)
{
if(n < 2) return false;
for(int i = 2; i < n; ++i)
if(n%i == 0) return false;
return true;
}

复制代码

时间复杂度O(n).

2. 改进, 去掉偶数的判断

bool isPrime(int n)
{
if(n < 2) return false;
if(n == 2) return true;
for(int i = 3; i < n; i += 2)
if(n%i == 0) return false;
return true;
}

复制代码

时间复杂度O(n/2), 速度提高一倍.

3. 进一步减少判断的范围

定理: 如果n不是素数, 则n有满足1<d<=sqrt(n)的一个因子d.
证明: 如果n不是素数, 则由定义n有一个因子d满足1<d<n.
如果d大于sqrt(n), 则n/d是满足1<n/d<=sqrt(n)的一个因子.

bool isPrime(int n)
{
if(n < 2) return false;
if(n == 2) return true;
for(int i = 3; i*i <= n; i += 2)
if(n%i == 0) return false;
return true;
}

复制代码

时间复杂度O(sqrt(n)/2), 速度提高O((n-sqrt(n))/2).

4. 剔除因子中的重复判断.
例如: 11%3 != 0 可以确定 11%(3*i) != 0.

定理: 如果n不是素数, 则n有满足1<d<=sqrt(n)的一个"素数"因子d.
证明: I1. 如果n不是素数, 则n有满足1<d<=sqrt(n)的一个因子d.
I2. 如果d是素数, 则定理得证, 算法终止.
I3. 令n=d, 并转到步骤I1.

由于不可能无限分解n的因子, 因此上述证明的算法最终会停止.

// primes[i]是递增的素数序列: 2, 3, 5, 7, ...
// 更准确地说primes[i]序列包含1->sqrt(n)范围内的所有素数
bool isPrime(int primes[], int n)
{
if(n < 2) return false;
for(int i = 0; primes[i]*primes[i] <= n; ++i)
if(n%primes[i] == 0) return false;
return true;
}

复制代码

假设n范围内的素数个数为PI(n), 则时间复杂度O(PI(sqrt(n))).

函数PI(x)满足素数定理: ln(x)-3/2 < x/PI(x) < ln(x)-1/2, 当x >= 67时.

因此O(PI(sqrt(n)))可以表示为O(sqrt(x)/(ln(sqrt(x))-3/2)),

O(sqrt(x)/(ln(sqrt(x))-3/2))也是这个算法的空间复杂度.

5. 构造素数序列primes: 2, 3, 5, 7, ...

由4的算法我们知道, 在素数序列已经被构造的情况下, 判断n是否为素数效率很高;

但是, 在构造素数序列本身的时候, 是否也可是达到最好的效率呢?

事实上这是可以的! -- 我们在构造的时候完全可以利用已经被构造的素数序列!

假设我们已经我素数序列: p1, p2, .. pn

现在要判断pn+1是否是素数, 则需要(1, sqrt(pn+1)]范围内的所有素数序列,

而这个素数序列显然已经作为p1, p2, .. pn的一个子集被包含了!

// 构造素数序列primes[]

void makePrimes(int primes[], int num)

{

int i, j, cnt;

primes[0] = 2;

primes[1] = 3;

for(i = 5, cnt = 2; cnt < num; i += 2)

{

      int flag = true;

      for(j = 1; primes[j]*primes[j] <= i; ++j)

      {

         if(i%primes[j] == 0)

         {

            flag = false; break;

         }

      }

      if(flag) primes[cnt++] = i;

}

}

复制代码

makePrimes的时间复杂度比较复杂, 而且它只有在初始化的时候才被调用一次.

在一定的应用范围内, 我们可以把近似认为makePrimes需要常数时间.

在后面的讨论中, 我们将探讨一种对计算机而言更好的makePrimes方法.

6. 更好地利用计算机资源...

当前的主流PC中, 一个整数的大小为2^32. 如果需要判断2^32大小的数是否为素数,

则可能需要测试[2, 2^16]范围内的所有素数(2^16 == sqrt(2^32)).

由4中提到的素数定理我们可以大概确定[2, 2^16]范围内的素数个数.

由于2^16/(ln(2^16)-1/2) = 6138, 2^16/(ln(2^16)-3/2) = 6834,

我们可以大概估计出[2, 2^16]范围内的素数个数6138 < PI(2^16) < 6834.

在对[2, 2^16]范围内的素数进行统计, 发现只有6542个素数:

p_6542: 65521, 65521^2 = 4293001441 < 2^32, (2^32 = 4294967296)
p_6543: 65537, 65537^2 = 4295098369 > 2^32, (2^32 = 4294967296)

在实际运算时unsigned long x = 4295098369;将发生溢出, 为131073.

在程序中, 我是采用double类型计算得到的结果.

分析到这里我们可以看到, 我们只需要缓冲6543个素数, 我们就可以采用4中的算法

高效率地判断[2, 2^32]如此庞大范围内的素数!

(原本的2^32大小的问题规模现在已经被减小到6543规模了!)

虽然用现在的计算机处理[2, 2^16]范围内的6542个素数已经没有一点问题,

虽然makePrimes只要被运行一次就可以, 但是我们还是考虑一下是否被改进的可能?!

我想学过java的人肯定想把makePrimes作为一个静态的初始化实现, 在C++中也可以

模拟java中静态的初始化的类似实现:

#define NELEMS(x) ((sizeof(x)) / (sizeof((x)[0])))

static int primes[6542+1];
static struct _Init { _Init(){makePrimes(primes, NELEMS(primes);} } _init;

如此, 就可以在程序启动的时候自动掉用makePrimes初始化素数序列.

但, 我现在的想法是: 为什么我们不能在编译的时候调用makePrimes函数呢?

完全可以!!! 代码如下:

// 这段代码可以由程序直接生成

const static int primes[] =

{

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,

...

65521, 65537

};

复制代码

有点不可思议吧, 原本makePrimes需要花费的时间复杂度现在真的变成O(1)了!

(我觉得叫O(0)可能更合适!)

7. 二分法查找

现在我们缓存了前大约sqrt(2^32)/(ln(sqrt(2^32)-3/2))个素数列表, 在判断2^32级别的

素数时最多也只需要PI(sqrt(2^32))次判断(准确值是6543次), 但是否还有其他的方式判断呢?

当素数比较小的时候(不大于2^16), 是否可以直接从缓存的素数列表中直接查询得到呢?

答案是肯定的! 由于primes是一个有序的数列, 因此我们当素数小于2^16时, 我们可以直接

采用二分法从primes中查询得到(如果查询失败则不是素数).

// 缺少的代码请参考前边

#include <stdlib.h>

static bool cmp(const int *p, const int *q)

{

return (*p) - (*q);

}

bool isPrime(int n)

{

if(n < 2) return false;

if(n == 2) return true;

if(n%2 == 0) return false;

if(n >= 67 && n <= primes[NELEMS(primes)-1])

{

      return NULL !=

         bsearch(&n, primes, NELEMS(primes), sizeof(n), cmp);

}

else

{

      for(int i = 1; primes[i]*primes[i] <= n; ++i)

         if(n%primes[i] == 0) return false;

      return true;

}

}

复制代码

时间复杂度:

  if(n <= primes[NELEMS(primes)-1] && n >= 67): O(log2(NELEMS(primes))) < 13;
  if(n >  primes[NELEMS(primes)-1]): O(PI(sqrt(n))) <= NELEMS(primes).

8. 素数定理+2分法查找

在7中, 我们对小等于primes[NELEMS(primes)-1]的数采用2分法查找进行判断.

我们之前针对2^32缓冲的6453个素数需要判断的次数为13次(log2(1024*8) == 13).

对于小的素数而言(其实就是2^16范围只内的数), 13次的比较已经完全可以接受了.

不过根据素数定理: ln(x)-3/2 < x/PI(x) < ln(x)-1/2, 当x >= 67时, 我们依然

可以进一步缩小小于2^32情况的查找范围(现在是0到NELEMS(primes)-1范围查找).

我们需要解决问题是(n <= primes[NELEMS(primes)-1):

如果n为素数, 那么它在素数序列可能出现的范围在哪?

---- (n/(ln(n)-1/2), n/(ln(n)-3/2)), 即素数定理!

上面的代码修改如下:

bool isPrime(int n)

{

if(n < 2) return false;

if(n == 2) return true;

if(n%2 == 0) return false;

int hi = (int)ceil(n/(ln(n)-3/2));

if(n >= 67 && hi < NELEMS(primes))

{

      int lo = (int)floor(n/(ln(n)-1/2));

      return NULL !=

         bsearch(&n, primes+lo, hi-lo, sizeof(n), cmp);

}

else

{

      for(int i = 1; primes[i]*primes[i] <= n; ++i)

         if(n%primes[i] == 0) return false;

      return true;

}

}

复制代码

时间复杂度:

  if(n <= primes[NELEMS(primes)-1] && n >= 67): O(log2(hi-lo))) < ???;
  if(n >  primes[NELEMS(primes)-1]): O(PI(sqrt(n))) <= NELEMS(primes).

9. 打包成素数库(给出全部的代码)

到目前为止, 我已经给出了我所知道所有改进的方法(如果有人有更好的算法感谢告诉我).

这里需要强调的一点是, 这里讨论的素数求法是针对0-2^32范围的数而言, 至于像寻找

成百上千位大小的数不在此讨论范围, 那应该算是纯数学的内容了.

代码保存在2个文件: prime.h, prime.cpp.

// file: prime.h

#ifndef PRIME_H_2006_10_27_

#define PRIME_H_2006_10_27_

extern int  Prime_max(void);       // 素数序列的大小

extern int  Prime_get (int i);       // 返回第i个素数, 0 <= i < Prime_max

extern bool Prime_test(int n);       // 测试是否是素数, 1 <= n < INT_MAX

#endif

///////////////////////////////////////////////////////

// file: prime.cpp

#include <assert.h>

#include <limits.h>

#include <math.h>

#include <stdlib.h>

#include "prime.h"

// 计算数组的元素个数

#define NELEMS(x) ((sizeof(x)) / (sizeof((x)[0])))

// 素数序列, 至少保存前6543个素数!

static const int primes[] =

{

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,

...

65521, 65537

};

// bsearch的比较函数

static int cmp(const void *p, const void *q)

{

return (*(int*)p) - (*(int*)q);

}

// 缓冲的素数个数

int Prime_max()

{

return NELEMS(primes);

}

// 返回第i个素数

int Prime_get(int i)

{

assert(i >= 0 && i < NELEMS(primes));

return primes[i];

}

// 测试n是否是素数

bool Prime_test(int n)

{

assert(n > 0);

if(n < 2) return false;

if(n == 2) return true;

if(!(n&1)) return false;

// 如果n为素数, 则在序列hi位置之前

int lo, hi = (int)ceil(n/(log(n)-3/2.0));

if(hi < NELEMS(primes))

{

      // 确定2分法查找的范围

      // 只有n >= 67是才满足素数定理

      if(n >= 67) lo = (int)floor(n/(log(n)-1/2.0));

      else { lo = 0; hi = 19; }

      // 查找成功则为素数

      return NULL !=

         bsearch(&n, primes+lo, hi-lo, sizeof(n), cmp);

}

else

{

      // 不在保存的素数序列范围之内的情况

      for(int i = 1; primes[i]*primes[i] <= n; ++i)

         if(n%primes[i] == 0) return false;

      return true;

}

}

复制代码

10. 回顾, 以及推广

到这里, 关于素数的讨论基本告一段落. 回顾我们之前的求解过程, 我们会发现

如果缺少数学的基本知识会很难设计好的算法; 但是如果一味地只考虑数学原理,

而忽律了计算机的本质特征, 也会有同样的问题.

一个很常见的例子就是求Fibonacci数列. 当然方法很多, 但是在目前的计算机中

都没有实现的必要!

因为Fibonacci数列本身是指数增长的, 32位的有符号整数所能表示的位置只有前46个:

static const int Fibonacci[] =

{

0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,

2584,4181,6765,10946,17711,28657,46368,75025,121393,196418,

317811,514229,832040,1346269,2178309,3524578,5702887,9227465,

14930352,24157817,39088169,63245986,102334155,165580141,267914296,

433494437,701408733,1134903170,1836311903,-1323752223

};

复制代码

因此, 我只需要把前46个Fibonacci数保存到数组中就可以搞定了!

比如: F(int i) = {return Fibonacci;} 非常简单, 效率也非常好.

同样的例子如求阶乘n!, 虽然也有很多数学上的描述, 但是计算机一般都表示不了,

我们直接把能用的计算好放到内存中就可以了.

总之, 许多东西本身是好的, 但是不要被它束缚了!

(完)

chai2010 · 发表于 2006-10-27 16:49:21

大家是否有其他好的方法贴出来探讨一下啊

realtang · 发表于 2006-10-27 22:49:58

真的挺不错的。

soloforce · 发表于 2006-10-27 23:20:06

不错。素数判定是个很有趣的话题~

可以参考一下 http://bbs.linuxsir.cn/showthread.php?t=267468

chai2010 · 发表于 2006-10-30 09:04:04

Post by Lolita
不错。素数判定是个很有趣的话题~

可以参考一下 http://bbs.linuxsir.cn/showthread.php?t=267468

今天又学了不少东西:-)
感谢Lolita斑竹!!

chai2010 · 发表于 2006-11-2 17:08:39

最近看了一些资料
又想了一个笨的方法，
还是测试2^32以内的素数：

先根据费马小定理测试，
为了避免伪素数的情况，
我们保存1-2^32以内的所有Carmichael数，
具体有1118个（c1118: 4277982241）。

这样可能会有更好的效率:-)

chai2010 · 发表于 2006-11-3 17:36:21

Post by chai2010
最近看了一些资料
又想了一个笨的方法，
还是测试2^32以内的素数：

先根据费马小定理测试，
为了避免伪素数的情况，
我们保存1-2^32以内的所有Carmichael数，
具体有1118个（c1118: 4277982241）。

这样可能会有更好的效率:-)

现在给出具体的代码：
(转载请署名作者)

// 函数: bool isPrime(unsigned n);
// 描述: 判断n是否为素数
// 算法: 费马小定理-卡尔麦克数
// 作者: 柴树杉
// 主页: [url]chaishushan.googlepages.com[/url]
// 时间: 2006-11-04
#include <assert.h>
#include <stdlib.h>
// 计算数组的元素个数
#define NELEMS(x) ((sizeof(x)) / (sizeof((x)[0])))
// 2^32范围内的所有卡尔麦克数
static const unsigned CarmichaeNumbers[] =
{
561,1105,1729,2465,2821,6601,8911,10585,15841,29341,41041,46657,52633,62745,
63973,75361,101101,115921,126217,162401,172081,188461,252601,278545,294409,
314821,334153,340561,399001,410041,449065,488881,512461,530881,552721,656601,
658801,670033,748657,825265,838201,852841,997633,1024651,1033669,1050985,
1082809,1152271,1193221,1461241,1569457,1615681,1773289,1857241,1909001,
2100901,2113921,2433601,2455921,2508013,2531845,2628073,2704801,3057601,
3146221,3224065,3581761,3664585,3828001,4335241,4463641,4767841,4903921,
4909177,5031181,5049001,5148001,5310721,5444489,5481451,5632705,5968873,
6049681,6054985,6189121,6313681,6733693,6840001,6868261,7207201,7519441,
7995169,8134561,8341201,8355841,8719309,8719921,8830801,8927101,9439201,
9494101,9582145,9585541,9613297,9890881,10024561,10267951,10402561,10606681,
10837321,10877581,11119105,11205601,11921001,11972017,12261061,12262321,
12490201,12945745,13187665,13696033,13992265,14469841,14676481,14913991,
15247621,15403285,15829633,15888313,16046641,16778881,17098369,17236801,
17316001,17586361,17812081,18162001,18307381,18900973,19384289,19683001,
20964961,21584305,22665505,23382529,25603201,26280073,26474581,26719701,
26921089,26932081,27062101,27336673,27402481,28787185,29020321,29111881,
31146661,31405501,31692805,32914441,33302401,33596641,34196401,34657141,
34901461,35571601,35703361,36121345,36765901,37167361,37280881,37354465,
37964809,38151361,38624041,38637361,39353665,40160737,40280065,40430401,
40622401,40917241,41298985,41341321,41471521,42490801,43286881,43331401,
43584481,43620409,44238481,45318561,45877861,45890209,46483633,47006785,
48321001,48628801,49333201,50201089,53245921,53711113,54767881,55462177,
56052361,58489201,60112885,60957361,62756641,64377991,64774081,65037817,
65241793,67371265,67653433,67902031,67994641,68154001,69331969,70561921,
72108421,72286501,74165065,75151441,75681541,75765313,76595761,77826001,
78091201,78120001,79411201,79624621,80282161,80927821,81638401,81926461,
82929001,83099521,83966401,84311569,84350561,84417985,87318001,88689601,
90698401,92625121,93030145,93614521,93869665,94536001,96895441,99036001,
99830641,99861985,100427041,101649241,101957401,102090781,104404861,104569501,
104852881,105117481,105309289,105869401,106041937,107714881,109393201,109577161,
111291181,114910489,115039081,115542505,116682721,118901521,119327041,120981601,
121247281,122785741,124630273,127664461,128697361,129255841,129762001,130032865,
130497361,132511681,133205761,133344793,133800661,134809921,134857801,135556345,
136625941,139592101,139952671,140241361,144218341,145124785,146843929,150846961,
151530401,151813201,153589801,153927961,157731841,158404141,158864833,159492061,
161035057,161242705,161913961,163954561,167979421,168659569,169057801,169570801,
170947105,171454321,171679561,172290241,172430401,172947529,173085121,174352641,
175997185,176659201,178451857,178482151,178837201,180115489,181154701,182356993,
184353001,186393481,186782401,187188001,188516329,188689501,189941761,193708801,
193910977,194120389,194675041,196358977,200753281,206955841,208969201,212027401,
213835861,214850881,214852609,216821881,221884001,225745345,226509361,227752993,
228842209,230630401,230996949,231194965,237597361,238244041,238527745,241242001,
242641153,246446929,247095361,250200721,252141121,255160621,256828321,257495641,
258634741,266003101,270857521,271481329,271794601,273769921,274569601,275283401,
277241401,278152381,279377281,280067761,280761481,288120421,289766701,289860481,
291848401,292244833,292776121,295643089,295826581,296559361,299736181,300614161,
301704985,302751505,306871201,311388337,318266641,321197185,321602401,325546585,
328573477,329769721,333065305,333229141,334783585,338740417,346808881,348612265,
354938221,357277921,357380101,358940737,360067201,362569201,364590721,366532321,
366652201,367804801,367939585,368113411,382304161,382536001,390489121,392099401,
393122521,393513121,393716701,395044651,395136505,396262945,399906001,403043257,
403317421,405739681,413058601,413138881,413631505,416964241,417241045,419520241,
426821473,429553345,434330401,434932961,438359041,440306461,440707345,455106601,
458368201,461502097,461854261,462199681,471441001,471905281,473847121,477726145,
481239361,483006889,484662529,490099681,490503601,492559141,496050841,499310197,
503758801,507726901,509033161,510825601,511338241,516684961,517937581,518117041,
518706721,527761081,529782121,530443201,532758241,533860309,540066241,542497201,
544101481,545363281,545570641,547652161,548871961,549117205,549333121,549538081,
551672221,552894301,555465601,556199281,556450777,557160241,557795161,558570961,
558977761,561481921,561777121,564651361,568227241,569332177,573896881,577240273,
579606301,580565233,590754385,593234929,595405201,597717121,600892993,602074585,
602426161,602593441,606057985,609865201,611397865,612347905,612816751,616463809,
618068881,620169409,621101185,625060801,625482001,629692801,631071001,633639097,
638959321,642708001,652969351,656187001,662086041,672389641,683032801,683379841,
684106401,686059921,689880801,697906561,698548201,702683101,703995733,704934361,
705101761,707926801,710382401,710541481,711374401,713588401,717164449,721244161,
722923201,727083001,739444021,743404663,744866305,745864945,746706961,752102401,
759472561,765245881,771043201,775368901,775866001,776176261,781347841,784966297,
790020001,790623289,794937601,798770161,804978721,809702401,809883361,814056001,
822531841,824389441,824405041,829678141,832060801,833608321,834244501,834720601,
836515681,839275921,841340521,842202361,843704401,847491361,849064321,851703301,
851934601,852729121,854197345,855734401,860056705,863984881,867800701,868234081,
876850801,882796321,885336481,888700681,891706861,897880321,902645857,914801665,
918661501,928482241,930745621,931694401,934784929,935794081,939947009,940123801,
941056273,945959365,947993761,954732853,955134181,957044881,958735681,958762729,
958970545,962442001,962500561,963163201,963168193,968553181,968915521,975303121,
977737321,977892241,981567505,981789337,985052881,986088961,990893569,993420289,
993905641,1001152801,1018928485,1027334881,1030401901,1031750401,1035608041,
1038165961,1055384929,1070659201,1072570801,1074363265,1079556193,1090842145,
1093916341,1100674561,1103145121,1125038377,1131222841,1132988545,1134044821,
1136739745,1138049137,1140441121,1150270849,1152793621,1162202581,1163659861,
1177195201,1177800481,1180398961,1183104001,1189238401,1190790721,1193229577,
1194866101,1198650961,1200456577,1200778753,1206057601,1207252621,1210178305,
1213619761,1214703721,1216631521,1223475841,1227220801,1227280681,1232469001,
1251295501,1251992281,1254318481,1256855041,1257102001,1260332137,1264145401,
1268604001,1269295201,1271325841,1295577361,1299963601,1309440001,1312114945,
1312332001,1316958721,1317828601,1318126321,1321983937,1330655041,1332521065,
1337805505,1348964401,1349671681,1362132541,1376844481,1378483393,1382114881,
1384157161,1394746081,1394942473,1404111241,1404228421,1407548341,1410833281,
1419339691,1420379065,1422477001,1423668961,1428966001,1439328001,1439492041,
1441316269,1442761201,1445084173,1452767521,1481619601,1483199641,1490078305,
1490522545,1504651681,1505432881,1507746241,1515785041,1520467201,1521835381,
1528936501,1529544961,1534274841,1540454761,1545387481,1571503801,1573132561,
1574601601,1576826161,1583582113,1588247851,1592668441,1597009393,1597821121,
1615335085,1618206745,1626167341,1632785701,1641323905,1646426881,1648076041,
1657700353,1659935761,1672719217,1674309385,1676203201,1676641681,1678569121,
1685266561,1688214529,1688639041,1689411601,1690230241,1696572001,1698623641,
1699279441,1701016801,1705470481,1708549501,1726372441,1746692641,1750412161,
1752710401,1760460481,1769031901,1772267281,1773486001,1776450565,1778382541,
1784291041,1784975941,1785507361,1795216501,1801558201,1803278401,1817067169,
1825568641,1828377001,1831048561,1833328621,1836304561,1841034961,1845871105,
1846817281,1848112761,1848681121,1849811041,1854001513,1855100017,1858395529,
1879480513,1887933601,1894344001,1896961801,1899525601,1913016001,1916729101,
1918052065,1919767681,1932608161,1942608529,1943951041,1949646601,1950276565,
1954174465,1955324449,1958102641,1962804565,1976295241,1984089601,1988071801,
1992841201,1999004365,2000436751,2023528501,2029554241,2048443501,2048751901,
2049293401,2064236401,2064373921,2067887557,2073560401,2080544005,2097317377,
2101170097,2105594401,2107535221,2111416021,2111488561,2117725921,2126689501,
2140538401,2140699681,2159003281,2170282969,2176838049,2178944461,2199700321,
2199931651,2201169601,2212935985,2215407601,2216430721,2217951073,2223876601,
2224519921,2230305949,2232385345,2239622113,2244932281,2246916001,2258118721,
2265650401,2272748401,2278677961,2295209281,2301745249,2302419601,2309027281,
2313774001,2320224481,2320690177,2322648901,2323147201,2332627249,2335640077,
2339165521,2342644921,2353639681,2359686241,2361232477,2367379201,2377166401,
2391137281,2396357041,2407376665,2414829781,2430556381,2436691321,2438403661,
2443829641,2444950561,2456536681,2457411265,2467813621,2470348441,2470894273,
2479305985,2480343553,2489462641,2494465921,2494621585,2497638781,2509860961,
2510085721,2519819281,2523947041,2527812001,2529410281,2539024741,2544590161,
2547621973,2560104001,2560600351,2561945401,2564889601,2573686441,2574243721,
2575260241,2586927553,2588653081,2597928961,2598933481,2601144001,2602378721,
2605557781,2607162961,2607237361,2616662881,2617181281,2630374741,2642025673,
2657502001,2665141921,2677147201,2685422593,2690867401,2693939401,2702470861,
2709611521,2723859001,2733494401,2735309521,2766172501,2766901501,2770560241,
2776874941,2787998641,2797002901,2801124001,2806205689,2811315361,2815304401,
2832480001,2833846561,2842912381,2858298301,2867755969,2942952481,2943556201,
2965085641,2998467901,3001561441,3007991701,3022354401,3024774901,3025708561,
3030758401,3034203361,3035837161,3044238121,3044970001,3068534701,3069196417,
3072080089,3072094201,3077802001,3078386641,3086434561,3088134721,3090578401,
3102234751,3104207821,3105567361,3112974481,3119101921,3138302401,3159422785,
3159939601,3164207761,3180288385,3180632833,3188744065,3190894201,3193414093,
3203895601,3215031751,3222053185,3232450585,3240392401,3245477761,3246238801,
3248891101,3249390145,3263564305,3264820001,3270933121,3277595665,3281736601,
3284630713,3307322305,3313196881,3313744561,3314111761,3319323601,3328437481,
3345878017,3347570941,3348463105,3353809537,3378014641,3380740301,3411338491,
3413656441,3429457921,3438721441,3441837421,3480174001,3504570301,3508507801,
3521441665,3534510001,3555636481,3574014445,3575798785,3576804001,3600918181,
3618244081,3630291841,3637405045,3682471321,3697952401,3711456001,3712280041,
3713287801,3715938721,3722793481,3727589761,3754483201,3767865601,3776698801,
3787491457,3799111681,3800513761,3801823441,3805181281,3832646221,3834444901,
3835537861,3858853681,3863326897,3880251649,3891338101,3892568065,3901730401,
3907357441,3922752121,3928256641,3951813601,3981047941,3998554561,4015029061,
4030864201,4034969401,4059151489,4065133501,4077957961,4115677501,4127050621,
4134273793,4138747921,4146685921,4160472121,4162880401,4167038161,4169092201,
4169867689,4189909501,4199202001,4199529601,4199932801,4202009461,4210922233,
4212413569,4215885697,4216799521,4277982241 // 1119: 4295605861, 溢出
};
// hash表的改进方法, 只需要保存下标
// 经测试, 最多有6次冲突, b[i][0]保存数目
// 注意buckets的个维大小是人工测试得到!!
static short buckets[1021][6+1];
// 自动初始化hash表
static struct init
{
init(void)
{
for(int i = 0; i < NELEMS(CarmichaeNumbers); ++i)
{
// buckets[v]类似一个栈
// buckets[v][0]记录栈顶位置
// 最多有6次冲突, 不会溢出
int key = CarmichaeNumbers[i]%NELEMS(buckets);
assert(buckets[key][0] < NELEMS(buckets[0])-1);
buckets[key][++buckets[key][0]] = i;
}
}
} _init;
// bsearch的判断函数
static int cmp(const void *x, const void *y)
{
return *((unsigned*)x) - *((unsigned*)y);
}
// 判断是否是卡尔麦克数
static bool isCarmichaelNum(unsigned n)
{
// 改进: 采用hash表查找
int i, key = n%NELEMS(buckets);
for(i = 1; i <= buckets[key][0]; ++i)
if(n == CarmichaeNumbers[buckets[key][i]]) return true;
return false;
// 二分法查表CarmichaeNumbers表
return NULL !=
bsearch(&n, CarmichaeNumbers, NELEMS(CarmichaeNumbers), sizeof(n), cmp);
}
// 判断是否是费马(伪)素数
static bool isFermatNum(unsigned p)
{
unsigned m = p--;
unsigned r = 2%m;
unsigned k = 1;
// 蒙格马利快速幂模算法
while(p > 1)
{
if(p&1) k = (k*r) % m;
r = (r*r)%m; p >>= 1;
}
return (r*k)%m == 1;
}
// 根据费马小定理测试
bool isPrime(unsigned n)
{
// 处理比较特殊的情况
if(n < 2) return false;
if(n == 2) return true;
// 是费马(伪)素数但不是卡尔麦克数则必定为素数
return isFermatNum(n) && !isCarmichaelNum(n);
}
// 测试: hash的效果
#include <stdio.h>
main()
{
printf("num = %d\n", NELEMS(CarmichaeNumbers));
return 0;
}

复制代码

chai2010 · 发表于 2006-11-3 17:53:14

hash函数关键在于buckets[1021][6+1]
桶的大小是根据测试效果设置(最大冲突次数为6)

sizeof(bucket) + sizeof(CarmichaeNumbers) =

= sizeof(short)*1021*7 + sizeof(unsigned)*1118

= sizeof(int)*4691.5 > sizeof(int)*1024*4

chai2010 · 发表于 2006-11-6 09:26:24

makePrimes小的改进：
可以在计算商的时候得到余数，
因为很机器支持这种功能，
因此div函数可能比/%更高效些

// 构造素数序列primes[]
#include <stdlib.h>
void makePrimes(int primes[], int num)
{
primes[0] = 2;
for(int p = 3, cnt = 1; cnt < num; p += 2)
{
for(int k = 0; ; ++k)
{
// 可以在div时得到余数
div_t dt = div(p, primes[k]);
if(dt.rem == 0) break;
if(dt.quot <= primes[k]) { primes[cnt++] = p; break; }
}
}
}

复制代码

chai2010 · 发表于 2006-11-6 09:31:37

Post by chai2010

// 判断是否是费马(伪)素数

static bool isFermatNum(unsigned p)

{

unsigned m = p--;

unsigned r = 2%m;

unsigned k = 1;

// 蒙格马利快速幂模算法

while(p > 1)

{

if(p&1) k = (k*r) % m;

r = (r*r)%m; p >>= 1;

}

return (r*k)%m == 1;

}

复制代码

这里还有个问题：(r*r)%m可能会溢出！
解决的方法是采用64位整数%32位整数，

现在代码还有些问题，以后再给出

		自动登录	找回密码
密码			注册