# 【算法之高效求素数】浅析求素数算法

1. 根据概念判断:

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

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; }

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

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; }

4. 剔除因子中的重复判断.

I2. 如果d是素数, 则定理得证, 算法终止.
I3. 令n=d, 并转到步骤I1.

// 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; }

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

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

// 构造素数序列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的时间复杂度比较复杂, 而且它只有在初始化的时候才被调用一次.

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

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

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

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

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

// 这段代码可以由程序直接生成 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 };

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

7. 二分法查找

// 缺少的代码请参考前边 #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分法查找

---- (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. 打包成素数库(给出全部的代码)

// 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. 回顾, 以及推广

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 };