/* prime.cpp : Search for 11 digit primes 101234xyzuv, where the
 *   five letters are a permutation on {5, 6, 7, 8, 9}. These are
 *   the smallest 11-digit primes containing all 10 digits [there
 *   are no 10-digit primes containing all 10 digits by the rule
 *   of 9 : 9 divides a number iff 9 divides the sum of its digits].
 *
 *   The program's output follows [where pi(x) = #primes <= x]:
 *
 *          pi(101000) = 9673
 *          1 : 10123485679 is prime
 *          2 : 10123465789 is prime
 *          3 : 10123457689 is prime
 *          4 : 10123496857 is prime
 *          5 : 10123485769 is prime
 *          6 : 10123465897 is prime
 *
 *   The smallest prime including all 10 digits is 10 123 457 689.
 */
#include "stdio.h"
#include "math.h"     // sqrt
#include "stdlib.h"   // malloc, atoi

typedef int BOOL;

#define FALSE     0
#define TRUE      1
#define N    101000

BOOL CalculatePermutations(char letters[], int n, char **ppPerms);
void ShiftPerm(int x[], int j);
void LoadCurPerm(int x[], int n, char *pCurPerm, char letters[]);
BOOL CalculatePrimes(int n, int **primeList, int *numPrimes);
BOOL IsPrime(__int64 m, int primes[], int primeLimit);

void main(void)
{
  int *pPrimes;    // ptr to prime table
  char *pPerms;    // ptr to permutation strings
  int numPrimes;   // number of primes in table
  char letters[5]; // symbols to be permuted
  __int64 x;       // candidate for primality (11 digits)
  int count;       // to count the special primes found
  int i;

  // Init to NULL; will be non-NULL iff functions loaded them.
  pPrimes = NULL;
  pPerms = NULL;

  // Calculate primes <= N, exiting in case of failure.
  if (!CalculatePrimes(N, &pPrimes, &numPrimes))
    return;

  printf("pi(%d) = %d\n", N, numPrimes);

  // Load symbols to permute = {'5', '6', '7', '8', '9'}.
  for (i = 0; i < 5; i++)
    letters[i] = '5' + i;

  // Load permutation strings, exiting in case of failure.
  if (!CalculatePermutations(letters, 5, &pPerms))
    {
    free(pPrimes);
    return;
    }

  // Check candidates for primality, last 5 digits a permutation of 56789.
  for (i = 0, count = 1; i < 120; i++)
    {
    x = 10123400000 + atoi(pPerms + 6*i);
    if (IsPrime(x, pPrimes, N))
      printf("%2d : %I64d is prime\n", count++, x);
    }

  // Free permutation and prime memory.
  free(pPerms);
  free(pPrimes);
}


BOOL CalculatePermutations(char letters[], int n, char **ppPerms)
/*
   USE:  Calculate all permutations on n symbols and load into memory
   IN:   letters[] = symbols to permute
         n = number of symbols to permute
   OUT:  *ppPerms = ptr to the permutations (NULL if malloc failure)
   NOTE: Load the permutations one after the other in the order
         generated, with binary 0 after each one : ABC\0BCA\0...
 */
{
  int i;
  int factorial;     // for n! (n factorial)
  char *pPerms;      // ptr to permutations to return to caller
  int permNum;       // to count permutations as generated
  int s[20];         // to hold permutation on {0, 1, 2, ... n-1}
  int k;             // shift position in permutation

  // Returned ptr will be NULL, unless all successful here.
  *ppPerms = NULL;

  // Calc n! = 1 x 2 x....x (n-1) x n = #permutations on n symbols.
  for (i = factorial = 1; i <= n; i++)
    factorial *= i;

  // Malloc for permutation strings, each of size n+1 (+1 for final 0).  
  if ((pPerms = (char *) malloc(factorial * (n + 1))) == NULL)
    return FALSE;

  // Init permutation counter.
  permNum = 0;

  // Init s[] with identity permutation.
  for (i = 0; i < n; i++)
    s[i] = i;

  // Load initial permutation and shift.
  LoadCurPerm(s, n, pPerms, letters);
  ShiftPerm(s, k = n-1);

  while (TRUE)
    {
    // Roll down to where kth element not fixed and shift.
    while ((s[k] == k) && (k > 1))
      ShiftPerm(s, --k);

    // If have rolled down to k == 1, done; else shift k times.
    if ((k == 1) && (s[k] == k))
      break;  // from while
    else
      while (s[k] != k)
        {
        LoadCurPerm(s, n, pPerms + ++permNum * (n + 1), letters);
        ShiftPerm(s, k = n-1);
        }
    }

  // Set returned permutation ptr and return ok.
  *ppPerms = pPerms;
  return TRUE;
}


void ShiftPerm(int x[], int j)
/*
   USE:  Shift permutation's left elements from 0th to jth.
   IN:   x[] = permutation to shift
         j = cut-off point for shifting.
*/
{
  int save, i;

  // Save 0th item.
  save = x[0];

  // Shift 1st thru jth items to left.
  for (i = 0; i < j; i++)
    x[i] = x[i+1];

  // Set jth to 0th, completing the circuit.
  x[j] = save;
}


void LoadCurPerm(int x[], int n, char *pCurPerm, char letters[])
/*
   USE:  Load current permutation into memory.
   IN:   x[] = current permutation on {0, 1, 2, ... n-1}
         n = number of elements in permutation
         pCurPerm = memory address to copy permutation into
         letters[] = symbols for final permutation (mapped to x[])
*/
{
  int i;

  // Copy permutation to memory. Note x[] is permutation on the
  // first n integers (0 -- n-1), letters[x[]] is corresponding
  // permutation on letters[].
  for (i = 0; i < n; i++)
    *(pCurPerm + i) = letters[x[i]];

  // Put 0 at end for string handling.
  *(pCurPerm + n) = 0;
}


BOOL CalculatePrimes(int n, int **ppPrimeList, int *numPrimes)
/*
   USE:  Calculate all the primes up to (and including) n.
   IN:   n = upper limit for primes to produce
   OUT:  *ppPrimeList = ptr to the primes (NULL if malloc failure)
         *numPrimes = pi(n) = number of primes <= n.
   RET:  TRUE if table successfully built, FALSE if not.
   NOTE: Implement the Sieve of Eratosthenes by 'crossing off'
         all the multiples of 2, then all the multiples of 3,
         then 5, and so on, using multiples of all numbers not
         previously crossed off (you needn't cross off multiples
         of 4, for example, because it and all of its multiples
         are already multiples of 2). Use an array to represent
         all the integers <= n. Initialize all elements to 1 and 
         'cross off' by placing a 0 in the slot. At the end of
         the process, all remaining 1s are primes. You need only
         cross off multiples up to the square root of n, because
         a composite (non-prime) number must have a factor <=
         its square root.

         If the function returns successfully, then the prime
         list was malloc'ed and filled in, and the caller must
         free it.
 */
{
  int *pList;        // ptr to prime list to return to caller
  char *primeFlags;  // flags to tell if each int a prime (prime = 1)
  int i;
  int sqRoot;        // sqrt(n) = highest prime to sieve multiples of
  int p;             // primes to sieve multiples of
  int q;             // multiples of base prime p
  int prime;         // primes, ints flagged with 1 after sieving

  // Returned ptr will be NULL, unless all successful here.
  *ppPrimeList = NULL;

  // Malloc flags, one per integer up to (and including) n.
  if ((primeFlags = (char *) malloc(n + 1)) == NULL)
    return FALSE;

  // Set 0 and 1 non-prime and ...
  primeFlags[0] = primeFlags[1] = 0;

  // ... init the rest to prime.
  for (i = 2; i < n; i++)
    primeFlags[i] = 1;

  // sqrt(n) is the upper limit for sieving.
  sqRoot = (int) sqrt(n);

  // p will be the base prime, q its multiples.
  p = q = 2;

  // Outer loop on primes up to sqrt(n).
  while (p <= sqRoot)
    {
    // Sieve multiples of p.
    while ((q+=p) <= n)
      primeFlags[q] = 0;

    // Move p to next prime.
    while (!primeFlags[++p])
      ;

    // Reset q to current prime, p.
    q = p;
    } // while (p < sqRoot)

  // Count the number of primes (1's left in primeFlags[]).
  *numPrimes = 0;
  for (i = 0; i <= n; i++)
    if (primeFlags[i])
      (*numPrimes)++;

  // Malloc for the primes themselves, +1 since index starts at 1.
  if ((pList = (int *) malloc((*numPrimes + 1) * sizeof(int))) == NULL)
    {
    free(primeFlags);
    return FALSE;
    }

  // Copy primes to array (ints flagged with 1).
  for (i = 0, prime = 1; i <= n; i++)
    if (primeFlags[i])
      pList[prime++] = i;
  
  // Free temporary flags array, set returned list ptr, return ok.
  free(primeFlags);
  *ppPrimeList = pList;
  return TRUE;

} // CalculatePrimes()


BOOL IsPrime(__int64 m, int primes[], int numPrimes)
/*
   USE:  Say whether m is prime.
   IN:   m = integer to test for primality.
         primes = table of already existing primes
         numPrimes = number of primes in the table
   RET:  TRUE iff m is prime
   NOTE: Will get a (possibly) false negative if there aren't
         enough primes in the table to make the test.
*/
{
  int sqRoot;        // sqrt(m) = highest prime to divide by
  int i;

  // sqrt(m) is the largest prime needed to divide into m
  sqRoot = (int) sqrt((double)m);

  // Not enough primes on hand.
  if (sqRoot > numPrimes)
    return FALSE;

  // If some prime divides m, then it is composite.
  for (i = 1; primes[i] <= sqRoot; i++)
    if ((m % primes[i]) == 0)
      return FALSE;

  // If no prime divides m, then it is prime.
  return (TRUE);
}