- Added a trial division test using byte-sized primes, before other testing.
- Modified the Miller-Rabin test to use 2, 3, 5, 7, 11, 13, 17, 19, and 23, rather than random numbers. This Wikipedia link - http://en.wikipedia.org/wiki/Strong_pseudoprime#Examples - says that there's no number less than 3825123056546413051, 19 digits, about 1.6 x 2
^{61}, that will fool every single-pass Miller-Rabin tests with these bases. *[Edit: Add this]*Modified loop parameters to wield only integers = 6k + 1 and 6k - 1, for integer k.

Most of the speed improvement appears to come from eliminating Miller-Rabin tests for numbers that divide by small primes. That suggests that there may be a sweet spot for the number of prime divisors in the trial division test that depends on the speed of the Miller-Rabin test, and it may be different from the one I chose. Lesser improvement comes from eliminating selection of random BigNumbers, and reducing the number of Miller-Rabin tests from ten to nine.

*[Edit: Add this]*The 6K +- 1 trick doesn't do much, since it only eliminates multiples of 3, and those fail the trial division test almost immediately.

Average displayed time was about 11.5 mS for the primes shown in replies #44 and #45. It also ran about 3.5% slower than the values shown in reply #64 - calculation times were sometimes wildly different, but on average, they were almost as good.

Here's the code, based on the sketch in reply #43, minimally changed, auto-formatted:

Code:

#include <BigNumber.h>

PROGMEM prog_uchar primes[] = {

2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31,

37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79,

83, 89, 97,101,103,107,109,113,127,131,137,

139,149,151,157,163,167,173,179,181,191,193,

197,199,211,223,227,229,233,239,241,251

};

const byte nPrimes = sizeof(primes)/sizeof(primes[0]);

const char FIRST_PRIME_TO_USE [] = "1000000000000001";

BigNumber candidate;

BigNumber one (1);

BigNumber two (2);

byte bump = 2;

void setup()

{

BigNumber five (5);

BigNumber six (6);

Serial.begin(115200);

while (!Serial) {

}

BigNumber::begin ();

Serial.println ();

Serial.println ();

Serial.println ("Starting.");

candidate = BigNumber (FIRST_PRIME_TO_USE);

BigNumber rem = candidate % six;

while ((rem != one) && (rem != five)) {

rem = rem + one;

rem = rem % six;

candidate = candidate + one;

}

if (rem == one) {

bump = 4;

}

} // end of setup

void loop()

{

for ( ; candidate < BigNumber ("9999999999999999"); candidate += bump, bump = 6 - bump)

showIfPrime (candidate);

candidate = 3; // back to start!

} // end of loop

void rng (BigNumber & result)

{

result = BigNumber (random ());

} // end of rng

bool Miller(BigNumber source)

{

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;

}

const byte bases[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23 };

const byte nbases = sizeof(bases)/sizeof(bases[0]);

for(int i = 0; i < nbases; i++)

{

// do

// {

// rng (a);

// }

// while(a < two || a >= source - two);

BigNumber a = BigNumber (bases[i]);

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 (!check(num))

return;

if (!Miller (num))

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

bool check(BigNumber n) {

for (byte i = 0;i < nPrimes; i++) {

BigNumber f = BigNumber(pgm_read_byte_near(primes + i));

if ((n % f) == 0) {

return false;

}

}

return true;

}

PROGMEM prog_uchar primes[] = {

2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31,

37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79,

83, 89, 97,101,103,107,109,113,127,131,137,

139,149,151,157,163,167,173,179,181,191,193,

197,199,211,223,227,229,233,239,241,251

};

const byte nPrimes = sizeof(primes)/sizeof(primes[0]);

const char FIRST_PRIME_TO_USE [] = "1000000000000001";

BigNumber candidate;

BigNumber one (1);

BigNumber two (2);

byte bump = 2;

void setup()

{

BigNumber five (5);

BigNumber six (6);

Serial.begin(115200);

while (!Serial) {

}

BigNumber::begin ();

Serial.println ();

Serial.println ();

Serial.println ("Starting.");

candidate = BigNumber (FIRST_PRIME_TO_USE);

BigNumber rem = candidate % six;

while ((rem != one) && (rem != five)) {

rem = rem + one;

rem = rem % six;

candidate = candidate + one;

}

if (rem == one) {

bump = 4;

}

} // end of setup

void loop()

{

for ( ; candidate < BigNumber ("9999999999999999"); candidate += bump, bump = 6 - bump)

showIfPrime (candidate);

candidate = 3; // back to start!

} // end of loop

void rng (BigNumber & result)

{

result = BigNumber (random ());

} // end of rng

bool Miller(BigNumber source)

{

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;

}

const byte bases[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23 };

const byte nbases = sizeof(bases)/sizeof(bases[0]);

for(int i = 0; i < nbases; i++)

{

// do

// {

// rng (a);

// }

// while(a < two || a >= source - two);

BigNumber a = BigNumber (bases[i]);

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 (!check(num))

return;

if (!Miller (num))

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

bool check(BigNumber n) {

for (byte i = 0;i < nPrimes; i++) {

BigNumber f = BigNumber(pgm_read_byte_near(primes + i));

if ((n % f) == 0) {

return false;

}

}

return true;

}