Go Down

### Topic: Arduino generating prime numbers (Read 33199 times)previous topic - next topic

#### nickgammon

#15
##### Oct 21, 2013, 09:35 pm
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.
Please post technical questions on the forum, not by personal message. Thanks!

#### nickgammon

#16
##### Oct 21, 2013, 10:16 pm
The code below generates a 16-digit prime number using two daisy-chained displays, using long long arithmetic.

Code: [Select]
`// 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 MAX7219sconst int SHOW_EVERY = 100;  // how often to echo a prime to the serial port (and update EEPROM)const unsigned long TIME_TO_WAIT = 500;  // mSconst long long FIRST_PRIME_TO_USE = 99999999999LL;  // if ULONG_MAX then read from EEPROM// MAX7219 registersconst byte MAX7219_REG_NOOP        = 0x0;// codes 1 to 8 are digit positions 1 to 8const 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 9const 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 sendToAllvoid 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 LEDsvoid sendString (const char * s){  for (int pos = 0; *s; pos++)    sendChar (pos, *s++);}  // end of sendStringlong long candidate;long found = 5; // Number we foundint 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 loopvoid 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.

Please post technical questions on the forum, not by personal message. Thanks!

#### Jack Christensen

#17
##### Oct 23, 2013, 03:54 pm

I'm up to 8 digits!

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...

#### nickgammon

#18
##### Oct 23, 2013, 09:38 pm
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.
Please post technical questions on the forum, not by personal message. Thanks!

#### odometer

#19
##### Oct 27, 2013, 11:12 am
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.

#### nickgammon

#20
##### Oct 27, 2013, 09:21 pm
If you want to modify my code further up you are welcome to demonstrate this.

I should point out that this is using long long type (64 bits) which probably makes the processor work fairly hard. I don't have memory for a lot of inverses, and indeed the problem with a prime number like 100000000063 is that you don't get any quick matches, so it has to do a lot of divisions (obviously excluding even numbers) to make sure it is prime.

Code: [Select]
`?100000000063 / 4 = 79057`

So I suppose it is OK that it takes 23 seconds to do 79057 long long divisions. But I am certainly open to faster methods, bearing in mind such methods only have 2 kB of RAM to work in.
Please post technical questions on the forum, not by personal message. Thanks!

#### odometer

#21
##### Oct 27, 2013, 11:41 pm
First of all, I hope your square root is working properly.
If you are going to convert to a float before taking the square root, you lose precision, and possibly your square root will come out too low-- you know what that means.
In which case, you might want to get at your limit using something like:
Code: [Select]
`limit = (long long)sqrt((float)n)limit = max(limit, n/limit)`
What I meant by avoiding division was this:

Begin by assuming that the division does come out even, and go ahead and divide.

For example, let's say you are testing 21005 for divisibility by 5. You and I can at once see it divides, but remember, the Arduino is testing 0x520D for divisibility by 0x05.

We go ahead and divide.

Hex number ends in: 0 1 2 3 4 5 6 7 8 9 A B C D E F
Quintuple ends in: 0 5 A F 4 9 E 3 8 D 2 7 C 1 6 B

Which we can rearrange as:

Quintuple ends in: 0 1 2 3 4 5 6 7 8 9 A B C D E F
Hex number ends in: 0 D A 7 4 1 E B 8 5 2 F C 9 6 3

Note that this, as well, is an arithmetic progression modulo sixteen.
Multiply the "quintuple" by 0xD to get the "original" number.
So we don't need the whole table, just the entry for 1.

0x520D ends in D. 0xD times 0xD gives us 0xA9. So our quotient ends in 9.
0x520D - (0x9 * 0x5) = 0x51E0, and we slough off the final 0.
0x51E ends in E. 0xE times 0xD gives us 0xB6. So our next quotient digit is 6.
0x51E - (0x6 * 0x5) = 0x500, and we slough off the final 0.
0x50 ends in 0. 0x0 times 0xD gives us 0x0. So our next quotient digit is 0.
0x50 - (0x0 * 0x5) = 0x50, and we slough off the final 0.
0x5 ends in 5. 0x5 times 0xD gives us 0x41. So our next quotient digit is 1.
0x5 - (0x1 * 0x5) = 0x0, so we're done.
The quotient is 0x1069.
Check by multiplication: 0x1069 * 0x5 = 0x520D

Because our dividend has only four digits, we know we must terminate within four digits. If the algorithm had not terminated soon enough, we would know that we were getting a garbage result, and that the division did not come out even after all.

#### robtillaart

#22
##### Oct 28, 2013, 02:17 pm
Quote
But I am certainly open to faster methods, bearing in mind such methods only have 2 kB of RAM to work in.

you can check the 30K + i method as I referred to in #2  it should be a bit faster as it only tests 8 numbers of every 30 (is approx 1 in 4) where the 6k+i test 2 numbers in 6  (1 in 3) ~~> gain 10%?
Rob Tillaart

Nederlandse sectie - http://arduino.cc/forum/index.php/board,77.0.html -
(Please do not PM for private consultancy)

#### nickgammon

#23
##### Oct 28, 2013, 08:35 pm
Sure, but the version using longs does one every 1/2 seconds which I am fine with, and for the long long version saving 10% over 27 seconds isn't that impressive.

I was hoping for one of those "probability" primality tests to work, but my brain isn't quite up to the maths right now.
Please post technical questions on the forum, not by personal message. Thanks!

#### robtillaart

#24
##### Oct 29, 2013, 01:59 pm
A substantial gain could be made if you could get from long long (64)  to long (32)
Not tested, but just an idea

first one tests 6k-1   if (n % (k-1) == 0) return
then one test  6k+1  if (n % (k+1) ==0) retrun

One might do the following first:
let p = n /(k+1);
then  n%(k+1) <==>  n - p*(k+1);

This takes an equal amount of time, but with p we have an estimation for  q = n/(k-1):  q >= p

With this estimation we can simplify the #digits
then  n%(k-1)  <==>  (n-p*(k-1)) % (k-1)

the multiplication p*(k-1) is relative fast and  if  (n- p*(k-1)) would fit in a long  the remaining % would be 32 bit iso 64 bit.

Note a long % is 6 times faster than a long long  %

Rob Tillaart

Nederlandse sectie - http://arduino.cc/forum/index.php/board,77.0.html -
(Please do not PM for private consultancy)

#### odometer

#25
##### Oct 29, 2013, 04:22 pm
Do you even need to use the divide or modulo operators?
If the / and % operators were forbidden to you, how then would you go about this?

#### robtillaart

#26
##### Oct 29, 2013, 05:08 pm
a well known method without division is the sieve
1) make a bool  array on 2..n, set them all to true

2)  take next element from array (first elm == 2) and print it
3) remove all multiples of this element from the array
4) goto step 2

array = 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 ...
take next element (2)
remove all multiples
array = 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 ...
take next element (3)
remove all multiples
array = 5 7 11 13 15 17 19 23 25 29 31 35 ...
take next element (5)
remove all multiples
array = 7 11 13 15 17 19 23 29 31 ...... Note this partial list are all prime

the sieve is very slow in the begin and uses a lot of memory ( at least 1 bit per number). So an UNO with 2000 bytes RAM can hold max  16000 numbers
(we can skip the even ones) so 0.. 32000

Rob Tillaart

Nederlandse sectie - http://arduino.cc/forum/index.php/board,77.0.html -
(Please do not PM for private consultancy)

#### odometer

#27
##### Oct 29, 2013, 05:39 pm
I didn't mean "don't divide" or "don't calculate remainders".
What I meant was, don't use the divide or modulo operators provided as part of the programming language.
Again, see: http://www.steike.com/code/bits/division/
Also see: http://gmplib.org/manual/Exact-Division.html

#### robtillaart

#28
##### Oct 29, 2013, 07:32 pmLast Edit: Oct 29, 2013, 08:29 pm by robtillaart Reason: 1
Quote
Again, see: http://www.steike.com/code/bits/division/

This trick is well known,
If I understand it correctly one needs to know the inverse of the divider in advance. As the prime generator goes into the millions the compiler cannot generate all this numbers in advance. The program would not fit the memory. Calculating them runtime is almost as expensive as the division in itself (I assume).

Quote
Also see: http://gmplib.org/manual/Exact-Division.html

assume one test:  x % n == 0;
One could implement a quick check to see if it is a possible divider, just by checking if the last byte of x is a possible multiple of n. For this one a lookup table with the possible last bytes of multiples of n.  128 tables needed (odd numbers only) of 256 bits (=32bytes)  ==> 2^7 * 2^5 = 4096 bytes. Too much for UNO, but doable for a MEGA,

pseudo-code
p = lastbyte of X;
if (array[n/2][p] == false) ==> division not possible

fascinating...

update: a test showed that all outcomes occur when looking at the last n bits (iso last digit in decimal system). This is due to the fact that odd numbers (which are checked as dividers) are relative prime to powers of 2 (basic number theory - long ago) which is the anount of possible values for the last n bits.
see simple example for 2 bits below

Code: [Select]
`01 * 00 = 00 01 * 01 = 0101 *  10 = 1001 * 11 = 1111 * 00 = 0011 * 01 = 1111 * 10 = 1011 * 11 = 01`
All outcomes are in the last column, so it cannot be used as a selection.

In short this trick does not work for binary or hex representations, (Or convince me that I'm wrong
Rob Tillaart

Nederlandse sectie - http://arduino.cc/forum/index.php/board,77.0.html -
(Please do not PM for private consultancy)

#### nickgammon

#29
##### Oct 29, 2013, 08:48 pm

a well known method without division is the sieve
1) make a bool  array on 2..n, set them all to true

I trust you are jesting here. I have 2 kB of RAM, and I am trying to calculate a 16 decimal digit prime number. Even representing each number as a bit, they are not going to fit.
Please post technical questions on the forum, not by personal message. Thanks!