/* 
  Compute prime numbers using Sieve of Eratosthenes 
  and find the smallest difference between two primes 
  whose sum is 1 billion.
*/

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

/* maximum size of primes and it's square root plus 1*/ 
#define MAXP 1000000000  
#define MAXPSQR   31623  

/* info about characters for my i386 machine         */
#define CHAR_BIT 8            /* bits in a character */     
  /* mask to turn on odd bits and turn off even bits */
#define CHAR_BIT_ODD_MASK 0xaa

#define MIDP  MAXP / 2   /* The middle of the list */

#define MAXC  MAXP / CHAR_BIT + 2  /* number of chars for bit map */

#define BITMASK(b) (1 << ((b) % CHAR_BIT))
#define BITSLOT(b) ((b) / CHAR_BIT)
#define BITSET(a, b)   ((a)[BITSLOT(b)] |= BITMASK(b))
#define BITCLEAR(a, b) ((a)[BITSLOT(b)] &= ~BITMASK(b))
#define BITTEST(a, b)  ((a)[BITSLOT(b)] & BITMASK(b))

int main() {
  int j,k;
  int jmin,jmax;

/* this array has one bit for each number up to MAXP */
  static unsigned char map[MAXC]; 

  printf("Find the smallest difference between" 
         " two primes whose sum is %d\n", MAXP);

/* clear map of primes. Assume odd numbers are prime except 2 */
  for (j=0; j<MAXC; j++) 
    map[j] = (unsigned char) CHAR_BIT_ODD_MASK;
  BITSET(map,2);

/* starting with 3 go through map of primes and remove  */
/* the multiples (non-primes) */
  for (j=3; j<=MAXPSQR; j+=2)   /* up to square root of maximum */
    if (BITTEST(map,j))         /* if the number is a prime    */
      for (k=j*j; k<MAXP; k+=j) /* start with the square */
        BITCLEAR(map,k);     

/* search for a primes that add up to MAXP starting at the midpoint */
  jmax = 0;                        /* solution prime */
  for (j=MIDP+1; j<MAXP; j+=2)     /*  search after the midpoint */
    if (BITTEST(map,j)) {
      if (BITTEST(map,(MAXP-j))) { 
        jmax = j;
        break;
      }
    }

/* show results */
  if (jmax) {
    jmin = MAXP - jmax;
    printf("min prime=%d\n",jmin);
    printf("max prime=%d\n",jmax);
    printf("\n");
    k = jmax - jmin;
    printf("difference=%d\n",k);
  }
  else
    printf("no two primes found that add up to %d\n", MAXP);
  printf("\n");

  return 0;
} /* main() */