/******************************************************* 
/ 
/ fermat-little-theorem.c 
/ 
/ A C program that uses Fermat's Little Theorem 
/ to test for probable primes 
/ 
/ Paul Soper 
/ 
/ June, 2014 
/ 
/******************************************************/ 
 
/****************************************************** 
/ 
/ Fermat's little theorem, stated in C-type language,  
/ says that if p is prime, and 1 < a < p,  
/ ((a^p - a) % p) == 0). 
/ 
/ If ((a^p - a) % p) != 0), p is definitely composite, 
/ a is called the composite witness. 
/ 
/ If ((a^p - a) % p) == 0), p is probably prime.  The 
/ probability of p being prime increases as more 
/ values of a are tested and are found not to be 
/ composite witnesses. 
/ 
/ If ((a^p - a) % p) == 0), but p is not prime, a is 
/ called a Fermat liar. 
/ 
/ There are a small group of numbers - the Carmichael 
/ numbers, for which all a's 1 < a < p are Fermat 
/ liars. 561 is such a number, and the only one that
/ will be found by this program.
/ 
/*****************************************************/  
 
 
#include "stdio.h" 
#include "stdlib.h" /* because we're using malloc */ 
#include "time.h"   /* because we need to seed srand() */ 
 
/* a^p - a may be very, very big, so our 64 bit  
   arithmetic will likely not suffice.  We will use the 
   GNU large integer library, gmp */ 
 
#include "gmp.h"     
 
int main() { 
  /* seed the random number generator with the current 
     time */ 
 
  srand(time(NULL)); 
 
  /* we will look for primes less than 1000 */ 
 
  int max = 1000; 
  unsigned int p; 
 
  mpz_t a;  /* big number variable for a */ 
  mpz_t e;  /* big number variable for a^p - a */ 
  mpz_t r1; /* big number variable for (a^p - a) % p */ 
 
  /* initialize the large integer variables */ 
 
  mpz_init (a); 
  mpz_init (e); 
  mpz_init (r1); 
 
   /* run through the candidates, from 3 to max, looking 
     for probable primes */ 
 
  for (p = 3; p < max; ++p) { 
    /* For each candidate, we're going to test 20 a's,  
       counting them with trial */ 
 
    int trial = 0;   
 
    /* if p is found to be composite, that is if we  
       find an a such that ((a^p - a) mod p != 0), we 
       want to report that a as our composite witness */ 
 
    int composite_witness = 0; 
 
    /* if we find a fermat liar, report that as well */ 
 
    int fermat_liar = 0; 
 
    /* a variable for a */ 
 
    unsigned int a1 = 0; 
 
    while ((trial < 20) && (composite_witness == 0)) { 
      /* in case we discover that the old a1 was a 
	 fermat liar, keep track of the old a1 */ 
 
      unsigned int olda1 = a1; 
 
      /* choose a random integer in the range 1 < a1 < p 
	 to be our a for the test */ 
 
      a1 = rand() % (p-2) + 2; 
 
      /* put a1 into the bug number variable a */ 
 
      mpz_set_ui(a, a1); 
 
      /* raise a to the power p, and put the  
	 result in e */ 
 
      mpz_pow_ui(e, a, p); 
 
      /* now take that result, and subtract a from it, 
	 putting the new result back into e */ 
 
      mpz_sub(e, e, a); 
 
      /* get the remainder from dividing e by p */ 
 
      mpz_mod_ui (r1, e, p); 
 
      /* compare that big number remainder to 0 */ 
 
      unsigned int cmp1 = mpz_cmp_ui (r1, 0); 
 
      /* if the remainder is not 0, p is composite */ 
 
      if (cmp1 != 0) { 
	/* identify the witness that showed p to be 
	   composite, and store it in composite_witness -  
	   this will also end the loop, as no further 
	   tests are needed */ 
	
	composite_witness = a1; 
  
	/* but, if this is not the first pass through the  
	   loop, a previous a must have lied about the 
	   primality of p rather than given witness to 
	   the composite nature of p, so the previous a 
	   is a fermat liar */ 
     
	fermat_liar = olda1; 
      } 
      ++trial; 
    } 
 
    if (composite_witness == 0) { 
      printf("%d is a probable prime\n", p); 
    }  
    else { 
      printf("%d is composite - %d is a composite witness", 
	     p, composite_witness); 
      if (fermat_liar > 0) { 
	printf(" - %d is a fermat liar for %d",  
	       fermat_liar, p); 
      } 
      printf("\n"); 
    } 
  } 
 
  /* clear the big number variables */ 
 
  mpz_clear (r1); 
  mpz_clear (a); 
  mpz_clear (e); 
 
  return (0); 
}