Arduino generating prime numbers

I got inspired after watching the guys on Numberphile unbox a gadget that displays prime numbers:

I made a quick version that did a similar thing with an 8-digit LED display:

Totally useless of course, but strangely mesmerizing to watch.

Hardware: Arduino Uno, plus one of these 8-digit display boards from eBay for about $10.

References:

Code:

// Prime number generator
//
// Author: Nick Gammon
// Date:   8th October 2013

#include <SPI.h>

const int SHOW_EVERY = 500;  // how often to echo a prime to the serial port
const unsigned long TIME_TO_WAIT = 500;  // mS

// MAX7219 registers

const byte MAX7219_REG_NOOP        = 0x0;
// codes 1 to 8 are digit positions 1 to 8
const byte MAX7219_REG_DECODEMODE  = 0x9;
const byte MAX7219_REG_INTENSITY   = 0xA;
const byte MAX7219_REG_SCANLIMIT   = 0xB;
const byte MAX7219_REG_SHUTDOWN    = 0xC;
const byte MAX7219_REG_DISPLAYTEST = 0xF;

// 7-segments patterns for digits 0 to 9
const byte digits [10] = {
  0b1111110,  // 0
  0b0110000,  // 1
  0b1101101,  // 2
  0b1111001,  // 3
  0b0110011,  // 4
  0b1011011,  // 5
  0b1011111,  // 6
  0b1110000,  // 7
  0b1111111,  // 8
  0b1111011,  // 9
};


// send one byte via  SPI to MAX7219
void sendByte (const byte reg, const byte data)
  {    
  digitalWrite (SS, LOW);
  SPI.transfer (reg);
  SPI.transfer (data);
  digitalWrite (SS, HIGH); 
  }  // end of sendByte

// send one character (data) to position (pos) with or without decimal place
// pos is 0 to 7
// data can be '0' to '9' plus various other things as below
void sendChar (const byte pos, const char data, const bool dp = false)
  {
  byte converted;  
  switch (data)
    {
    case '0' ... '9' : converted = digits [data - '0']; break;
    case '-':          converted = 0b0000001;  break;  
    case 'A':          converted = 0b1110111;  break;  
    case 'b':          converted = 0b0011111;  break;  
    case 'c':          converted = 0b0001101;  break;  
    case 'C':          converted = 0b1001111;  break;  
    case 'd':          converted = 0b0111101;  break;  
    case 'E':          converted = 0b1001111;  break;  
    case 'F':          converted = 0b1000111;  break;  
    case 'h':          converted = 0b0010111;  break;  
    case 'H':          converted = 0b0110111;  break;  
    case 'L':          converted = 0b0001110;  break;  
    case 'o':          converted = 0b0011101;  break;  
    case 'P':          converted = 0b1100111;  break;  
    case 'r':          converted = 0b0000101;  break;  
    case 'S':          converted = 0b1011011;  break;  
    case 't':          converted = 0b0000111;  break;  
    case 'u':          converted = 0b0011100;  break;  
    case ' ':          converted = 0b0000000;  break;  
    case 'I':          converted = digits [1]; break;
    
    default:           converted = 0b0000001;  break;  // -
    } // end of switch  
  if (dp)
    converted |= 0b10000000;
  sendByte (8 - pos, converted);
  }  // end of sendChar
  
// write an entire null-terminated string to the LEDs
void sendString (const char * s)
{
  byte pos;
  
  for (pos = 0; pos < 8 && *s; pos++)
    {
    boolean dp = s [1] == '.';
    sendChar (pos, *s++, dp);   // turn decimal place on if next char is a dot
    if (dp)  // skip dot
      s++;
    }
    
  // space out rest
  while (pos < 8)
    sendChar (pos++, ' ');
}  // end of sendString


long candidate;
long found = 5; // Number we found
int count = found - 1;


void setup() {
  Serial.begin(115200); 
  
  while (!Serial) { }
  
  Serial.println ();
  Serial.println ();
  
  SPI.begin ();

  sendByte (MAX7219_REG_SCANLIMIT, 7);      // show 8 digits
  sendByte (MAX7219_REG_DECODEMODE, 0);     // use bit patterns
  sendByte (MAX7219_REG_DISPLAYTEST, 0);    // no display test
  sendByte (MAX7219_REG_INTENSITY, 7);      // character intensity: range: 0 to 15
  sendString ("");                          // clear display
  sendByte (MAX7219_REG_SHUTDOWN, 1);       // not in shutdown mode (ie. start it up)

}

unsigned long start;
unsigned long lastShown;

void loop() 
  {

  for (candidate = 3; candidate < 99999999; candidate += 2)
    showIfPrime (candidate);
 
  }  // end of loop

void showIfPrime (long num) 
  {
  
  // we are already skipping odd numbers, now skip if divisible by 3
  if (num % 3 == 0) 
    return; 

  // Only have to check for divisible for the sqrt(number)
  long upper = sqrt(num) + 6;
  
  // Check if the number is divisible (start at 6 going up)
  for (long cnum = 6; cnum <= upper; cnum += 6)
    {
    if (num % (cnum - 1) == 0) 
      return; 
    if (num % (cnum + 1) == 0) 
      return; 
    }  // end of checking cnum-1, cnum+1

 // echo to serial port for validating
 if (++count >= SHOW_EVERY)
  {
  Serial.print(candidate);
  Serial.print(" is prime. Took ");
  Serial.print(millis () - start);
  Serial.print (" mS. Found ");
  Serial.print (found);
  Serial.println (" primes.");
  start = millis ();
  count = 0;
  }

  found++;

  // delay until interval is up
  // (this absorbs the calculation time)
  while (millis () - lastShown < SHOW_EVERY)
    { }
    
  char buf [10];
  sprintf (buf, "%8ld", candidate);
  sendString (buf);
  lastShown = millis ();
  
}  // end of showIfPrime

The code was written in a bit of a hurry. It doesn't show the first few primes (well, we know what they are). I think it starts at 7, because of the nature of the algorithm. When I timed it at around the 128201 mark (the prime number 128201) it was taking about 10 mS per prime. Of course, this makes the display tick over unreadably fast, so it is throttled back to show a number every 500 mS.

The algorithm skips even numbers (naturally) and also numbers divisible by 3. Then it tests in jumps of 6, plus or minus 1, to see if they are primes, up to the square root of the number we are testing for. It isn't particularly fast, but also doesn't use a lot of memory.

The next step is to start launching them into space towards candidate planets taken from the Kepler telescope's web site.

Hi Nick,
Nice prime machine, as every other programmer did some prime-finding too, see the 30k +i optimization here

The nice part of this trick (of the 30K method) is that you can encode the set of eight possible primes into one byte. imagine an array in which the arrayindex indicates the numberx30 and every bit indicates one of the eight mentioned offsets. a 1 bit means prime, a 0 bit means no prime.

Just as an update, it's been running all night (I'm guessing around 15 hours) and it is up to 1442723. So still the last digit to go.

Once again though, it is throttled back to update every 500 mS, and just by eye it looks like it is able to compute the next number in that time.

If you want to try it at home (without the LED display) there is also serial output. For debugging it prints every 500th prime (you can change that near the start). And you can change the delay between generating primes.

Program size:

Binary sketch size: 6,286 bytes (of a 32,256 byte maximum)

RAM usage unknown, presumably not much.

I was right, it seems. Using WolframAlpha to calculate:

PrimePi(1442723) = 110117

Therefore 1,442,723 is the 110,117th prime number. So at two primes a second it should have taken:

110117 / 2 / 60 / 60 = 15.29 hours

No such thing! I will admit to a penchant for "useless" projects myself but I'd also argue that it's the rare project that is really useless. Usually there is something to be learned or some fun to be had.

At work we used to say that no project is useless because it can always be used as a bad example.

But I like this one! Can't wait to try it later.

I thought you might want to know ...

That's at 5:20 pm the day after it started. So, around 24 hours later.

As you can see, it will be around 5 days before it ticks over to the 8th digit.

This sort of thing might look good in a shop window. :wink:

About how long between numbers at this point?

Next morning now. I am up to 4020503 and the time between numbers seems to be 500 mS still.

By my reckoning the number of tests it has to do in the worst case is:
(?4020503) / 4 = 501

Then in most cases it will quickly find the number divisible by 3, 5, 7 and so on. So probably testing all 501 would be something fairly infrequent.

To test that timing, I wrote this:

volatile long foo = 4020503;
volatile long bar;

void setup ()
  {
  Serial.begin (115200);
  Serial.println ();
  unsigned long start = micros ();
  for (long i = 0; i < 501; i++)
    {
    bar = foo % i;  
    }
  Serial.println (micros () - start);
  }  // end of setup

void loop () { }

Output:

20924

So, 21 mS to do 501 modulus operations. It's no wonder it can keep up, since I am giving it 500 mS.

Of course, that would only apply for adjacent primes. Still, the gap between primes doesn't seem to often be much greater than 50 or so.

The average gap between primes is apparently ln(p) where p is the prime we are looking at. So:

ln (4020503) = 15.2

Thus I would expect an average gap of 15. And of course, not all of those are being tested (like even numbers).

Sweet. I need to get one of those 8-digit displays. I'm sure I have the parts to put one together but it'd be nice to have it on a compact board like that.

I'm up to 8 digits!

Still going strong, it's still doing one every 500 mS.

Those display modules aren't chainable, are they?

How is your chip at dividing long longs? Too bad there isn't a 48-bit integer datatype.

Absolutely they are chainable. I had a test set up yesterday. Code here:

I'm not sure about the speed of doing long longs. I could find out, I suppose. :slight_smile:

The code below generates a 16-digit prime number using two daisy-chained displays, using long long arithmetic.

// Prime number generator using long long
//
// Author: Nick Gammon
// Date:   22nd October 2013

#include <SPI.h>
#include <EEPROM.h>
#include <EEPROMAnything.h>
#include <limits.h>

const int CHIP_COUNT = 2;  // how many MAX7219s

const int SHOW_EVERY = 100;  // how often to echo a prime to the serial port (and update EEPROM)
const unsigned long TIME_TO_WAIT = 500;  // mS
const long long FIRST_PRIME_TO_USE = 99999999999LL;  // if ULONG_MAX then read from EEPROM

// MAX7219 registers

const byte MAX7219_REG_NOOP        = 0x0;
// codes 1 to 8 are digit positions 1 to 8
const byte MAX7219_REG_DECODEMODE  = 0x9;
const byte MAX7219_REG_INTENSITY   = 0xA;
const byte MAX7219_REG_SCANLIMIT   = 0xB;
const byte MAX7219_REG_SHUTDOWN    = 0xC;
const byte MAX7219_REG_DISPLAYTEST = 0xF;

// 7-segments patterns for digits 0 to 9
const byte digits [10] = {
  0b1111110,  // 0
  0b0110000,  // 1
  0b1101101,  // 2
  0b1111001,  // 3
  0b0110011,  // 4
  0b1011011,  // 5
  0b1011111,  // 6
  0b1110000,  // 7
  0b1111111,  // 8
  0b1111011,  // 9
};


void sendToAll (const byte reg, const byte data)
  {    
  digitalWrite (SS, LOW);
  for (int chip = 0; chip < CHIP_COUNT; chip++)
    {
    SPI.transfer (reg);
    SPI.transfer (data);
    }
  digitalWrite (SS, HIGH); 
  }  // end of sendToAll


void sendChar (const int which, const char c)
  {
  byte i;
  
  // segment is in range 1 to 8
  const byte segment = 8 - (which % 8);
  // for each daisy-chained display we need an extra NOP
  const byte nopCount = which / 8;
  // start sending
  digitalWrite (SS, LOW);
  // send extra NOPs to push the data out to extra displays
  for (byte i = 0; i < nopCount; i++)
    {
    SPI.transfer (MAX7219_REG_NOOP);
    SPI.transfer (MAX7219_REG_NOOP);  // need 16 bits of NOP
    }    
  // send the segment number and data
  SPI.transfer (segment);
  SPI.transfer (c);
  // start with enough NOPs so later chips don't update
  for (int i = 0; i < CHIP_COUNT - nopCount - 1; i++)
    {
    SPI.transfer (MAX7219_REG_NOOP);
    SPI.transfer (MAX7219_REG_NOOP);  // need 16 bits of NOP
    }    
  // all done!
  digitalWrite (SS, HIGH); 
  }  // end of sendChar
  
// write an entire null-terminated string to the LEDs
void sendString (const char * s)
{
  for (int pos = 0; *s; pos++)
    sendChar (pos, *s++);
}  // end of sendString


long long candidate;
long found = 5; // Number we found
int count = found - 1;


void setup() {
  Serial.begin(115200); 
  
  while (!Serial) { }
  
  Serial.println ();
  Serial.println ();
  
  SPI.begin ();

  sendToAll (MAX7219_REG_SCANLIMIT, 7);      // show 8 digits
  sendToAll (MAX7219_REG_DECODEMODE, 0xFF);  // use digits (not bit patterns)
  sendToAll (MAX7219_REG_DISPLAYTEST, 0);    // no display test
  sendToAll (MAX7219_REG_INTENSITY, 15);      // character intensity: range: 0 to 15
  sendToAll (MAX7219_REG_SHUTDOWN, 1);       // not in shutdown mode (ie. start it up)
  
  candidate = FIRST_PRIME_TO_USE;
  
  if (candidate == ULONG_MAX)
    EEPROM_readAnything (0, candidate);
  else
    EEPROM_writeAnything (0, candidate);

}

unsigned long start;
unsigned long lastShown;

void loop() 
  {

  for ( ; candidate < 9999999999999999LL; candidate += 2)
    showIfPrime (candidate);
 
  candidate = 3;  // back to start!
  }  // end of loop

void showIfPrime (long long num) 
  {
  
  // we are already skipping odd numbers, now skip if divisible by 3
  if (num % 3 == 0) 
    return; 

  // Only have to check for divisible for the sqrt(number)
  long long upper = sqrt(num) + 6;
  
  // Check if the number is divisible (start at 6 going up)
  for (long long cnum = 6; cnum <= upper; cnum += 6)
    {
    if (num % (cnum - 1) == 0) 
      return; 
    if (num % (cnum + 1) == 0) 
      return; 
    }  // end of checking cnum-1, cnum+1

  char buf [20] = { 0 };
  long long temp = candidate;
  byte pos = 0;
  
  // make printable
  for (pos = 0; pos < 16; pos++)
    {
    if (temp == 0)
      break;
    char digit = temp % 10;
    buf [15 - pos] = digit | '0';
    temp /= 10;
    }

// insert leading spaces
  for (; pos < 16; pos++)
    buf [15 - pos] = 15;
    
 // echo to serial port for validating
 if (++count >= SHOW_EVERY)
  {
  Serial.print(buf);
  Serial.print(" is prime. Took ");
  Serial.print(millis () - start);
  Serial.print (" mS. Found ");
  Serial.print (found);
  Serial.println (" primes.");
  start = millis ();
  count = 0;
  
  for (byte i = 0; i < 5; i++)
    {
    sendString ("\x0F\n\x0F\n\x0F\n\x0F\n\x0F\n\x0F\n\x0F\n\x0F\n\x0F");  // warn not to turn off
    delay (500);
    sendString ("\x0F\x0F\x0F\x0F\x0F\x0F\x0F\x0F\x0F\x0F\x0F\x0F\x0F\x0F\x0F\x0F");  // blink
    delay (500);
    }  // end of for
  
  EEPROM_writeAnything (0, candidate);  // for restarting next time
  }  // end of if time to write to serial and update EEPROM

  found++;

  // delay until interval is up
  // (this absorbs the calculation time)
  while (millis () - lastShown < SHOW_EVERY)
    { }
  
  sendString (buf);
  lastShown = millis ();
  
}  // end of showIfPrime

This has slowed it down somewhat, as you might expect. The time between calculating 100000000063 and 100000000069 was 28 seconds.

Nick, is the 8-digit version still running? Wondering if it's rolled over yet. I finally breadboarded one up here. I took the delay out, it's been running 10 or 12 hours, it's up to 13,000,000 and change. It will also save the time to EEPROM when it rolls over. The longer it runs, the slower it goes...

The 8-digit one is up to 444xxxxx. I'm saving the current number to EEPROM every 500 numbers (around every 4 minutes). That's just in case the power goes. It still seems to be keeping up with the rate of two every second.

How is the compiler with division? From what I heard, it leaves something to be desired.

You might want to optimize. See:
http://www.steike.com/code/bits/division/

What I am suggesting is this:
You don't need to do modulus the normal way. Instead, since you only want to test for divisibility, you can try a trick: assume that the division is exact, try to divide from right to left, and see whether you get a valid answer or just garbage. For this, it would help to have a table of multiplicative inverses modulo 256.