确定一个数是否为质数

我正在尝试提出一个方法,它接受一个整数,并返回一个布尔值来表示这个数字是否为质数,我对 C 知之甚少; 有人愿意给我一些指示吗?

基本上,我会这样用 C # 做:

static bool IsPrime(int number)
{
for (int i = 2; i < number; i++)
{
if (number % i == 0 && i != number)
return false;
}
return true;
}
203879 次浏览

好,那么忘掉 C 吧,假设我给你一个数,然后让你确定它是不是质数。你是怎么做到的?清楚地写下这些步骤,那么担心把它们翻译成代码。

一旦你确定了算法,你就可以更容易地找出如何编写程序,而且其他人也可以帮助你。

编辑: 这是你发布的 C # 代码:

static bool IsPrime(int number) {
for (int i = 2; i < number; i++) {
if (number % i == 0 && i != number) return false;
}
return true;
}

这是 非常接近有效的 C; C 中没有 bool类型,也没有 truefalse,所以需要稍微修改一下(编辑: Christopher Johnson 正确地指出 C99添加了 stdboool.h 头)。因为有些人不能访问 C99环境(但是您应该使用它!),让我们做一个非常小的改变:

int IsPrime(int number) {
int i;
for (i=2; i<number; i++) {
if (number % i == 0 && i != number) return 0;
}
return 1;
}

这是一个完全有效的 C 程序,可以做你想做的事情。我们可以不费吹灰之力地稍微改进一下。首先,请注意 i总是小于 number,因此 i != number总是成功的检查; 我们可以去掉它。

另外,您实际上不需要一直尝试到 number - 1的除数; 当您到达 sqrt (number)时,您可以停止检查。由于 sqrt是一个浮点运算,它带来了一大堆细微之处,因此我们实际上不会计算 sqrt(number)。相反,我们可以检查 i*i <= number:

int IsPrime(int number) {
int i;
for (i=2; i*i<=number; i++) {
if (number % i == 0) return 0;
}
return 1;
}

但是,还有最后一件事: 在您的原始算法中有一个小 bug!如果 number是负数,或者0,或者1,这个函数就会声明这个数是质数。您可能希望正确地处理这个问题,并且可能希望使 number为无符号的,因为您更可能只关心正值:

int IsPrime(unsigned int number) {
if (number <= 1) return 0; // zero and one are not prime
unsigned int i;
for (i=2; i*i<=number; i++) {
if (number % i == 0) return 0;
}
return 1;
}

这绝对不是检查一个数字是否为素数的最快方法,但它确实有效,而且非常简单。我们几乎不需要修改你的代码!

  1. 建立一个小素数表,并检查他们是否除以你的输入数。
  2. 如果数字存活到1,尝试使用基数增加的伪素性测试。
  3. 如果你的数字存活到2,你可以得出结论,它是素数,如果它低于一些众所周知的界。否则你的答案只会是“可能是质数”。您将在 wiki 页面中找到这些边界的一些值。

检查每个整数的模数,从2到你要检查的数字的根。

如果模等于零,那就不是质数。

伪代码:

bool IsPrime(int target)
{
for (i = 2; i <= root(target); i++)
{
if ((target mod i) == 0)
{
return false;
}
}


return true;
}

我很惊讶竟然没人提起这件事。

使用 埃拉托斯特尼筛法

详情:

  1. 基本上,除了1和它们自己之外,非素数还可以被另一个数整除
  2. 因此: 一个非素数将是素数的乘积。

埃拉托斯特尼筛法找到一个质数并存储起来。当检查一个新数是否为素数时,所有以前的素数都会对照已知素数列表进行检查。

原因:

  1. 这个算法/问题被称为“ 令人尴尬的平行
  2. 它创建了一个素数集合
  3. 这是一个动态规划问题的例子
  4. 很快的!

我要补充的是,偶数(第2栏)不能是质数。这将导致 for 循环之前的另一个条件。所以结束代码应该是这样的:

int IsPrime(unsigned int number) {
if (number <= 1) return 0; // zero and one are not prime
if ((number > 2) && ((number % 2) == 0)) return 0; //no even number is prime number (bar 2)
unsigned int i;
for (i=2; i*i<=number; i++) {
if (number % i == 0) return 0;
}
return 1;
}
int is_prime(int val)
{
int div,square;


if (val==2) return TRUE;    /* 2 is prime */
if ((val&1)==0) return FALSE;    /* any other even number is not */


div=3;
square=9;    /* 3*3 */
while (square<val)
{
if (val % div == 0) return FALSE;    /* evenly divisible */
div+=2;
square=div*div;
}
if (square==val) return FALSE;
return TRUE;
}

2和偶数的处理被排除在主循环之外,主循环只处理奇数除以奇数的问题。这是因为奇数模偶数总是给出一个非零的答案,这使得这些测试变得多余。或者,换句话说,一个奇数 可以被另一个奇数整除,但是 永远不会可以被一个偶数整除(E * E = > E,E * O = > E,O * E = > E 和 O * O = > O)。

在 x86架构上,除法/模数的开销确实很大,尽管开销有所不同(参见 http://gmplib.org/~tege/x86-timing.pdf)。另一方面,乘法相当便宜。

斯蒂芬 · 卡农回答得很好!

但是

  • 该算法可以通过观察除2和3以外的所有素数都是6k ± 1形式来进一步改进。
  • 这是因为对于某个整数 k 和 i =-1,0,1,2,3或4,所有整数都可以表示为(6k + i) ; 2除(6k + 0) ,(6k + 2) ,(6k + 4) ; 3除(6k + 3)。
  • 所以一个更有效的方法是检验 n 是否可以被2或3整除,然后检查所有6k ± 1≤√ n 的数。
  • 这是测试所有 m 直到√ n 速度的3倍。

    int IsPrime(unsigned int number) {
    if (number <= 3 && number > 1)
    return 1;            // as 2 and 3 are prime
    else if (number%2==0 || number%3==0)
    return 0;     // check if number is divisible by 2 or 3
    else {
    unsigned int i;
    for (i=5; i*i<=number; i+=6) {
    if (number % i == 0 || number%(i + 2) == 0)
    return 0;
    }
    return 1;
    }
    }
    

这个程序是非常有效的检查一个数字的素数检查。

bool check(int n){
if (n <= 3) {
return n > 1;
}


if (n % 2 == 0 || n % 3 == 0) {
return false;
}
int sq=sqrt(n); //include math.h or use i*i<n in for loop
for (int i = 5; i<=sq; i += 6) {
if (n % i == 0 || n % (i + 2) == 0) {
return false;
}
}


return true;
}

在阅读了这个问题之后,我被这样一个事实所吸引: 有些答案通过运行一个倍数为2 * 3 = 6的循环来提供优化。

所以我用同样的思想创建了一个新函数,但是它的倍数是2 * 3 * 5 = 30。

int check235(unsigned long n)
{
unsigned long sq, i;


if(n<=3||n==5)
return n>1;


if(n%2==0 || n%3==0 || n%5==0)
return 0;


if(n<=30)
return checkprime(n); /* use another simplified function */


sq=ceil(sqrt(n));
for(i=7; i<=sq; i+=30)
if (n%i==0 || n%(i+4)==0 || n%(i+6)==0 || n%(i+10)==0 || n%(i+12)==0
|| n%(i+16)==0 || n%(i+22)==0 || n%(i+24)==0)
return 0;


return 1;
}

通过同时运行函数和检查次数,我可以说明这个函数真的比较快。让我们看两个有两个不同素数的测试:

$ time ./testprimebool.x 18446744069414584321 0
f(2,3)
Yes, its prime.
real    0m14.090s
user    0m14.096s
sys     0m0.000s


$ time ./testprimebool.x 18446744069414584321 1
f(2,3,5)
Yes, its prime.
real    0m9.961s
user    0m9.964s
sys     0m0.000s


$ time ./testprimebool.x 18446744065119617029 0
f(2,3)
Yes, its prime.
real    0m13.990s
user    0m13.996s
sys     0m0.004s


$ time ./testprimebool.x 18446744065119617029 1
f(2,3,5)
Yes, its prime.
real    0m10.077s
user    0m10.068s
sys     0m0.004s

所以我想,如果推广的话,会不会有人获益太多?我想出了一个函数,它将首先清理一个给定的原始质数列表,然后使用这个列表来计算更大的一个。

int checkn(unsigned long n, unsigned long *p, unsigned long t)
{
unsigned long sq, i, j, qt=1, rt=0;
unsigned long *q, *r;


if(n<2)
return 0;


for(i=0; i<t; i++)
{
if(n%p[i]==0)
return 0;
qt*=p[i];
}
qt--;


if(n<=qt)
return checkprime(n); /* use another simplified function */


if((q=calloc(qt, sizeof(unsigned long)))==NULL)
{
perror("q=calloc()");
exit(1);
}
for(i=0; i<t; i++)
for(j=p[i]-2; j<qt; j+=p[i])
q[j]=1;


for(j=0; j<qt; j++)
if(q[j])
rt++;


rt=qt-rt;
if((r=malloc(sizeof(unsigned long)*rt))==NULL)
{
perror("r=malloc()");
exit(1);
}
i=0;
for(j=0; j<qt; j++)
if(!q[j])
r[i++]=j+1;


free(q);


sq=ceil(sqrt(n));
for(i=1; i<=sq; i+=qt+1)
{
if(i!=1 && n%i==0)
return 0;
for(j=0; j<rt; j++)
if(n%(i+r[j])==0)
return 0;
}
return 1;
}

我想我没有优化代码,但这是公平的。现在,测试。因为有这么多的动态内存,我希望列表235比硬编码的列表235慢一点。不过还可以,因为你可以看到下面。在那之后,时间变得越来越少,最终的最佳列表是:

235711131719

用了8.6秒。因此,如果有人要创建一个硬编码程序,使用这样的技术,我会建议使用清单23和5,因为增益不是很大。但是,如果愿意编写代码,这个列表也是可以的。问题是,如果没有循环,就不能说明所有的情况,否则代码会非常大(有1658879个 ORs,即各自内部 if中的 ||)。下一个清单:

23571113171923

时间变得越来越长,只有13秒,整个测试过程是这样的:

$ time ./testprimebool.x 18446744065119617029 2 3 5
f(2,3,5)
Yes, its prime.
real    0m12.668s
user    0m12.680s
sys     0m0.000s


$ time ./testprimebool.x 18446744065119617029 2 3 5 7
f(2,3,5,7)
Yes, its prime.
real    0m10.889s
user    0m10.900s
sys     0m0.000s


$ time ./testprimebool.x 18446744065119617029 2 3 5 7 11
f(2,3,5,7,11)
Yes, its prime.
real    0m10.021s
user    0m10.028s
sys     0m0.000s


$ time ./testprimebool.x 18446744065119617029 2 3 5 7 11 13
f(2,3,5,7,11,13)
Yes, its prime.
real    0m9.351s
user    0m9.356s
sys     0m0.004s


$ time ./testprimebool.x 18446744065119617029 2 3 5 7 11 13 17
f(2,3,5,7,11,13,17)
Yes, its prime.
real    0m8.802s
user    0m8.800s
sys     0m0.008s


$ time ./testprimebool.x 18446744065119617029 2 3 5 7 11 13 17 19
f(2,3,5,7,11,13,17,19)
Yes, its prime.
real    0m8.614s
user    0m8.564s
sys     0m0.052s


$ time ./testprimebool.x 18446744065119617029 2 3 5 7 11 13 17 19 23
f(2,3,5,7,11,13,17,19,23)
Yes, its prime.
real    0m13.013s
user    0m12.520s
sys     0m0.504s


$ time ./testprimebool.x 18446744065119617029 2 3 5 7 11 13 17 19 23 29
f(2,3,5,7,11,13,17,19,23,29)
q=calloc(): Cannot allocate memory

附言。我没有刻意释放(r) ,把这个任务交给操作系统,因为程序一退出,内存就会被释放,以获得一些时间。但是,如果您打算在计算之后继续运行代码,那么释放它将是明智的。


额外奖励

int check2357(unsigned long n)
{
unsigned long sq, i;


if(n<=3||n==5||n==7)
return n>1;


if(n%2==0 || n%3==0 || n%5==0 || n%7==0)
return 0;


if(n<=210)
return checkprime(n); /* use another simplified function */


sq=ceil(sqrt(n));
for(i=11; i<=sq; i+=210)
{
if(n%i==0 || n%(i+2)==0 || n%(i+6)==0 || n%(i+8)==0 || n%(i+12)==0 ||
n%(i+18)==0 || n%(i+20)==0 || n%(i+26)==0 || n%(i+30)==0 || n%(i+32)==0 ||
n%(i+36)==0 || n%(i+42)==0 || n%(i+48)==0 || n%(i+50)==0 || n%(i+56)==0 ||
n%(i+60)==0 || n%(i+62)==0 || n%(i+68)==0 || n%(i+72)==0 || n%(i+78)==0 ||
n%(i+86)==0 || n%(i+90)==0 || n%(i+92)==0 || n%(i+96)==0 || n%(i+98)==0 ||
n%(i+102)==0 || n%(i+110)==0 || n%(i+116)==0 || n%(i+120)==0 || n%(i+126)==0 ||
n%(i+128)==0 || n%(i+132)==0 || n%(i+138)==0 || n%(i+140)==0 || n%(i+146)==0 ||
n%(i+152)==0 || n%(i+156)==0 || n%(i+158)==0 || n%(i+162)==0 || n%(i+168)==0 ||
n%(i+170)==0 || n%(i+176)==0 || n%(i+180)==0 || n%(i+182)==0 || n%(i+186)==0 ||
n%(i+188)==0 || n%(i+198)==0)
return 0;
}
return 1;
}

时间:

$ time ./testprimebool.x 18446744065119617029 7
h(2,3,5,7)
Yes, its prime.
real    0m9.123s
user    0m9.132s
sys     0m0.000s

使用埃拉托斯特尼筛法,计算速度比“已知范围”的素数算法要快得多。

通过使用来自它的 wiki (https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes)的伪代码,我能够在 C # 上得到解决方案。

public bool IsPrimeNumber(int val) {
// Using Sieve of Eratosthenes.
if (val < 2)
{
return false;
}


// Reserve place for val + 1 and set with true.
var mark = new bool[val + 1];
for(var i = 2; i <= val; i++)
{
mark[i] = true;
}


// Iterate from 2 ... sqrt(val).
for (var i = 2; i <= Math.Sqrt(val); i++)
{
if (mark[i])
{
// Cross out every i-th number in the places after i (all the multiples of i).
for (var j = (i * i); j <= val; j += i)
{
mark[j] = false;
}
}
}


return mark[val];
}

IsPrimeNumber (100000000)需要21秒758ms。

注意 : 值可能因硬件规格而异。

避免溢出错误

unsigned i, number;
...
for (i=2; i*i<=number; i++) {  // Buggy
for (i=2; i*i<=number; i += 2) {  // Buggy
// or
for (i=5; i*i<=number; i+=6) { // Buggy

number是素数且 i*i接近类型的最大值时,这些表单是不正确的。

问题存在于所有的整数类型,signed, unsigned和更广泛的。

例如:

UINT_MAX_SQRT作为最大整数值的平方根的底数。例如,当 unsigned为32位时,65535。

对于 for (i=2; i*i<=number; i++),发生这个10年之久的故障是因为当 UINT_MAX_SQRT*UINT_MAX_SQRT <= numbernumber是素数时,下一次迭代将导致乘法溢出。如果类型是 签了类型,则溢出为 UB。对于 未签名类型,这本身不是 UB,但逻辑已经崩溃。反复直到 被截断了产品超过 number。可能会出现不正确的结果。使用32位 unsigned,尝试使用质数4,294,967,291。

如果 some_integer_type_MAX梅森至尊,那么 i*i<=number就是 永远不会


为了避免这个错误,考虑到 number%inumber/i在许多编译器上都是高效的,因为商和余数的计算是一起完成的,因此两者兼顾不会产生额外的成本,而只需要1。

一个简单的全方位解决方案:

bool IsPrime(unsigned number) {
for(unsigned i = 2; i <= number/i; i++){
if(number % i == 0){
return false;
}
}
return number >= 2;
}

为了检查一个数字是否为质数,我使用了 米勒/拉宾算法。

#include <stdlib.h>


typedef size_t positive_number; // also try __uint128_t


static inline positive_number multiplication_modulo(positive_number lhs, positive_number rhs, positive_number mod) {
positive_number res = 0; // we avoid overflow in modular multiplication
for (lhs %= mod, rhs %= mod; rhs; (rhs & 1) ? (res = (res + lhs) % mod) : 0, lhs = (lhs << 1) % mod, rhs >>= 1);
return res; // <= (lhs * rhs) % mod
}


static int is_prime(positive_number n, int k) {
positive_number a = 0, b, c, d, e, f, g; int h, i;
if ((n == 1) == (n & 1)) return n == 2;
if (n < 51529) // fast constexpr check for small primes (removable)
return (n & 1) & ((n < 6) * 42 + 0x208A2882) >> n % 30 && (n < 49 || (n % 7 && n % 11 && n % 13 && n % 17 && n % 19 && n % 23 && n % 29 && n % 31 && n % 37 && (n < 1369 || (n % 41 && n % 43 && n % 47 && n % 53 && n % 59 && n % 61 && n % 67 && n % 71 && n % 73 && ( n < 6241 || (n % 79 && n % 83 && n % 89 && n % 97 && n % 101 && n % 103 && n % 107 && n % 109 && n % 113 && ( n < 16129 || (n % 127 && n % 131 && n % 137 && n % 139 && n % 149 && n % 151 && n % 157 && n % 163 && n % 167 && ( n < 29929 || (n % 173 && n % 179 && n % 181 && n % 191 && n % 193 && n % 197 && n % 199 && n % 211 && n % 223))))))))));
for (b = c = n - 1, h = 0; !(b & 1); b >>= 1, ++h);
for (; k--;) {
for (g = 0; g < sizeof(positive_number); ((char*)&a)[g++] = rand()); // random number
do for (d = e = 1 + a % c, f = n; (d %= f) && (f %= d););
while (d > 1 && f > 1);
for (d = f = 1; f <= b; f <<= 1);
for (; f >>= 1; d = multiplication_modulo(d, d, n), f & b && (d = multiplication_modulo(e, d, n)));
if (d == 1) continue;
for (i = h; i-- && d != c; d = multiplication_modulo(d, d, n));
if (d != c) return 0;
}
return 1;
}

测试是打印第一个质数:

#include <stdio.h>
int main() {
// C Fast Iterative Algorithm
// The First 10,000 Primes
for (int i = 0 ; i < 104730 ; ++i)
if (is_prime(i, 20))
printf("%d %c", i, (i+1) % 10 ? ' ' : '\n');


if (is_prime(9223372036854775783UL, 12))
if (is_prime(9223372036854775643UL, 12))
if (!is_prime(3037000493ULL * 3037000453ULL, 12))
printf("Done.\n");
}

您可以将其放入 primes.c文件,然后编译 + 执行:

gcc -O3 -std=c99 -Wall -pedantic primes.c ; ./a.out ;

费马测试的这个变体在 log(n)中有一个多项式运行时。
_ _ uint128 _ t类型可以与 128位整数 GCC 扩展一起使用。