CODEWITHSUNDEEP
ctestingprimes
CharisAlex
CharisAlex

Reputation: 91

Fermat Primality Test failure in C

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

So the problem is that the algorithm finds in every try 5000000 prime numbers when it should find around 466646 with some deviation. Which means that in every try it should find a number of primes close to the one mentioned above.

Upvotes: 3

Views: 397

Answers (1)

r3mainer
r3mainer

Reputation: 24587

It looks like the main problem is caused by integer overflows in the modulo() function. Specifically, result=(result*b) is going to overflow quite regularly. You need to store these variables in 64-bit unsigned integers, and calculate the modulus of this result every time.

This will work (with a few minor corrections elsewhere):

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

Result:

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

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

Upvotes: 2

Related Questions