c - C 中的费马素数测试失败

Question

#include <stdio.h>
#include <time.h>
#include <math.h>
#include <stdlib.h> 

#define MAXNUM 2000000000
#define MINNUM 1990000001
#define MAXTRIES 10

unsigned long long b, e, m, result; 

int modulo(b, e, m) 
{
    result = 1;

    while(e > 0)
    {
        if(e % 2 == 1)
        {
            result = (result * b);
        } 

        b = (b * b) % m;
        e = e / 2;
    }

    return result % m;
}

int isPrime(n) 
{
    unsigned long long a; 

    int i; 

    for(i = 1; i <= 10; i++)
    {
        a = rand() % (n - 1) + 1;
        if(modulo(a, n - 1, n) != 1)
        {
            return 0;
        }
    }

    return 1;
}

int main() 
{
    unsigned int prime = 0;
    unsigned int flag = 0;
    unsigned int tries;
    unsigned int start;
    long curtime;
    unsigned long long p;

    curtime = time(NULL);
    srand((unsigned int) curtime);
    printf("Checking range [1990000001, 2000000000] for prime numbers.\n");
    if(MINNUM % 2 == 0)
    {            
        start = MINNUM + 1;      
    }
    else
    {
        start = MINNUM;    
    }

    printf("Trying Fermat test with seed %ld \n\n",curtime);
    prime = 0;

    for(tries = 1; tries <= MAXTRIES; tries++)
    {
        clock_t tic = clock();
        for(p = start; p <= MAXNUM; p += 2)
        {
            if(isPrime(p))
                prime++;
        } 

        clock_t toc = clock();
        printf("Probabilistic algorithm: Found %ld primes in %f seconds.(tries = %d)\n", prime, (double)(toc - tic) / CLOCKS_PER_SEC,tries);
        prime = 0;
    } 

    return 0;
}

所以问题是算法在每次尝试中找到 5000000 个素数，而它应该找到大约 466646 个，但有一些偏差。这意味着在每次尝试中它都应该找到接近上述数量的素数。

Answer 1

看起来主要问题是由 modulo() 函数中的整数溢出引起的。具体来说，result=(result*b) 会经常溢出。你需要将这些变量存储在64位无符号整数中，并且每次都计算这个结果的模数。

这将起作用(在其他地方进行一些小的更正):

#include <inttypes.h>

#define MAXNUM 2000000000
#define MINNUM 1990000001
#define MAXTRIES 10


uint64_t modulo(uint64_t b, uint64_t e, uint64_t m){
    uint64_t result=1;
    while(e>0){
        if(e%2==1){
            result=(result*b)%m;
        }
        b=(b*b)%m;
        e=e/2;
    }
    return result%m;
}

结果:

Checking range [1990000001, 2000000000] for prime numbers.
Trying Fermat test with seed 1416322197 

Probabilistic algorithm: Found 466646 primes in 5.157485 seconds.(tries=1)