I decided to have a go at encoding the
Miller-Rabin primality test in EvalDraw.
It's a probabilistic test: it can tell you a number is either composite or "probably prime".
While the floating-point math is a bit overkill for the test, the graphing feature provides a very useful way of seeing how reliable it is, for any given prime and witness.
For comparison, I've provided a simple deterministic test for when the witness is less than 2, although it's easy to see from a graphical representation which numbers are prime and which aren't.
Also, because it only needed one extra line of code, I've also set it to show up (in a lighter colour) which numbers the
Fermat test finds as probably prime. The Miller-Rabin test is directly superior to the Fermat test, and Fermat can't detect any composites that aren't detected by Miller-Rabin.
(x,y)
{
// Press Home for a 1:1 grid view
if(keystatus[0xc7]) setgrid(-xres/2,yres/2,xres/2,-yres/2);
if( x < 0 ) return 0;
if(y<2||x<2)
// Run simple trial-factoring test
return isprime(x) + 0.3;
else
// Run Miller-Rabin test on number x, using witness y
return (millrab(x, y)) * 0.5 + 0.5;
}
millrab(n, a)
{
// Miller-Rabin probabilistic primality test
n = abs(n);
if(n<3) return (n>=2); // Handle n<3 separately
n = int(n); a = int(abs(a));
if(a%n==0) return 1; // Can't tell, return probably prime
m = n-1;
b = 2^floor(log(m,2));
p = 1;
while(m)
{
p=(p*p)%n;
if(m >= b)
{
m -= b;
p = (p*a)%n;
}
b/=2;
}
if(p==1) return 1; //probably prime
while(b>=1)
{
if(p==n-1) return 1; //probably prime
p=(p*p)%n;
b/=2;
}
if(p==1) return 0.5; //composite, but passes the Fermat test
return 0; //definitely composite
}
isprime(n)
{
// Simple deterministic trial-factoring test
n = abs(n);
if(n < 4) return (n >= 2);
if(n%2 < 1) return 0;
if(n%3 < 1) return 0;
rootn = sqrt(n);
i = 5; di = 2;
while(i <= rootn)
{
if(n%i < 1) return 0;
i += di; di = 6 - di;
}
return 1;
}