Andrew Que Sites list Photos
Projects Contact
Main

July 05, 2021

C11 Prime Number Counter

Back in February of 2010 I wrote a simple program to count all the prime numbers between 2 and 232.  It use pthreads so I could fully utilize the duel processor Red-Dragon.  The C11 standard introduced threads right into the language.  Several years ago I thought I'd port my prime counting program to use this because any system that supported the full C11 standard would be able to compile and run it.  Sadly I discovered that gcc had not implemented C11 threads at that time.  However, as of version GCC 9.3.0 C11 threads are present—you just have to link with the pthreads library as it uses pthreads under the hood. 

Porting didn't take too long.  However, there was one item I did not have from the start that I would need: counting semaphores.  Counting semaphores are used by the program to keep some specified number of running threads—typically one thread for each core of the CPU.  The C11 threads implementation does have mutexes and c conditional waits and those can be used to make a counting semaphore.  So I wrote a very simple header file with inline functions for counting semaphores.

//-----------------------------------------------------------------------------
// Name: primeCount.c
// Uses: Calculate the number of prime number between some range of numbers.
// Date: 2010-02-02
// Author: Andrew Que (https://www.DrQue.net/)
// Revisions:
//  1.0 - 2010-02-02 - Creation.
//  1.1 - 2014-08-22 - Bug fix.  Prime table needs to be initialized
//    outside of threads or problems could occur.
//  1.2 - 2021-07-05 - Converted to C11 threads.
//
// Build instructions:
//   This software was designed to compile and run in with strict C11
//   compliance with threads (pthreads).
//
//   Linux:
//     gcc -Wall -Wextra -Werror -pedantic -std=c11 -O2 primeCount.c -o primeCount -lpthread
//
//                 (C) Copyright 2010,2014,2021 by Andrew Que
//                         Released as public domain.
//-----------------------------------------------------------------------------
#include <stdio.h>
#include <stdint.h>
#include <stdbool.h>
#include <time.h>
#include <threads.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 = 16 };

// What range of numbers to check.
static uint32_t const START_NUMBER = UINT32_C( 0x2 );
static uint32_t const END_NUMBER   = UINT32_C( 0xFFFFFFFF );

// Semaphore used to dispatch work to threads.
static Semaphore 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).
} WorkType;

//-----------------------------------------------------------------------------
// Uses:
//   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.
  // Count by 2 since all even numbers are divisible by two.
  for ( index = 3; index < UINT32_C( 0x10000 ); index += 2 )
  {
    // Assume the number is prime until it is determined to be otherwise.
    bool isPrime = true;

    // 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;
      }

      ++subIndex;
    }

    // If this number doesn't divide evenly, it is prime...
    if ( isPrime )
      // Add numbe to our prime list.
      primeNumbers[ primeNumberCount++ ] = index;
  }

} // generatePrimeNumberLookup

//-----------------------------------------------------------------------------
// Uses:
//   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

//-----------------------------------------------------------------------------
// Uses:
//   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 )
{
  // 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

//-----------------------------------------------------------------------------
// Uses:
//   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 int primeThread( void * argumentPointer )
{
  // Get the work data passed to the thread.
  WorkType * data = (WorkType *)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.
  semaphoreRelease( &semaphore );

  // End this thread.
  thrd_exit( 0 );

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

} // primeThread

//-----------------------------------------------------------------------------
// Uses:
//   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( "============================ " );
  printf( "Prime number count " );
  printf( "============================ " );
  printf
  (
    "Counting the number of primes between %u and %u ",
    (unsigned)START_NUMBER, (unsigned)END_NUMBER
  );

  generatePrimeNumberLookup();

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

  // Worker threads.
  thrd_t threads[ NUMBER_OF_THREADS ];

  // data storage for threads.
  WorkType data[ NUMBER_OF_THREADS ];

  // Setup this dispatch semaphore such that it can handle NUMBER_OF_THREADS counts
  // before it blocks the request.
  semaphoreInit( &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 "(unsigned)number, (unsigned)numberOfPrimes );
    fflush( stdout )// <- Make sure the screen is updated.

    // Wait for a free worker thread.
    semaphoreWait( &semaphore );

    // Was this thread running?
    if ( data[ threadIndex ].numberToCheck )
    {
      // Rejoin the thread--should be finished now.
      thrd_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;
    thrd_create
    (
      &threads[ threadIndex ],
      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...               " );
  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.
      thrd_join( threads[ threadIndex ], NULL );

      // Accumulate the number of primes found in this thread.
      numberOfPrimes += data[ threadIndex ].numberFound;
    }
  }

  // Let go of dispatch semaphore.
  semaphoreDestroy( &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. ",
    (unsigned)numberOfPrimes,
    (unsigned)START_NUMBER,
    (unsigned)END_NUMBER,
    elapsedTime
  );

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

} // main

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

The implementation is nearly identical, just using the C11 thread structures and my semaphore unit.  I don't know the speeds of the original Red Dragon, but I ran a speed test in May of 2017.  My fastest machine at the time, a duel-core 4-thread Intel i7, required 28.26 minutes.  My AMD Ryzen 7 1700, 8-core/16-thread CPU needs 16.58 minutes.

Prime Counter v1.2

SHA-256: cfa16f3ca5c38d0c99f8181eeeac440c84ac4a942e5b07ad5d2afe39b4e67692

   Noah and Mary came to visit today and we decided to have Indian food at my favorite place across the street.  Turns out that of all the people I have taken out for Indian, Mary was not one of them.  She had never had Indian food before.  She ordered my recommendation and quickly fell in love.  Indian food will do that.

July 03, 2021

Eve Liberty Goes to Her New Home

Aislinn and Eve Liberty

Aislinn and Eve Liberty

   Today Eve Liberty leaves me to begin her new life with Aislinn as Aislinn's first car.  Eve has over 300,000 miles, but refuses to count past 299,999.  I understood that she would be unable to take any more +4,000 mile road trips with me.  However, Eve does just fine around town and still has plenty of life left in her.  So I decided Eve needed to go to a good home. 
   I met Aislinn shortly after she was born and have known her mom for over 20 years.  I taught Aislinn how to roller skate and she skated in the Beloit rink until they closed.  Knowing Aislinn didn't yet have a car, I thought she would be able appreciate Eve.  So there you go Aislinn.  And I hope you end up loving Eve as much as I did.