Programs with Proofs: Finding if an Integer is Prime


/* Method is_prime() finds whether the input integer n (n >= 2) is a composite
   (returns 0) or a prime (returns 1). It uses a simple algorithm, which is not
   efficient for large n.

   To know if n is a composite, we can check sequentially if there is any divisor
   of n in S = {2, 3, ..., n-1}. If any integer i > sqroot(n) (sqroot(n) is the
   square-root of n) divides n, there must exist another integer j = n/i which
   also divides n. But we must have j < sqroot(n). So, it will suffice to check
   the integers in S only upto sqroot(n).

   Further, if integer i in S divides n and i is a composite, there must exist a
   prime p in S such that p < i and p divides n. So, to find if n is a composite,
   it will suffice to check if any prime in S divides n. But we don't already know
   which integers in S are prime.

   We can make use of the fact that all prime numbers greater than 3 are of the
   form (6k+1) or (6k-1), where k = 1, 2, 3,... . For the proof, you may check
   https://primes.utm.edu/notes/faq/six.html.

   (Note that the inverse is not necessarily true; e.g. 25 and 35 are not prime).

   Say S1 = {2, 3}, S2 = {integers (6k+1)/(6k-1) upto sqroot(n); k = 1, 2, 3, ...}.
   So, to find if n is a composite, we can check if the set (S1 U S2) (considering
   integers only upto sqroot(n)) contains any divisor of n.
 */

/* expects (n >= 2) */
int is_prime(int n)
{
  int d;

  if(n == 2 || n == 3)
    return 1;

  /* (n >= 4) holds */

  /* check divisors of n in S1 (since n >= 4 here) */
  if(n % 2 == 0 || n % 3 == 0)
    return 0;

  /* R1: (No divisors of n in S1) holds */

  d = 5;
  /* k = 1; (k is just a representative integer variable) */
  /* (d = 6k-1) holds */

  /* loop invariants:
       P: (d = 6k-1)
       Q: (No divisors of n in S2 upto 6(k-1)+1) */

  while((long)d*d <= n)
  {
    /* since d >= 5 here, (d+2 < d*d); so we must always have (d+2) < n */

    /* check whether 6k-1 (d) or 6k+1 (d+2) divides n */
    if(n % d == 0 || n % (d+2) == 0)
      return 0;

    d = d + 6;
    /* k = k + 1; */
  }

  /* R2: (Q applied for the complete S2, i.e. No divisors of n in S2) holds */

  /* R1 and R2 imply n is prime */

  return 1;
}