素数判断

版本1

素数，即除1以外，只能被1和其本身整除的数。因此根据定义，最初的判断素数的方法为版本1。

//judge prime

//program 1

bool prime_1(int n)

{

int i;

for(i=2;i

{

if(n%i==0)

{

return false;

}

}

return true;

}

版本2

由于一个数的因式分解可表示为x = a * b, a<=b，因此在判断整数n是不是素数时并不需要整数从2 ～ n，只需要整除2 ～ sqrt(n) 即可。

//judge prime

//program 2

bool prime_2(int n)

{

int i,n2;

n2=sqrt(n);

for(i=2;i<=n2;i++)

{

if(n%i==0)

{

return false;

}

}

return true;

}

但是需要注意的是，int类型只能表示2^-31 ~ 2^31 范围，并不能完整表示10位大小的整数n，因此需要使用能够long long int类型。其次，求平方根的函数原型为double sqrt(double)，进行sqrt(n)计算时并不确定能保证其精度。最后，对e小数点后，每个10位数字大小的整数n进行判断，需要对每个2 ～ sqrt(n)的整数进行判断，单个判断耗时间太久，所有判断重复运算太多。因此次方法必然不行。

版本3

任何一个合数都可以表现为适当个素数的乘积的形式，所以我们只用素数去除要判断的数即可，比如要判断100以内的素数，只用10以内的2，3，5，7就够了。那么首先就要构建一个素数表，以空间换时间。

可以使用一种被称为“埃拉托色尼筛”的算法。该算法一开始初始化一个2～n的连续整数序列。从2开始，2的倍数肯定不是素数，去除，剩下的数中下一个是3。再把3的倍数去除，再下一个是5(4已经去除了)，以此类推，最后剩下的数必然是素数，因为它没有被筛，不是任何一个数的倍数(除了1和自己)。这样只需要很多次加法就可以了。

//create prime less then n

//program 3

#include

#include

#include //Linux, gcc -lm

#include

#define Long64 long long

#define MAX_N 8000000

Long64 prime_table[MAX_N/2]={0};

void create_prime_table(Long64 n)

{

Long64 i,j,n2;

bool prime[MAX_N+1]={0};

n2=sqrt(n);

i=2;

for(i=2;i<=n2;i++)

{

if(prime[i]==0)

{

for(j=i+i;j<=n;j+=i)

{

prime[j]=1;

}

}

}

j=0;

for(i=2;i<=n;i++)

{

if(prime[i]==0)

{

prime_table[j]=i;

j++;

printf("%lld : %lld\n",j,i);

}

}

}

void main()

{

create_prime_table(MAX_N);

}

这算是基本的筛法了，但它还有很多优化的空间。首先内存中存的表不只是素数，还有合数的空间，合数比素数多得多。 其次对于合数可能会被筛选多次，比如6，在2的时候筛掉一次。在3的时候又筛掉一次。

版本4

采用线性筛选方式，即对一个合数，只进行一次筛选，极大程度的减少重复计算。 同时，采用堆结构存放素数表。

//create prime less then n

//program 4

#include

#include

#include //Linux, gcc -lm

#include

#define Long64 long long

#define MAX_N 1000000000

void create_prime_table2(Long64 n)

{

Long64 i,j,k;

bool *prime;

Long64 *prime_table;

prime=(bool *)calloc(MAX_N+1,sizeof(bool));

prime_table=(Long64 *)calloc(MAX_N/100,sizeof(Long64));

i=2;

k=0;

for(i=2;i<=n;i++)

{

if(prime[i]==0)

{

prime_table[k++]=i;

//printf("%lld : %lld \n",k,i);

}

for(j=0;j

{

prime[prime_table[j]*i]=1;

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

break;

}

}

printf("k = %lld\n",k);

if(prime_table)

free(prime_table);

if(prime)

free(prime);

}

void main()

{

create_prime_table2(MAX_N);

}

这样做之后，程序运行速度快多了，但是只能计算到1x10^9....再大时，运行就提示段错误了。还要再在其他地方想办法。。。 暂时先记录到这里。

自然数e

在上上一篇博客中，虽然实现了计算n为自然数e的程序，但是GetLength 函数实现的数学原理一直没有想清楚。其实是这样的。 计算自然数e时，我们使用的时泰勒公式

e^x = 1 + x + x^2/2! +... + x^n/n!...

当考虑到泰勒公式的余项时，余项可以表示为

R(x) = (e^esp * x^(n+1))/(n+1)!

其中esp在0～x之间。再取x=1,得到

R(1) = e^esp/(n+1)! <= e/(n+1)! < 10^(-M)

10^(-M)就是我们要达到的误差限，也就是小数点后的后M位。 因此，若要计算到10000位的精度，只需计算(n+1)!/e > 10^10000即可。 程序中GetLength的实现原理就是这样。

素数判断

在上一篇博客中，只计算到1x10^9内的素数，并没有达到要求，主要时卡在了数据结构上，再大的话在我2G内存32位系统上就已经提示段错误了。现在回过头了想想，其实还可以在算法上再改进。 我们知道a*b = c，一个M位的数乘上另外一个M位的数，得到结果位数范围在M+1~2M的范围内。我们只要求出5位数以内的素数，就可以用来判断10位数的整数是否为素数了。 那么素数判断的程序就可以修改成这样，也就是版本2和版本4的结合。

版本5

#include

#include

#include //Linux, gcc -lm

#include

#define Long64 long long

#define MAX_N 100000

//create prime less then n

//program 5

Long64 create_prime_table3(bool *prime, Long64 *prime_table, Long64 n)

{

Long64 i,j,k;

i=2;

k=0;

for(i=2;i<=n;i++)

{

if(prime[i]==0)

{

prime_table[k++]=i;

//printf("%lld : %lld \n",k,i);

}

for(j=0;j

{

prime[prime_table[j]*i]=1;

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

break;

}

}

return k;

}

bool prime3(Long64 prime_table[], Long64 length, Long64 n)

{

Long64 i;

for(i=0;i

{

if(n%(prime_table[i])==0)

{

return false;

}

}

return true;

}

void main()

{

bool *prime;

Long64 *prime_table,i=0,j=0,length=0;

prime=(bool *)calloc(MAX_N+1,sizeof(bool));

prime_table=(Long64 *)calloc(MAX_N,sizeof(Long64));

length = create_prime_table3(prime,prime_table,MAX_N);

for(i=1000000000; i<10000000000; i++)

if(prime3(prime_table,length,i)!=false)

printf("%lld : %lld\n",j++,i);

}

收尾：判断e中首个连续10位素数

有了以上问题解决后，这道题终于可以解决了。不仅可以找到第一个10位长度的素数，还可以找到前10个！前10个结果如下，格式为“第几个：在e中的起始位置：要找的素数”

1:99:7427466391

2:123:7413596629

3:149:6059563073

4:171:3490763233

5:182:2988075319

6:201:1573834187

7:214:7021540891

8:218:5408914993

9:254:6480016847

10:295:9920695517

完整实现代码如下

#include

#include

#include //Linux, gcc -lm

#include

#include

#define MAX_E_BITE_N 1000

#define PRIME_N 1000000

#define Long64 long long

typedef unsigned int UINT;

int get_length(int nbite)

{

int m = 1;

double x = 2.0/2.72;

while(nbite)

{

while(x<1.0)

{

m++;

x *= (double)(m+1);

}

while(x >= 1.0 && nbite)

{

x /= 10.0;

--nbite;

}

}

while(x < 1.0)

{

m++;

x *= (double)(m+1);

}

return m;

}

int create_e(Long64 *euler, const UINT base, int m)

{

UINT *r = NULL;

UINT y = 0;

int i = 0,j = 0,k = 0;

r=(UINT *)calloc(m+1,sizeof(UINT));

assert(r != NULL);

assert(euler != NULL);

for(i = 0;i <= m; ++i)

{

r[i] = 1;

euler[i] = 0;

}

j = 1;

euler[0] = 2;

for(k = MAX_E_BITE_N; k > 0; k -= 5)

{

y = 0;

for(i = m; i >= 2; --i)

{

y = y + r[i]*base;

r[i] = y%i;

y /= i;

}

if(k < 5)

{

euler[j++] = y%base;

}

else

{

if(y < base)

euler[j++] = y;

else

{

if(r)

free(r);

return 0;

}

}

}

if(r)

free(r);

printf("\n\nE=%lld.\n",euler[0]);

for(i = 1; i < j; i++)

{

printf("%.5lld ",euler[i]);

}

printf("\n\neulur number is ok j=%d\n",j);

return j;

}

Long64 create_prime_table(bool *prime, Long64 *prime_table, Long64 n)

{

Long64 i,j,k;

i = 2;

k = 0;

for(i = 2;i <= n;i++)

{

if(prime[i] == 0)

{

prime_table[k++] = i;

// printf("%lld : %lld \n",k,i);

}

for(j = 0; j < k && (prime_table[j]*i) <= n; j++)

{

prime[prime_table[j] * i] = 1;

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

break;

}

}

return k;

}

bool is_prime(Long64 prime_table[], Long64 length, Long64 n)

{

Long64 i;

for(i = 0; i < length; i++)

{

if(n%(prime_table[i]) == 0)

{

return false;

}

}

return true;

}

void main()

{

const UINT base = 100000;

int m = 0,elength = 0,j = 0,find_number = 10;

Long64 *euler = NULL;

bool *prime = NULL;

Long64 *prime_table = NULL;

Long64 i = 0,length = 0,data = 0;

Long64 ten[10] = {1,10,100,1000,10000,100000,1000000,10000000,100000000,1000000000};

m = get_length(MAX_E_BITE_N);

euler = (Long64 *)calloc(MAX_E_BITE_N,sizeof(Long64));

elength = create_e(euler,base,m);

prime = (bool *)calloc(PRIME_N+1,sizeof(bool));

prime_table = (Long64 *)calloc(PRIME_N,sizeof(Long64));

length = create_prime_table(prime,prime_table,PRIME_N);

if(prime)

free(prime);

m = 0;

for(i = 1; (i+2) < elength; i++)

{

for(j = 0; j < 5; j++)

{

data = euler[i]%ten[5-j]*ten[5+j] + euler[i+1] * ten[j] + euler[i+2]/ten[5-j];

// printf("data : %lld \n",data);

if(is_prime(prime_table,length,data) == true)

{

printf("find it! %d:%lld:%lld\n",++m,5*(i-1)+j+1,data);

if(m == find_number)

break;

}

}

if(m == find_number)

break;

}

if(euler)

free(euler);

if(prime_table)

free(prime_table);

return;

}