Go Down

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

#### robtillaart

#30
##### Oct 30, 2013, 09:28 am

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.

It was not a solution for you , it is an prime find algorithm that doesn't use division.

However you could use RAM to hold as much primes as possible  to test less numbers. An array of 900 primes would take 1800bytes RAM -
assume this fits your sketch you could do

Code: [Select]
bool isPrime(uint64_t number){   for (int i=0; i<900; i++)  {    if (number % primeArray[i] == 0) return false;  }  for (uint32_t  cnum = 7001 ; cnum < 0xFFFFFFF0; cnum += 6)   {    if (number % (cnum-1) == 0) return false;     if (number % (cnum+1) == 0) return false;   }return true;}
The 900th prime is 6997.  // http://primes.utm.edu/lists/small/1000.txt
to test until 6997 with the 6+-1 method would  2332 steps, now it take 900. => factor 2.5 faster  for that part.

Another optimization is coming soon (have to clean the code first)
Rob Tillaart

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

#### robtillaart

#31
##### Oct 30, 2013, 01:09 pmLast Edit: Oct 30, 2013, 08:59 pm by robtillaart Reason: 1
The % operator is one of the slowest operators, especially for uint64_t.
Here a timing table of different (unsigned) types for 1000 % operations in micros (UNO).
Code: [Select]
u64 % u64: 414200u64 % u32: 409600u64 % u16: 409412u64 % u8 : 416396u32 % u32 : 42276u32 % u16 : 41964u32 % u8 : 43468u16 % u16 : 16492u16 % u8 : 16564u8 % u8 : 16552
from the table we see
1) the biggest type is what counts.
2) 64 bit % is 10x slower than 32 bit
3) 32 bit %  is 2.5x slower than 16 bit
4) 16 and 8 bit is equally fast.

For the prime generator / testing it therefor helps to keep the code 32 bit as far as possible and only go to 64 bit after 32 bits is exhausted.

This raised the question if it was possible to improve the 64bit % math?
I found a trick that will speed up the % operation if the denominator is in the uint16_t range (0..65535) with a factor 3++

It works as follows (explained with a decimal example)

43211234 % 13  The long number 43211234 is split in 3 overlapping parts which are % separately.
=>  4321 % 13 =  05  replace 4321 with 05
=>  0512 % 13 = 05  replace 0512 with 05
=>  0534 %13 =  1

if we do this trick on a uint64_t  we can split it in three overlapping uint32_t parts.
To ensure results fit in the max denominator is an uint16_t

IN theory the timing changes from  414 uSec  to  3x 42 usec (126) which is a factor 3.3

Proof of concept code (does also contain the modulo timing above)
Code: (experimental) [Select]
//    FILE: faster_mod_uint64_t.ino//  AUTHOR: Rob Tillaart// VERSION: 0.1.00// PURPOSE: proof of concept//    DATE: 2013-10-30//     URL://// Released to the public domain//void setup() {  Serial.begin(115200);  Serial.println("Start ");  volatile uint64_t num = 12345676123123123LL;  volatile uint64_t cnum = 6666;  volatile uint64_t x;  x = num % (cnum - 1);  Serial.println(x);  uint32_t start = micros();  for (int i=0; i<1000; i++)  {    x = num % (cnum - 1);  }  uint32_t stop = micros();  Serial.print("u64 % u64:\t");  Serial.println(stop - start);  start = micros();  for (int i=0; i<1000; i++)  {    x = num % ((uint32_t)(cnum - 1));  }  stop = micros();  Serial.print("u64 % u32:\t");  Serial.println(stop - start);  start = micros();  for (int i=0; i<1000; i++)  {    x = num % ((uint16_t)(cnum - 1));  }  stop = micros();  Serial.print("u64 % u16:\t");  Serial.println(stop - start);  start = micros();  for (int i=0; i<1000; i++)  {    x = num % ((uint8_t)(cnum - 1));  }  stop = micros();  Serial.print("u64 % u8 :\t");  Serial.println(stop - start);  start = micros();  for (int i=0; i<1000; i++)  {    x = ((uint32_t) num) % ((uint32_t)(cnum-1));  }  stop = micros();  Serial.print("\nu32 % u32 :\t");  Serial.println(stop - start);  start = micros();  for (int i=0; i<1000; i++)  {    x = ((uint32_t) num) % ((uint16_t)(cnum-1));  }  stop = micros();  Serial.print("u32 % u16 :\t");  Serial.println(stop - start);  start = micros();  for (int i=0; i<1000; i++)  {    x = ((uint32_t)num) % ((uint8_t)(cnum - 1));  }  stop = micros();  Serial.print("u32 % u8 :\t");  Serial.println(stop - start);  start = micros();  for (int i=0; i<1000; i++)  {    x = ((uint16_t)num) % ((uint16_t)(cnum - 1));  }  stop = micros();  Serial.print("\nu16 % u16 :\t");  Serial.println(stop - start);  start = micros();  for (int i=0; i<1000; i++)  {    x = ((uint16_t)num) % ((uint8_t)(cnum - 1));  }  stop = micros();  Serial.print("u16 % u8 :\t");  Serial.println(stop - start);  start = micros();  for (int i=0; i<1000; i++)  {    x = ((uint8_t)num) % ((uint8_t)(cnum - 1));  }  stop = micros();  Serial.print("\nu8 % u8 :\t");  Serial.println(stop - start);  // THE TRICK NEEDS AN UNION  union  {    uint64_t ll;    struct    {      uint32_t lb;      uint32_t lh;    }    a;    struct    {      uint16_t ia;      uint32_t lm;      uint16_t ib;    }    b;  }  tt;  start = micros();  for (int i=0; i<1000; i++)  {    if (cnum < 0xFFFF)    {      tt.ll = num;      uint16_t t = cnum - 1;      if (tt.a.lh >= t) tt.a.lh %= t;      if (tt.b.lm >= t) tt.b.lm %= t;      if (tt.a.lb >= t) x = tt.a.lb % t;      else x = tt.a.lb;    }    else    {      x =  num % (cnum -1);    }      }  stop = micros();  Serial.print("\nsmart:\t");  Serial.println(stop - start);  tt.ll = num;  uint16_t t = cnum - 1;  if (tt.a.lh >= t) tt.a.lh %= t;  if (tt.b.lm >= t) tt.b.lm %= t;  if (tt.a.lb >= t) x = tt.a.lb % t;  else x = tt.a.lb;  Serial.println(x);}void loop() {}
The example sketch runs the modulo calculation in 123 usec which matches the theory quite well.

This way at least all numbers from 0..65535 can be tested faster.

(not investigated if extended to byte is faster)

update: a similar trick speeds up division: u64 /  u16 from 408 millis to 91 millis
Code: [Select]
////    FILE: faster_div_uint64_t.ino//  AUTHOR: Rob Tillaart// VERSION: 0.1.00// PURPOSE: proof of concept//    DATE: 2013-10-30//     URL://// Released to the public domain//void setup() {  Serial.begin(115200);  Serial.println("Start ");  volatile uint64_t num = 12345676123123123LL;  volatile uint64_t cnum = 6666;  volatile uint64_t x;  x = num / (cnum - 1);   Serial.println(x);  uint32_t start = micros();  for (int i=0; i<1000; i++)  {    x = num / (cnum - 1);   }  uint32_t stop = micros();  Serial.print("u64 / u64:\t");  Serial.println(stop - start);  start = micros();  for (int i=0; i<1000; i++)  {    x = num / ((uint32_t)(cnum - 1));   }  stop = micros();  Serial.print("u64 / u32:\t");  Serial.println(stop - start);  start = micros();  for (int i=0; i<1000; i++)  {    x = num / ((uint16_t)(cnum - 1));   }  stop = micros();  Serial.print("u64 / u16:\t");  Serial.println(stop - start);  start = micros();  for (int i=0; i<1000; i++)  {    x = num / ((uint8_t)(cnum - 1));   }  stop = micros();  Serial.print("u64 / u8 :\t");  Serial.println(stop - start);  start = micros();  for (int i=0; i<1000; i++)  {    x = ((uint32_t) num) / ((uint32_t)(cnum-1));  }  stop = micros();  Serial.print("\nu32 / u32 :\t");  Serial.println(stop - start);  start = micros();  for (int i=0; i<1000; i++)  {    x = ((uint32_t) num) / ((uint16_t)(cnum-1));  }  stop = micros();  Serial.print("u32 / u16 :\t");  Serial.println(stop - start);  start = micros();  for (int i=0; i<1000; i++)  {    x = ((uint32_t)num) / ((uint8_t)(cnum - 1));   }  stop = micros();  Serial.print("u32 / u8 :\t");  Serial.println(stop - start);  start = micros();  for (int i=0; i<1000; i++)  {    x = ((uint16_t)num) / ((uint16_t)(cnum - 1));   }  stop = micros();  Serial.print("\nu16 / u16 :\t");  Serial.println(stop - start);  start = micros();  for (int i=0; i<1000; i++)  {    x = ((uint16_t)num) / ((uint8_t)(cnum - 1));   }  stop = micros();  Serial.print("u16 / u8 :\t");  Serial.println(stop - start);  start = micros();  for (int i=0; i<1000; i++)  {    x = ((uint8_t)num) / ((uint8_t)(cnum - 1));   }  stop = micros();  Serial.print("\nu8 / u8 :\t");  Serial.println(stop - start);  union   {    uint64_t ll;    struct    {      uint32_t lb;      uint32_t lh;    }     a;    struct    {      uint16_t ia;      uint32_t lm;      uint16_t ib;    }     b;  }   tt;  uint64_t quot;  start = micros();  for (int i=0; i<1000; i++)  {    quot = 0;    if (cnum < 0xFFFF)    {      tt.ll = num;      uint16_t t = cnum - 1;      uint32_t q = tt.a.lh / t;      quot = q;      tt.a.lh -= q*t;      quot <<= 32;      q = tt.b.lm / t;      quot += ((uint64_t)q) << 16;      tt.b.lm -= q*t;      q = tt.a.lb / t;      quot += q;      //tt.a.lb -= q*t;    }    else    {      x =  num / (cnum -1);    }      }  stop = micros();  Serial.print("\nsmart:\t");  Serial.println(stop - start);  quot = 0;  tt.ll = num;  uint16_t t = cnum - 1;  uint32_t q = tt.a.lh / t;  tt.a.lh -= q*t;  quot = ((uint64_t)q) << 32;  q = tt.b.lm / t;  tt.b.lm -= q*t;  quot += ((uint64_t)q) << 16;  q = tt.a.lb / t;  tt.a.lb -= q*t;  quot += q;  Serial.println(quot);}void loop() {}
Rob Tillaart

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

#### Jack Christensen

#32
##### Oct 30, 2013, 04:49 pm

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

Without a delay in the code, it took nearly 7.4 days to roll over.

#### nickgammon

#33
##### Oct 30, 2013, 09:28 pm
Quote

This way at least all numbers from 0..65535 can be tested faster.

I tried a simple test of doing the calculations up to 2^32 using unsigned long and then switching to long long. It didn't make any noticeable difference (still took 27 seconds per prime, roughly).

Checking how long it takes to reach 89999999 shows that the early calculations are fast anyway (around 800 mS) so improving that 800 mS a bit doesn't change the overall result. It's the remaining 26 seconds we are worried about.

So what I think is needed is one of those primality probability tests, that we run half-a-dozen times, hopefully taking somewhat less than 27 seconds overall.

Let's put it this way, testing 2 out of 6 numbers up to 2^32 =

Code: [Select]
 (sqrt(2^32) ) / 4 = 16384 tests

Testing up to the target prime of 100000000063 requires:

Code: [Select]
 (sqrt (100000000063) ) / 4 = 79057  tests

Difference is 62673 tests. So there are still a lot needed doing with long long. Still that is only around 4 times as many tests. But we have to test more than one number (otherwise we wouldn't need any tests!). The average gap between primes is ln(p) namely:

Code: [Select]
 ln (100000000063) / 4 = 6.33

So given we need to test 6 numbers, and they take 4 times as long each, taking 27 seconds is in the ball park.

(I think I'm wrong about dividing by 4 above, but I'm waiting for someone to challenge me on that)

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

#### robtillaart

#34
##### Oct 30, 2013, 10:18 pmLast Edit: Oct 30, 2013, 11:31 pm by robtillaart Reason: 1
What my above code does is speed up a part of the primality test when you are in the long long domain for the number under test. so every long long tested will benefit from that.
Testing ..65535 will reject a lot of numbers much faster, but not all that's true.

How much?:
for every prime long long  found there will be 21845  faster tests (1/3 of 65535   using the 6k+-1 method)
savings per calculation in this range is 414-126 uSec = 288 usec
21845 * 288 usec = 6.291.360 usec  or about 6.2 seconds per prime found.  (>20% of the 27 seconds)

As there is also a gain - far more difficult to calculate - for all the non primes.
I expect this to be at least 10%, probably substantial more, as most non primes are filtered by the lower numbers.
but the proof is in the pudding test.
Rob Tillaart

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

#### robtillaart

#35
##### Oct 30, 2013, 11:33 pm

some first numbers of the pudding test

[font=Courier]

1) reference version (I removed the EEPROMWriteAnything)
======================
100000000003 is prime. Took 27744 mS. Found 5 primes.
100000000019 is prime. Took 33283 mS. Found 6 primes.
100000000057 is prime. Took 37525 mS. Found 7 primes.
100000000063 is prime. Took 32755 mS. Found 8 primes.

2) optimized division
========================
100000000003 is prime. Took 23964 mS. Found 5 primes.
100000000019 is prime. Took 29147 mS. Found 6 primes.
100000000057 is prime. Took 30536 mS. Found 7 primes.
100000000063 is prime. Took 28974 mS. Found 8 primes.

3) 32bit  cnum loop 1.0.4  (as the sqrt(64bit) is maximum a 32 bit)
===========================
100000000003 is prime. Took 23510 mS. Found 5 primes.
100000000019 is prime. Took 28685 mS. Found 6 primes.
100000000057 is prime. Took 30016 mS. Found 7 primes.
100000000063 is prime. Took 28521 mS. Found 8 primes.
100000000069 is prime. Took 28520 mS. Found 9 primes.
100000000073 is prime. Took 28520 mS. Found 10 primes.
100000000091 is prime. Took 28524 mS. Found 11 primes.
100000000103 is prime. Took 28533 mS. Found 12 primes.
100000000129 is prime. Took 28593 mS. Found 13 primes.

4) removed 2 lines delay(500)
===============================
100000000003 is prime. Took 23510 mS. Found 5 primes.
100000000019 is prime. Took 23685 mS. Found 6 primes.
100000000057 is prime. Took 25017 mS. Found 7 primes.
100000000063 is prime. Took 23520 mS. Found 8 primes.
100000000069 is prime. Took 23521 mS. Found 9 primes.
[/font]
Seems not to be the predicted 6 seconds in practice, but worth the diff??
Rob Tillaart

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

#### nickgammon

#36
##### Oct 31, 2013, 04:49 am
It says "Found 5 primes" because I was compensating for not starting at 2.

I appreciate the work you are putting into this, but to be honest an order of magnitude would excite me. That is, getting it down to 2 or so seconds.

So I really think the repeated division/modulus has to go out the window.
Please post technical questions on the forum, not by personal message. Thanks!

#### nickgammon

#37
##### Oct 31, 2013, 07:15 am

http://community.topcoder.com/tc?module=Static&d1=tutorials&d2=primalityTesting

Code: [Select]
const long long FIRST_PRIME_TO_USE = 961748941LL;  long long candidate;void setup() {  Serial.begin(115200);     while (!Serial) { }    Serial.println ();  Serial.println ();  Serial.println ("Starting.");    candidate = FIRST_PRIME_TO_USE;  }void loop()   {  for ( ; candidate < 9999999999999999LL; candidate += 2)    showIfPrime (candidate);   candidate = 3;  // back to start!  }  // end of loop/* this function calculates (a*b)%c taking into account that a*b might overflow */long long mulmod(long long a,long long b,long long c){    long long x = 0,y=a%c;    while(b > 0){        if(b%2 == 1){            x = (x+y)%c;        }        y = (y*2)%c;        b /= 2;    }    return x%c;}/* This function calculates (ab)%c */int modulo(int a,int b,int c){    long long x=1,y=a; // long long is taken to avoid overflow of intermediate results    while(b > 0){        if(b%2 == 1){            x=(x*y)%c;        }        y = (y*y)%c; // squaring the base        b /= 2;    }    return x%c;}/* Miller-Rabin primality test, iteration signifies the accuracy of the test */bool Miller(long long p,int iteration){    if(p<2){        return false;    }    if(p!=2 && p%2==0){        return false;    }    long long s=p-1;    while(s%2==0){        s/=2;    }    for(int i=0;i<iteration;i++){        long long a=rand()%(p-1)+1,temp=s;        long long mod=modulo(a,temp,p);        while(temp!=p-1 && mod!=1 && mod!=p-1){            mod=mulmod(mod,mod,p);            temp *= 2;        }        if(mod!=p-1 && temp%2==0){            return false;        }    }    return true;}const byte DIGITS = 9;void showIfPrime (long long num)   {    if (!Miller (num, 50))    return;      char buf [20] = { 0 };  long long temp = num;  byte pos = 0;    // make printable  for (pos = 0; pos < DIGITS; pos++)    {    if (temp == 0)      break;    char digit = temp % 10;    buf [DIGITS - 1 - pos] = digit | '0';    temp /= 10;    }// insert leading spaces  for (; pos < DIGITS; pos++)    buf [DIGITS - 1 - pos] = ' '; //  = 15; static int count = 0; Serial.print (' '); Serial.print (buf); if (++count >= 8)   {   Serial.println (' ');   count = 0;     }}  // end of showIfPrime

It runs quite fast (a line every second or so) so that's great, right? Except that not all of these are prime:

Code: [Select]
Starting. 961748941 961748943 961748947 961748949 961748951 961748955 961748957 961748959 961748961 961748963 961748965 961748967 961748971 961748973 961748975 961748979 961748981 961748983 961748987 961748989 961748991 961748993 961748995 961748997 961748999 961749003 961749005 961749007 961749011 961749013 961749015 961749019 961749021 961749023 961749025 961749027 961749029 961749031 961749035 961749037 961749039 961749043 961749045 961749047 961749051 961749053 961749055 961749059 961749061 961749063 961749067 961749069 961749071 961749075 961749077 961749079 961749083 961749085 961749087 961749089 961749091 961749093 961749095 961749099 961749101 961749103 961749107 961749109 961749111 961749115 961749117 961749119 961749121 961749123 961749125 961749127 961749131 961749133 961749135 961749139 961749141 961749143 961749147 961749149 961749151 961749153 961749155 961749157 961749159 961749163 961749165 961749167 961749171 961749173 961749175 961749179 961749181 961749183 961749187 961749189 961749191

For example, the second one (961748943) is not prime.

Compare to:

Code: [Select]
               The Fiftieth 1,000,000 Primes (from primes.utm.edu) 961748941 961748947 961748951 961748969 961748987 961748993 961749023 961749037 961749043 961749067 961749079 961749091 961749097 961749101 961749121 961749157 961749167 961749193 961749199 961749221 961749227 961749247 961749253 961749293

I thought that the data types in modulo were probably too small, so I made them larger, eg.

Code: [Select]
int modulo(long a,long b,long c){

Now, for no particularly obvious reason, it outputs nothing at all, ie.

Code: [Select]
Starting.
Please post technical questions on the forum, not by personal message. Thanks!

#### robtillaart

#38
##### Oct 31, 2013, 09:51 amLast Edit: Oct 31, 2013, 10:12 pm by robtillaart Reason: 1

I appreciate the work you are putting into this, but to be honest an order of magnitude would excite me. That is, getting it down to 2 or so seconds.
So I really think the repeated division/modulus has to go out the window.

Hi Nick,

I agree, but I did not found the silver bullet yet (but I like optimizing / squeezing code as you know)

I'm still thinking that the result of num %(cnum-1)  can be used for num % (cnum +1) or other way around as the results are close.
Moving from the 64 bit to 32 bit domain gains the factor of 10 6 you're looking for.

Quote

MR is a probability test, not an exact test and the reason why it is used is that it detects numbers that are easy to factorize . That is the key to modern encryption to have a number that are hard to factorize. Ideally these are the product of two large primes. In practice the number nerds found that the product of any two hard to factorize numbers work pretty well too. Stil it is a probability, not a proof for primality.

In the meantime found some more small optimizations resulting in:
[font=Courier]5) %3 in faster version
===========================
100000000003 is prime. Took 23503 mS. Found 5 primes.
100000000019 is prime. Took 23677 mS. Found 6 primes.
100000000057 is prime. Took 25001 mS. Found 7 primes.
100000000063 is prime. Took 23514 mS. Found 8 primes.
100000000069 is prime. Took 23515 mS. Found 9 primes.

6) cnum uint32_t
--------------------
100000000003 is prime. Took 23409 mS. Found 5 primes.
100000000019 is prime. Took 23584 mS. Found 6 primes.
100000000057 is prime. Took 24914 mS. Found 7 primes.
100000000063 is prime. Took 23420 mS. Found 8 primes.
100000000069 is prime. Took 23420 mS. Found 9 primes.

7) void showIfPrime (uint64_t num)
-------------------------------------
100000000003 is prime. Took 21912 mS. Found 5 primes.
100000000019 is prime. Took 22087 mS. Found 6 primes.
100000000057 is prime. Took 23417 mS. Found 7 primes.
100000000063 is prime. Took 21922 mS. Found 8 primes.
100000000069 is prime. Took 21923 mS. Found 9 primes.[/font]

to be continued...
Rob Tillaart

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

#### odometer

#39
##### Oct 31, 2013, 08:15 pm
I would have written code myself, but I was busy working on another project.

Since both of you gentlemen seem to be missing the point:

Yes, a lookup table will come in handy.
No, it will not need to have zillions of entries: 128 entries will do.

The lookup table will have entries for odd indices only: 0x01, 0x03, 0x05, and so on, up to 0xFF. For each entry, the index multiplied by the entry will end in ...01 (hexadecimal). The table will begin 0x01, 0xAB, 0xCD, 0xB7, ...

Let's try "dividing" 0x18979 by 0x107.
We begin by "dividing" 1 by 0x107. From our table we see the quotient ends in ...B7. In a sort of "mirror image" of how division is taught in school, we multiply 0x107 by 0xB7, giving us 0xBC01, which we then subtract from 0x0001, giving us 0x...FFFF4400. We drop the two zeros from the end, giving us 0x...FFFFFF44. We then divide this by 0x107. We multiply 0x44 by 0xB7, giving 0x309C, which tells us our next two hex digits are 9C. We multiply 0x107 by 0x9C, giving 0xA044, which we then subtract from 0x...FFFFFF44, giving 0x...FFFF5F00. (Our quotient so far is 0x...9CB7.) Again, we drop two zeros, this time from 0x...FFFF5F00, giving us 0x...FFFFFF5F. Now we can grab even more digits at a time. We will now grab four digits. We multiply 0xFF5F by 0x9CB7, giving 0x9C5470E9, so the four hex digits we grab are 70E9. We multiply 0x107 by 0x70E9, giving 0x73FF5F, and... I think we have enough digits for now. We now know that 1 "divided" by 0x107 equals 0x...70E99CB7.
Now, we multiply 0x18979 by 0x70E99CB7. We get 0xAD8C,0000,017F. Since we were only working to 8 hex digits of "precision", we can only trust the last 8 hex digits, so we have 0x0000017F. This is a reasonably small answer, given the size of our operands, and 0x17F times 0x107 does indeed equal 0x18979.
If we had tried dividing a different number, let's say 0x1897D, by 0x107 using this method, we would have got 0x1897D multiplied by 0x70E99CB7 equals 0xAD8D,C3A6,745B, which is obvious garbage: this tells us that dividing 0x1897D by 0x107 is itself a "garbage" operation (i.e. the result is not an integer, and therefore 0x1897D is not a multiple of 0x107).

#### nickgammon

#40
##### Oct 31, 2013, 08:33 pm
There is no rush, but I would be grateful for a sketch that demonstrates these ideas, not only demonstrates, but shows that primes around the value 100000000003 or greater can be found in a couple of seconds. Also that the found list is accurate. Don't worry about outputting to LEDs or anything, like the sketch just above, just output to serial will be fine.
Please post technical questions on the forum, not by personal message. Thanks!

#### robtillaart

#41
##### Oct 31, 2013, 09:44 pm
last results (no big step but a number of small ones) leading to these times
[font=Courier]
12)
-----------------------
100000000003 is prime. Took 17503 mS. Found 5 primes.
100000000019 is prime. Took 17649 mS. Found 6 primes.
100000000057 is prime. Took 18748 mS. Found 7 primes.
100000000063 is prime. Took 17514 mS. Found 8 primes.
100000000069 is prime. Took 17514 mS. Found 9 primes.
[/font]

modified code
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 = 1;  // 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 = 3;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() {  byte cnt = 0;  for ( ; candidate < 9999999999999999LL; candidate += 2)  {    showIfPrime(candidate);  }  candidate = 3;  // back to start!}  // end of loopvoid showIfPrime(uint64_t num) {  if (fastModU64(num, 3) == 0) return;  if (fastModU64(num, 5) == 0) return;  if (fastModU64(num, 7) == 0) return;  if (fastModU64(num, 11) == 0) return;  if (fastModU64(num, 13) == 0) return;  if (fastModU64(num, 17) == 0) return;  if (fastModU64(num, 19) == 0) return;  if (fastModU64(num, 23) == 0) return;  if (fastModU64(num, 29) == 0) return;  uint32_t upper = sqrt(num) + 6;  // 30k + i method  for (uint32_t cnum = 30; cnum <= upper+6; cnum += 30)  {    if (cnum < 0xFFE0)    {      if (fastModU64(num, cnum+1) == 0) return;      if (fastModU64(num, cnum+7) == 0) return;      if (fastModU64(num, cnum+11) == 0) return;      if (fastModU64(num, cnum+13) == 0) return;      if (fastModU64(num, cnum+17) == 0) return;      if (fastModU64(num, cnum+19) == 0) return;      if (fastModU64(num, cnum+23) == 0) return;      if (fastModU64(num, cnum+29) == 0) return;    }    else    {      if (num % (cnum + 1) == 0) return;      if (num % (cnum + 7) == 0) return;      if (num % (cnum + 11) == 0) return;      if (num % (cnum + 13) == 0) return;      if (num % (cnum + 17) == 0) return;      if (num % (cnum + 19) == 0) return;      if (num % (cnum + 23) == 0) return;      if (num % (cnum + 29) == 0) return;    }  }  // end of checking  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 showIfPrimeuint32_t fastModU64(uint64_t num, uint16_t cnum){  union  {    uint64_t ll;    struct    {      uint32_t lb;      uint32_t lh;    }    a;    struct    {      uint16_t ib;      uint32_t lm;      uint16_t ih;    }    b;  }  tt;  tt.ll = num;  uint16_t t = cnum;  if (tt.b.ih != 0) tt.a.lh %= t;  if (tt.a.lh != 0) tt.b.lm %= t;  if (tt.a.lb >= t) tt.a.lb %= t;  return tt.a.lb;}
Rob Tillaart

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

#### robtillaart

#42
##### Oct 31, 2013, 10:02 pmLast Edit: Nov 01, 2013, 08:27 pm by robtillaart Reason: 1
@odometer,
sketch to generate the table   (UPDATE: fixed last number)
Code: [Select]
void setup(){  Serial.begin(115200);  Serial.println("uint8_t divider[]  = {");  for (int i= 0x01; i< 0x100; i+=2)  {    for (int j=0; j <= 0xFF; j++)    {      if (((i * j) & 0xFF) == 0x01)      {        Serial.print(j);        Serial.print(",\t");      }    }    if (((i+1) % 20) == 0) Serial.println();  }  Serial.print("};");}void loop(){}

Code: (generated) [Select]
uint8_t divider[]  = {1, 171, 205, 183, 57, 163, 197, 239, 241, 27, 61, 167, 41, 19, 53, 223, 225, 139, 173, 151, 25, 131, 165, 207, 209, 251, 29, 135, 9, 243, 21, 191, 193, 107, 141, 119, 249, 99, 133, 175, 177, 219, 253, 103, 233, 211, 245, 159, 161, 75, 109, 87, 217, 67, 101, 143, 145, 187, 221, 71, 201, 179, 213, 127, 129, 43, 77, 55, 185, 35, 69, 111, 113, 155, 189, 39, 169, 147, 181, 95, 97, 11, 45, 23, 153, 3, 37, 79, 81, 123, 157, 7, 137, 115, 149, 63, 65, 235, 13, 247, 121, 227, 5, 47, 49, 91, 125, 231, 105, 83, 117, 31, 33, 203, 237, 215, 89, 195, 229, 15, 17, 59, 93, 199, 73, 51, 85, 255, };
Rob Tillaart

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

#### nickgammon

#43
##### Nov 01, 2013, 01:07 am

This sketch uses Big Numbers and the Miller-Rabin technique:

Code: [Select]
#include <BigNumber.h>const char FIRST_PRIME_TO_USE [] =  "100000000003";BigNumber candidate;BigNumber one (1);BigNumber two (2);void setup()  {  Serial.begin(115200);  while (!Serial) {  }  BigNumber::begin ();    Serial.println ();  Serial.println ();  Serial.println ("Starting.");  candidate = BigNumber (FIRST_PRIME_TO_USE);}  // end of setupvoid loop() {  for ( ; candidate < BigNumber ("9999999999999999"); candidate += 2)    showIfPrime (candidate);  candidate = 3;  // back to start!}  // end of loopvoid rng (BigNumber & result){  result = BigNumber (random ());  }  // end of rngbool Miller(BigNumber source, int certainty){  if(source == 2 || source == 3)    return true;  if(source < two || source % two == 0)    return false;  BigNumber d = source - one;  int s = 0;  while(d % two == 0)  {    d /= two;    s += one;  }  BigNumber a;  for(int i = 0; i < certainty; i++)  {    do    {      rng (a);    }    while(a < two || a >= source - two);    BigNumber x = a.powMod (d, source);    if(x == one || x == source - one)      continue;    for(int r = 1; r < s; r++)    {      x = x.powMod(two, source);      if(x == one)        return false;      if(x == source - one)        break;    }    if(x != source - one)      return false;  }  return true;}  // end of Millerconst byte DIGITS = 16;unsigned long start;unsigned long lastShown;long found = 1; // Number we foundvoid showIfPrime (BigNumber num) {  if (!Miller (num, 10))    return;  char buf [20] = {  0  };  long long temp = num;  byte pos = 0;  // make printable  for (pos = 0; pos < DIGITS; pos++)  {    if (temp == 0)      break;    char digit = temp % 10;    buf [DIGITS - 1 - pos] = digit | '0';    temp /= 10;  }  // insert leading spaces  for (; pos < DIGITS; pos++)    buf [DIGITS - 1 - pos] = ' '; //  = 15;  Serial.print(num);  Serial.print(" is prime. Took ");  Serial.print(millis () - start);  Serial.print (" mS. Found ");  Serial.print (found);  Serial.println (" primes.");  start = millis ();  found++;  lastShown = millis ();}  // end of showIfPrime

You need an updated BigNumber library (below) which adds powMod to the list of functions.

http://www.gammon.com.au/Arduino/BigNumber.zip

Output:

Code: [Select]
Starting.100000000003 is prime. Took 4079 mS. Found 1 primes.100000000019 is prime. Took 7162 mS. Found 2 primes.100000000057 is prime. Took 12018 mS. Found 3 primes.100000000063 is prime. Took 5243 mS. Found 4 primes.100000000069 is prime. Took 5028 mS. Found 5 primes.100000000073 is prime. Took 4550 mS. Found 6 primes.100000000091 is prime. Took 7831 mS. Found 7 primes.100000000103 is prime. Took 6512 mS. Found 8 primes.100000000129 is prime. Took 9747 mS. Found 9 primes.100000000171 is prime. Took 13154 mS. Found 10 primes.100000000183 is prime. Took 6597 mS. Found 11 primes.100000000193 is prime. Took 6241 mS. Found 12 primes.100000000211 is prime. Took 7847 mS. Found 13 primes.100000000223 is prime. Took 6688 mS. Found 14 primes.100000000237 is prime. Took 7080 mS. Found 15 primes.100000000253 is prime. Took 7617 mS. Found 16 primes.100000000283 is prime. Took 10479 mS. Found 17 primes.100000000319 is prime. Took 11952 mS. Found 18 primes.100000000363 is prime. Took 13784 mS. Found 19 primes.100000000367 is prime. Took 4891 mS. Found 20 primes.100000000379 is prime. Took 6733 mS. Found 21 primes.100000000393 is prime. Took 6937 mS. Found 22 primes.100000000403 is prime. Took 6107 mS. Found 23 primes.100000000411 is prime. Took 5687 mS. Found 24 primes.100000000417 is prime. Took 5289 mS. Found 25 primes.100000000427 is prime. Took 6158 mS. Found 26 primes.100000000447 is prime. Took 8605 mS. Found 27 primes.100000000487 is prime. Took 13129 mS. Found 28 primes.100000000519 is prime. Took 11158 mS. Found 29 primes.100000000567 is prime. Took 14525 mS. Found 30 primes.100000000579 is prime. Took 6495 mS. Found 31 primes.100000000621 is prime. Took 13315 mS. Found 32 primes.100000000631 is prime. Took 6294 mS. Found 33 primes.100000000637 is prime. Took 5332 mS. Found 34 primes.100000000669 is prime. Took 11160 mS. Found 35 primes.100000000699 is prime. Took 10780 mS. Found 36 primes.100000000703 is prime. Took 4982 mS. Found 37 primes.100000000721 is prime. Took 8077 mS. Found 38 primes.100000000739 is prime. Took 8041 mS. Found 39 primes.100000000747 is prime. Took 5810 mS. Found 40 primes.100000000801 is prime. Took 16136 mS. Found 41 primes.100000000817 is prime. Took 7469 mS. Found 42 primes.100000000819 is prime. Took 4372 mS. Found 43 primes.100000000861 is prime. Took 13540 mS. Found 44 primes.100000000901 is prime. Took 13042 mS. Found 45 primes.100000000943 is prime. Took 13636 mS. Found 46 primes.100000000951 is prime. Took 5911 mS. Found 47 primes.100000001009 is prime. Took 17561 mS. Found 48 primes.100000001099 is prime. Took 24052 mS. Found 49 primes.100000001111 is prime. Took 6616 mS. Found 50 primes.100000001149 is prime. Took 12536 mS. Found 51 primes.100000001201 is prime. Took 15608 mS. Found 52 primes.100000001209 is prime. Took 5687 mS. Found 53 primes.100000001237 is prime. Took 10185 mS. Found 54 primes.100000001239 is prime. Took 4457 mS. Found 55 primes.

Many of the primes are around the 5 to 10 second mark, so this is an improvement. I'm not certain if they are all correct (my CAS calculator verifies some of them). However starting at a lower number (961748941) and changing the output routine, gives this:

Code: [Select]
?Starting. 961748941 961748947 961748951 961748969 961748987 961748993 961749023 961749037 961749043 961749067 961749079 961749091 961749097 961749101 961749121 961749157 961749167 961749193 961749199 961749221 961749227 961749247 961749253 961749293 961749307 961749317 961749331 961749337 961749343 961749361 961749379 961749427 961749433 961749443 961749497 961749527 961749563 961749583 961749589 961749631 961749643 961749647 961749683 961749689 961749697 961749703 961749743 961749781 961749793 961749829 961749839 961749907 961749937 961749953 961750019 961750039 961750043 961750063 961750121 961750147 961750159 961750183 961750193 961750199 961750211 961750219 961750241 961750249 961750271 961750343 961750351 961750397 961750411 961750463 961750477 961750499 961750507 961750511 961750523 961750567 961750633 961750637 961750663 961750679 961750687 961750729 961750739 961750787 961750847 961750859 961750901 961750903 961750931 961750997 961751009 961751051 961751059 961751107 961751117 961751123 961751137 961751171 961751179 961751191 961751207 961751209 961751243 961751257 961751261 961751267 961751321 961751339 961751369 961751377 961751407 961751423 961751449 961751491 961751513 961751543 961751551 961751599 961751611 961751663 961751723 961751731 961751753 961751771 961751783 961751821 961751837 961751851 961751881 961751887 961751893 961751909

This compares favourably to this reference file:

Code: [Select]
               The Fiftieth 1,000,000 Primes (from primes.utm.edu) 961748941 961748947 961748951 961748969 961748987 961748993 961749023 961749037 961749043 961749067 961749079 961749091 961749097 961749101 961749121 961749157 961749167 961749193 961749199 961749221 961749227 961749247 961749253 961749293 961749307 961749317 961749331 961749337 961749343 961749361 961749379 961749427 961749433 961749443 961749497 961749527 961749563 961749583 961749589 961749631 961749643 961749647 961749683 961749689 961749697 961749703 961749743 961749781 961749793 961749829 961749839 961749907 961749937 961749953 961750019 961750039 961750043 961750063 961750121 961750147 961750159 961750183 961750193 961750199 961750211 961750219 961750241 961750249 961750271 961750343 961750351 961750397 961750411 961750463 961750477 961750499 961750507 961750511 961750523 961750567 961750633 961750637 961750663 961750679 961750687 961750729 961750739 961750787 961750847 961750859 961750901 961750903 961750931 961750997 961751009 961751051 961751059 961751107 961751117 961751123 961751137 961751171 961751179 961751191 961751207 961751209 961751243 961751257 961751261 961751267 961751321 961751339 961751369 961751377 961751407 961751423 961751449 961751491 961751513 961751543 961751551 961751599 961751611 961751663 961751723 961751731 961751753 961751771 961751783 961751821 961751837 961751851 961751881 961751887 961751893 961751909

I'm not sure if the rng() function is acceptable. I've seen some implementations generate larger random numbers by banging together smaller random strings to get a larger one. However the results seem to speak for themselves.
Please post technical questions on the forum, not by personal message. Thanks!

#### nickgammon

#44
##### Nov 01, 2013, 01:12 am
And since I have a 16-digit LED display, starting at 1000000000000001 gives:

Code: [Select]
Starting.1000000000000037 is prime. Took 24783 mS. Found 1 primes.1000000000000091 is prime. Took 32116 mS. Found 2 primes.1000000000000159 is prime. Took 38772 mS. Found 3 primes.1000000000000187 is prime. Took 20750 mS. Found 4 primes.1000000000000223 is prime. Took 24499 mS. Found 5 primes.1000000000000241 is prime. Took 16200 mS. Found 6 primes.1000000000000249 is prime. Took 11772 mS. Found 7 primes.
Please post technical questions on the forum, not by personal message. Thanks!