//-----------------------------------------------------------------------------
// Name: prime_count.c
// Uses: Calculate the number of prime number between some range of numbers.
// Date: 02/02/2010
// Author: Andrew Que (http://www.DrQue.net/)
// Revisions:
//  1.0 - 02/02/2010- QUE - Creation.
//
// Build instructions:
//   This software was designed to compile and run in with strict C99 compliance
//   with POSIX threads (pthreads).
//
//   gcc -Wall -std=c99 -pedantic -O3 -pthread -lrt prime_count.c -o prime_count
//
//   Verified on:
//     Ubuntu desktop 9.10 (x86_64), gcc 4.4.1
//     Ubuntu server 9.04 (i686), gcc 4.3.3
//     SunOS Blue-Solaris 5.11; gcc 3.4.3
//     Cygwin 1.5.25; gcc 3.4.4
//
//   Verifion failed on:
//     NetBSD 4.0.1, gcc 4.1.2--crashes shortly after printing banner.
//
//                     (C) Copyright 2010 by Andrew Que
//                        Released as public domain.
//-----------------------------------------------------------------------------
#include <stdio.h>
#include <stdint.h>
#include <stdbool.h>
#include <time.h>
#include <pthread.h>
#include <semaphore.h>

// There happen to be 6542 primes between 2 and 65536.
enum { NUMBER_OF_LOOKUP_PRIMES = 6542 };
static unsigned PrimeNumbers[ NUMBER_OF_LOOKUP_PRIMES ];

// How many numbers to check in a thread.
enum { NUMBERS_PER_THREAD = 0x100000 };

// Number of threads used for calculations.
// (Set this to the number of CPU cores available on the target machine).
enum { NUMBER_OF_THREADS = 4 };

// What range of numbers to check.
static uint32_t const START_NUMBER = 0x2ul;
static uint32_t const END_NUMBER   = 0xFFFFFFFFul;

// Semaphore used to dispatch work to threads.
static sem_t Semaphore;

// Structure to hold worker thread information.
typedef struct
{
  uint32_t StartNumber;   // Where to start.
  uint32_t NumberToCheck; // Numbers to check--usually NUMBERS_PER_THREAD.
  uint32_t NumberFound;   // How many primes were found (return value).
} WORK_TYPE;

//-----------------------------------------------------------------------------
// FUNCTION INFORMATION:
//   Build a lookup table (PrimeNumbers) of all prime numbers between 2 and
// 65536 (or the square root of 2^32).  This function doesn't take long
// despite having to check 2^16-2 values.
//-----------------------------------------------------------------------------
static void GeneratePrimeNumberLookup()
{
  unsigned Index;
  unsigned PrimeNumberCount = 0;

  // First prime number in our list is 2.  We can get all the rest
  // knowing this.
  PrimeNumbers[ PrimeNumberCount++ ] = 2;

  // Get all remaining prime numbers between 3 and 2^16.
  for ( Index = 3; Index < 0x10000ul; ++Index )
  {
    // Assume the number is prime until it is determined to be otherwise.
    bool IsPrime = true;

    // First check: is the number even?
    if ( 0 == ( Index & 1 ) )
    {
      // No even numbers (except 2) are prime.
      IsPrime = false;
    }
    else
    {
      // Check to see if this number is prime...
      unsigned SubIndex = 0;
      while ( ( SubIndex < PrimeNumberCount )
           && ( IsPrime ) )
      {
        // Does it divide evenly by this prime number?
        if ( 0 == ( Index % PrimeNumbers[ SubIndex ] ) )
        {
          // If so, this number isn't prime.
          IsPrime = false;
          break; // <- We can stop checking.
        }

        ++SubIndex;
      }
    }

    // If this number doesn't divide evenly, it is prime...
    if ( IsPrime )
      // Add numbe to our prime list.
      PrimeNumbers[ PrimeNumberCount++ ] = Index;
  }

} // GeneratePrimeNumberLookup

//-----------------------------------------------------------------------------
// FUNCTION USAGE:
//   Return the integer square root of some unsigned 32-bit value.  This
// function uses a power-of-two bit trick such that the function will always
// have a value in 16 iterations.
//
// INPUT:
//    Argument - The number of which to find the square root.
//
// OUTPUT:
//    Integer portion of the square root of "Argument".
//-----------------------------------------------------------------------------
static inline uint16_t SquareRoot( uint32_t Argument )
{
  uint32_t Test;
  uint16_t Root    = 0;
  uint16_t BitMask = ( 1U << 15 );

  // 16 laps.
  while ( BitMask )
  {
    Test = Root + BitMask;

    // Argument >= Test^2?
    if ( Argument >= ( Test * Test ) )
      Root = Test; // <- Use result.

    BitMask >>= 1;
  }

  return Root;

} // SquareRoot

//-----------------------------------------------------------------------------
// FUNCTION USAGE:
//   Test to see if a number is a prime number.  Works on a 32-bit unsigned
// value by dividing it by all prime number up to the square root of the
// number.  This works because because of the nature of prime numbers.  A
// number (call it x) is prime if there are no two number (call them a and b)
// such that a * b = x.  All non-prime numbers can be expressed as the sum of
// two or more prime numbers.  For example, 125 can be made from 5 * 25, but
// 25 can be made from 5 * 5.  So 5 * 5 * 5 = 125, and represents the most
// factored version of 125.  This is true of any number.  Since it takes at
// least two prime numbers to create a factor, we only need to check the primes
// up to x^1/2 (or the square root) of the number.  This is because the if
// x = a * a, then x^1/2 = a.  If x = a * b, and b is greater a, then a must
// be less then x^1/2. Thus, we only need to check primes up to x^1/2.
//   Since the input is a 32-bit number, maximum number that can be represented
// is 2^32-1.  (2^32)^1/2 = 2^16.  So we need to check all primes up to 2^16.
// To do this, we keep a lookup table of all primes in this range.
//
// INPUT:
//   Number - A unsigned 32-bit value to test.
//
// OUTPUT:
//   Returns true if Number is prime, false if not.
//-----------------------------------------------------------------------------
static inline bool IsPrime( uint32_t Number )
{
  // Have we not yet build the prime number table?
  if ( 0 == PrimeNumbers[ 0 ] )
    GeneratePrimeNumberLookup();

  // Assume the number is prime until it is determined to be otherwise.
  bool IsPrime = true;

  // Is number even (and not the number 2)?
  if ( ( 0 == ( Number & 1 ) )
    && ( 2 != Number ) )
  {
    // No even numbers (except 2) are prime.
    IsPrime = false;
  }
  else
  {
    // We only need to check up to the square root of the number.
    uint16_t Root = SquareRoot( Number );
    unsigned Index = 1; // <- Start with 3.

    // While we still have prime numbers to test, the number is less
    // then the squre root of the number, and nothing so far has divided
    // evenly...
    while ( ( Index < NUMBER_OF_LOOKUP_PRIMES )
         && ( PrimeNumbers[ Index ] <= Root )
         && ( IsPrime ) )
    {
      // Does this prime divide into the number?
      if ( 0 == ( Number % PrimeNumbers[ Index ] ) )
        // Then the number is not prime.
        IsPrime = false;

      ++Index;
    }
  }

  // Return the results.
  return IsPrime;

} // IsPrime

//-----------------------------------------------------------------------------
// FUNCTION USAGE:
//   Thread used to count the number of primes in a given range of numbers.
// The range checked is from the 32-bit unsigned integer pointed to by
// ArgumentPointer to ArgumentPointer + NUMBERS_PER_THREAD.
//
// INPUT:
//   ArgumentPointer - A pointer to a 32-bit unsigned integer that contains
// the first number to check.
//
// OUTPUT:
//   The function itself returns nothing.  The unit global "NumberOfPrimes" is
// updated by the number of primes found.
//-----------------------------------------------------------------------------
static void * PrimeThread( void * ArgumentPointer )
{
  // Get the work data passed to the thread.
  WORK_TYPE * Data = (WORK_TYPE *)ArgumentPointer;
  uint32_t Number = Data->StartNumber;
  uint32_t Count = 0;
  unsigned Index;

  // For all the numbers to check...
  for ( Index = 0; Index < Data->NumberToCheck; ++Index )
  {
    // Is this number a prime?
    if ( IsPrime( Number ) )
      // Then count it.
      ++Count;

    // Next number.
    ++Number;
  }

  // Save results.
  Data->NumberFound = Count;

  // This thread is now complete.  Release one count from the dispatch
  // semaphore.
  sem_post( &Semaphore );

  // End this thread.
  pthread_exit( 0 );

  // Never reached--here for language consistency.
  return 0;

} // PrimeThread

//-----------------------------------------------------------------------------
// FUNCTION USAGE:
//   Program main function.  This function will setup the dispatch semaphore,
// and work threads that will count all the prime number in a range given by
// the unit globals START_NUMBER and END_NUMBER.  The total count is displayed when the
// program completes.
//
// OUTPUT:
//   This function (and program as a whole) always returns 0.
//-----------------------------------------------------------------------------
int main()
{
  // Print a header.
  printf( "============================\n" );
  printf( "Prime number count\n" );
  printf( "============================\n" );
  printf
  ( 
    "Counting the number of primes between %u and %u\n", 
    (unsigned)START_NUMBER, (unsigned)END_NUMBER 
  );

  // Mark the time this program began.
  time_t StartTime = time( NULL );

  // Worker threads.
  pthread_t Threads[ NUMBER_OF_THREADS ];

  // Data storage for threads.
  WORK_TYPE Data[ NUMBER_OF_THREADS ];

  // Setup this dispatch semaphore such that it can handle NUMBER_OF_THREADS counts
  // before it blocks the request.
  sem_init( &Semaphore, NUMBER_OF_THREADS, NUMBER_OF_THREADS );

  // Index for the next available thread.
  unsigned ThreadIndex;

  // Zero out thread data (used so we can tell if the thread has been used
  // yet or not).
  for ( ThreadIndex = 0; ThreadIndex < NUMBER_OF_THREADS; ++ThreadIndex )
    Data[ ThreadIndex ].NumberToCheck = 0;

  // Starting number to pass to the next work thread.
  uint32_t Number = START_NUMBER;
  uint32_t NumbersLeft = END_NUMBER - START_NUMBER;

  // Total number of primes found so far.
  unsigned NumberOfPrimes = 0;

  // Zero thread index.
  ThreadIndex = 0;

  // Loop until all the numbers have been checked...
  while ( NumbersLeft )
  {
    // Display progress.
    printf( "%08X => %u\r", (unsigned)Number, (unsigned)NumberOfPrimes );
    fflush( stdout ); // <- Make sure the screen is updated.

    // Wait for a free worker thread.
    sem_wait( &Semaphore );

    // Was this thread running?
    if ( Data[ ThreadIndex ].NumberToCheck )
    {
      // Rejoin the thread--should be finished now.
      pthread_join( Threads[ ThreadIndex ], NULL );

      // Accumulate the number of primes found in this thread.
      NumberOfPrimes += Data[ ThreadIndex ].NumberFound;
    }

    // How many number to check.
    if ( NumbersLeft > NUMBERS_PER_THREAD )
      Data[ ThreadIndex ].NumberToCheck = NUMBERS_PER_THREAD;
    else
      Data[ ThreadIndex ].NumberToCheck = NumbersLeft;

    // Create a worker thread to check this number set.
    Data[ ThreadIndex ].StartNumber = Number;
    pthread_create
    (
      &Threads[ ThreadIndex ],
      NULL,
      PrimeThread,
      (void *)&Data[ ThreadIndex ]
    );

    // Move to next number set.
    Number += Data[ ThreadIndex ].NumberToCheck;
    NumbersLeft -= Data[ ThreadIndex ].NumberToCheck;

    // Advance thread index with wrap around.
    ++ThreadIndex;
    if ( ThreadIndex >= NUMBER_OF_THREADS )
      ThreadIndex = 0;
  }

  // At this point, all work has been dispatched.  We just need to wait for
  // the worker threads to finish.

  // Denote the current state.
  printf( "Finishing...              \r" );
  fflush( stdout ); // <- Make sure the screen is updated.

  // Wait for each thread to finish.
  for ( ThreadIndex = 0; ThreadIndex < NUMBER_OF_THREADS; ++ThreadIndex )
  {
    // Was this thread running?
    if ( Data[ ThreadIndex ].NumberToCheck )
    {
      // Wait for thread to finish.
      pthread_join( Threads[ ThreadIndex ], NULL );

      // Accumulate the number of primes found in this thread.
      NumberOfPrimes += Data[ ThreadIndex ].NumberFound;
    }
  }

  // Let go of dispatch semaphore.
  sem_destroy( &Semaphore );

  // Calculate how long the program ran.
  unsigned ElapsedTime = (unsigned)difftime( time( NULL ), StartTime );

  // Display results.
  printf
  (
    "Found %u primes between %u and %u, %u seconds.\n",
    (unsigned)NumberOfPrimes,
    (unsigned)START_NUMBER,
    (unsigned)END_NUMBER,
    ElapsedTime
  );

  // Done, exit with no error.
  return 0;

} // main

//--------------------------------------=--------------------------------------