Primality test using Fermat little theorem

Question

code:

void prime()    
{    
    int i,N;    
    scanf("%d",&N);    
    for(i=2;i<N;i++)            
    {    
        if (((i^(N-1))%N )==1);     
        else{    
            printf("not prime");   
            return;
        }     
    }    
    printf("prime");    
    return;    
}    

This program is based on Fermat's Theorem on prime numbers. N is number to be tested as prime. This program is not showing correct result for '11'. Maybe due to some mistake which is not identified by me.

Answer 1

You are running into overflow if this is pseudo-code or
If C code, use of ^ as power operator is not valid.

Working with large integers quickly becomes a problem in C. The are various BigInt libraries available.

Using floating point is challenging with large integer computation. Recommend avoiding double, pow(), etc.

Since the problem is all >= 0, suggest using unsigned integers. Also use the largest integer type available - typically unsigned long long. As overflow is a real possibility, detect it.

unsigned long long upower(unsigned i, unsigned N) {
  unsigned long long power = 1;
  if (i <= 1) return i;
  while (N-- > 0) {
    unsigned long long power_before = power;
    power *= i;
    if (power < power_before) {
      printf("Overflow\n");
      return 0;
    }
  }
  return power;
}

void prime() {
  unsigned i, N;
  scanf("%u", &N);
  for (i = 2; i < N; i++) {
    if ((upower(i, N - 1) % N) != 1) {
      printf("not prime");
      return;
    }
  }
  printf("prime");
  return;
}

In lieu of huge integers, the Chinese remainder theorem may offer an alternative to (upower(i, N - 1) % N) != 1.

Answer 2

If I read your code as pseudo-code, You're overflowing.

10^10 is bigger that 2^31 -1 which is the max value for most int. You could solve this for N=11 by using longs, but that will not get you far, you'll start overflowing at some point as well.

That theorem, at least expressed like this, is very unpractical to use with finite length numbers.

Now, if your code is real C, note that ^ means XOR, not exponentiation. Exponentiation is pow(). Thanks to the commenters for pointing that out.

Answer 3

Modular mathematical rules and principles can be applied here to show that in order to compute

(i ^ (N-1)) % N,

you do not even need to compute i^(N-1) at the first place. You can easily break down (N-1) into powers of 2. Let's take an example to make it more clear.

Assume that the subject of our primality test, N = 58.

So,

N - 1 = 57

57 can be easily rewritten as:

57 = 1 + 8 + 16 + 32

or,

57 = 2^0 + 2^3 + 2^4 + 2^5

So, substituting this value for N-1, we need to compute

(i ^ (2^0 + 2^3 + 2^4 + 2^5))% 58

or,

((i^1) × (i^8) × (i^16) × (i^32))% 58

Which, using the Modular Multiplication identities, can be rewritten as:

((i^1)% 58 × (i^8)% 58 × (i^16)% 58 × (i^32)% 58) mod 58 ---(1)

Note that,

(i^1)% 58 = i%58

can be easily computed without worrying of any overflows.

Once again utilising the Modular Multiplication identities, we know that

(i^2)% 58 = ((i^1)% 58 × (i^1)% 58)% 58

Substitute the value of (i^1)% 58 to find (i^2)% 58.

You can continue in this fashion, computing (i^4)% 58 through (i^32)% 58. Once completed, you can finally substitute the values in (1) to finally find the required value, very efficiently avoiding any overflows.

---

Note that, other modular exponientation techniques exist too. This was just an example to showcase how modular mathematical techniques can be used in implementing Fermat's little primality test.

Cheers!

Answer 4

Sorry to change your code a little. Using the BigInteger class, you can calculate very quickly for much larger numbers. However, you can use this method not to get prime numbers in order, but to test if any numbers are prime.

using System;
using System.Numerics;
                    
public class Program
{
    public static void Main()
    {
        Console.WriteLine(2);
        for(var i = 3; i < 100000; i+=2) 
        {
            if(BigInteger.ModPow(2, i , i) == 2)
                Console.WriteLine(i);
        }
    }
}

https://dotnetfiddle.net/nwDP7h

This code will produce erroneous results when it falls into the following numbers.

https://oeis.org/a001567 https://oeis.org/a006935

To fix these errors, you need to edit the code as follows and make a binary search within these numbers to test whether the number is a pseudo prime.

public static bool IsPrime(ulong number)
{
    return number == 2
        ? true
        : (BigInterger.ModPow(2, number, number) == 2
            ? (number & 1 != 0 && BinarySearchInA001567(number) == false)
            : false)
}

public static bool BinarySearchInA001567(ulong number)
{
    // Is number in list?
    // todo: Binary Search in A001567 (https://oeis.org/A001567) below 2 ^ 64
    // Only 2.35 Gigabytes as a text file http://www.cecm.sfu.ca/Pseudoprimes/index-2-to-64.html
}