Pages: 1 2 [3] 4 5 ... 7   Go Down
Author Topic: Arduino generating prime numbers  (Read 9389 times)
0 Members and 1 Guest are viewing this topic.
Global Moderator
Netherlands
Offline Offline
Shannon Member
*****
Karma: 216
Posts: 13664
In theory there is no difference between theory and practice, however in practice there are many...
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

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:
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)
Logged

Rob Tillaart

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

Global Moderator
Netherlands
Offline Offline
Shannon Member
*****
Karma: 216
Posts: 13664
In theory there is no difference between theory and practice, however in practice there are many...
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

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:
u64 % u64: 414200
u64 % u32: 409600
u64 % u16: 409412
u64 % u8 : 416396
u32 % u32 : 42276
u32 % u16 : 41964
u32 % u8 : 43468
u16 % u16 : 16492
u16 % u8 : 16564
u8 % 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)
//    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.

Please give it a try

(not investigated if extended to byte is faster)

update: a similar trick speeds up division: u64 /  u16 from 408 millis to 91 millis
Code:
//
//    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()
{
}
« Last Edit: October 30, 2013, 02:59:12 pm by robtillaart » Logged

Rob Tillaart

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

Grand Blanc, MI, USA
Online Online
Faraday Member
**
Karma: 95
Posts: 4058
CODE is a mass noun and should not be used in the plural or with an indefinite article.
View Profile
WWW
 Bigger Bigger  Smaller Smaller  Reset Reset

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

MCP79411/12 RTC ... "One Million Ohms" ATtiny kit ... available at http://www.tindie.com/stores/JChristensen/

Global Moderator
Offline Offline
Brattain Member
*****
Karma: 480
Posts: 18721
Lua rocks!
View Profile
WWW
 Bigger Bigger  Smaller Smaller  Reset Reset

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:
(sqrt(2^32) ) / 4 = 16384 tests

Testing up to the target prime of 100000000063 requires:

Code:
(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:
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) smiley-razz


Logged


Global Moderator
Netherlands
Offline Offline
Shannon Member
*****
Karma: 216
Posts: 13664
In theory there is no difference between theory and practice, however in practice there are many...
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

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.
« Last Edit: October 30, 2013, 05:31:12 pm by robtillaart » Logged

Rob Tillaart

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

Global Moderator
Netherlands
Offline Offline
Shannon Member
*****
Karma: 216
Posts: 13664
In theory there is no difference between theory and practice, however in practice there are many...
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset


some first numbers of the pudding test




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.

Seems not to be the predicted 6 seconds in practice, but worth the diff??
Logged

Rob Tillaart

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

Global Moderator
Offline Offline
Brattain Member
*****
Karma: 480
Posts: 18721
Lua rocks!
View Profile
WWW
 Bigger Bigger  Smaller Smaller  Reset Reset

It says "Found 5 primes" because I was compensating for not starting at 2. smiley

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


Global Moderator
Offline Offline
Brattain Member
*****
Karma: 480
Posts: 18721
Lua rocks!
View Profile
WWW
 Bigger Bigger  Smaller Smaller  Reset Reset

I tried the Miller-Rabin Primality Test from this page:

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

Code:
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:
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:
               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:
int modulo(long a,long b,long c){

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

Code:
Starting.
Logged


Global Moderator
Netherlands
Offline Offline
Shannon Member
*****
Karma: 216
Posts: 13664
In theory there is no difference between theory and practice, however in practice there are many...
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

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
I tried the Miller-Rabin Primality Test from this page:
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:
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.


to be continued...
« Last Edit: October 31, 2013, 04:12:50 pm by robtillaart » Logged

Rob Tillaart

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

Offline Offline
Full Member
***
Karma: 2
Posts: 197
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

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

Global Moderator
Offline Offline
Brattain Member
*****
Karma: 480
Posts: 18721
Lua rocks!
View Profile
WWW
 Bigger Bigger  Smaller Smaller  Reset Reset

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


Global Moderator
Netherlands
Offline Offline
Shannon Member
*****
Karma: 216
Posts: 13664
In theory there is no difference between theory and practice, however in practice there are many...
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

last results (no big step but a number of small ones) leading to these times

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.


modified code
Code:
// 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 = 1;  // 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 = 3;
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()
{
  byte cnt = 0;
  for ( ; candidate < 9999999999999999LL; candidate += 2)
  {
    showIfPrime(candidate);
  }

  candidate = 3;  // back to start!
}  // end of loop


void 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 showIfPrime


uint32_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;
}
Logged

Rob Tillaart

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

Global Moderator
Netherlands
Offline Offline
Shannon Member
*****
Karma: 216
Posts: 13664
In theory there is no difference between theory and practice, however in practice there are many...
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

@odometer,
sketch to generate the table   (UPDATE: fixed last number)
Code:
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)
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, };
« Last Edit: November 01, 2013, 02:27:12 pm by robtillaart » Logged

Rob Tillaart

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

Global Moderator
Offline Offline
Brattain Member
*****
Karma: 480
Posts: 18721
Lua rocks!
View Profile
WWW
 Bigger Bigger  Smaller Smaller  Reset Reset

I've made some progress.

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

Code:
#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 setup

void loop()
{

  for ( ; candidate < BigNumber ("9999999999999999"); candidate += 2)
    showIfPrime (candidate);

  candidate = 3;  // back to start!
}  // end of loop

void rng (BigNumber & result)
{
  result = BigNumber (random ());  
}  // end of rng

bool 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 Miller

const byte DIGITS = 16;

unsigned long start;
unsigned long lastShown;
long found = 1; // Number we found

void 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:
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:
š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:
              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.
Logged


Global Moderator
Offline Offline
Brattain Member
*****
Karma: 480
Posts: 18721
Lua rocks!
View Profile
WWW
 Bigger Bigger  Smaller Smaller  Reset Reset

And since I have a 16-digit LED display, starting at 1000000000000001 gives:

Code:
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.
Logged


Pages: 1 2 [3] 4 5 ... 7   Go Up
Jump to: