One last tweak. This one includes a list of numbers that fool the Miller-Rabin test for bases 2, 3, 5, and 7. It compares each odd number to the next item in the list, and skips it if it finds them equal. The previous version, in reply #66, uses nine bases, as required to eliminate pseudoprimes below 1016. This one uses four and a table lookup, so it's a lot faster. It kind of feels like cheating, though - Nick's versions knew that 2 an 3 are prime, while this program has a lot of arcane knowledge about special numbers.
I was surprised to see BigNumber work so quickly compared to native number formats, as I'd expect it to carry substantial overhead. I suspect that its progenitors were developed, at least in part, for dealing with calculations of primes, and that may account for its fitness for this use.
My thanks to Nick for posting this project, and for coming up with the BigNumber library. Hunting for primes is easy; speeding up the process is really tricky. And, I learned how to spell, "pseudo-," I think, so my ability to eruditely describe things that aren't quite what they seem is much improved.
A list of strong psuedo fake primes to bases 2 and 3 is here - http://www.bell-labs.com/user/bleichen/diss/thesis.html. The CPAN module Math::Primality provided the engine for extracting the pseudoprimes that penetrate tests for bases 5 and 7. The included file containing the byte-wide bits of the list is attached - "Pseudoprimes.h".
The code has its oddities. Because there seems to be no way to store or access uint64_t values with PROGMEM, and because I couldn't get BigNumber to play nice with variables bigger than one byte, pseudoprimes are stored as 8 bytes, big-endian, and decoded into BigNumbers. Here's the sketch:
#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]);
#include "Pseudoprimes.h";
const uint16_t nPSP = (sizeof(psps)/sizeof(psps[0])) >> 3;
BigNumber candidate;
BigNumber one (1);
BigNumber two (2);
unsigned long start;
long found = 0; // Number of primes found
void setup()
{
Serial.begin(115200);
Serial.println ("Starting.");
BigNumber::begin ();
candidate = BigNumber ("100000000003");
} // end of setup
void loop()
{
BigNumber nextPSP;
int nextPSPindex;
for (int i=0;i<nPSP;i++) {
nextPSP = fetchPSP(i);
if ((nextPSP == candidate) || (nextPSP > candidate)) {
nextPSPindex = i;
break;
}
nextPSP = 0;
nextPSPindex = 0;
}
start = millis();
for ( ; candidate < BigNumber("9999999999999999"); candidate += 2) {
if (nextPSP != candidate) {
showIfPrime (candidate);
}
else {
nextPSPindex++;
if (nextPSPindex < nPSP) {
nextPSP = fetchPSP(nextPSPindex);
}
}
}
candidate = 3; // back to start!
} // end of loop
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;
}
BigNumber a;
const byte bases[] = {2, 3, 5, 7,};
const byte nbases = sizeof(bases)/sizeof(bases[0]);
for(int i = 0; i < nbases; i++)
{
a = BigNumber (bases[i]);
BigNumber x = a.powMod (d, source);
if(x == one || x == source - one)
continue;
bool stillOK = false;
for(int r = 1; r < s; r++)
{
x = x.powMod(two, source);
if(x == one)
return false;
if(x == source - one) {
stillOK = true;
break;
}
}
if(!stillOK)
return false;
}
return true;
} // end of Miller
void showIfPrime (BigNumber n)
{
if (!check(n))
return;
if (!Miller (n))
return;
found++;
Serial.print(n);
Serial.print(" is prime. Took ");
Serial.print(millis () - start);
Serial.print (" mS. Found ");
Serial.print (found);
Serial.println (" primes.");
start = millis ();
} // end of showIfPrime
bool check(BigNumber n) {
BigNumber f;
for (byte i = 0;i < nPrimes; i++) {
f = BigNumber(pgm_read_byte_near(primes+i));
if ((n % f) == 0) {
return false;
}
}
return true;
}
BigNumber fetchPSP(int i) {
// Progmem array psps contains strong pseudoprimes
// 8-byte numbers, arranged as bytes, big-endian
i <<= 3;
BigNumber x256 ("256");
BigNumber n = 0;
for (byte k=0;k<8;k++) {
n = n * x256;
uint8_t b = pgm_read_byte_near(psps+i+k);
n = n + BigNumber(b);
}
return n;
}
Here's some output:
100000000003 is prime. Took 1676 mS. Found 1 primes.
100000000019 is prime. Took 2798 mS. Found 2 primes.
100000000057 is prime. Took 3539 mS. Found 3 primes.
100000000063 is prime. Took 1770 mS. Found 4 primes.
Here's the output for a Due running the same code:
100000000003 is prime. Took 86 mS. Found 1 primes.
100000000019 is prime. Took 145 mS. Found 2 primes.
100000000057 is prime. Took 189 mS. Found 3 primes.
100000000063 is prime. Took 89 mS. Found 4 primes.
[Edit: Add this] And, the Due in the 16-digit zone:
1000000000000037 is prime. Took 416 mS. Found 1 primes.
1000000000000091 is prime. Took 236 mS. Found 2 primes.
1000000000000159 is prime. Took 529 mS. Found 3 primes.
1000000000000187 is prime. Took 261 mS. Found 4 primes.
Pseudoprimes.h (28.4 KB)
