Help implementing Miller-rabin algorithm

[Thaddy]

1. See https://en.wikipedia.org/w/index.php?title=Miller%E2%80%93Rabin_primality_test
That means that given 2,7,61 isprime  in the range ~4.8 billion scores no false positives.
2. The value of k is arbitrary, but a good estimate would be round(ln(range)).(over-values lower ranges a bit)'
This has a direct relation ship with
3. miller-rabin finds prime candidates faster than classic prime solvers. If it finds a prime candidate you can verify with classic prime like with this simple isprime:

Code: Pascal  [Select][+][-]

1. function IsPrime(n:Uint32):Boolean;inline;
2. var
3.   p:Uint32;
4. begin
5.   if n = 2 then exit(true);
6.   if (n < 2) or (n mod 2 = 0) then exit(false);
7.   p := 3;
8.   while p * p <= n do
9.   begin
10.     if n mod p = 0 then exit(false);
11.     inc(p,2);
12.   end;
13.   result:= true;
14. end;

Note that everything I wrote is limited to the range of native Uint32. For higher values you need a biginteger library, like the one by Rudy Velthuys.

With the chosen values, miller-rabin has been proven 100% correct for Cardinal (see wiki again);

[ad1mt]

I'm truly amazed at how quickly this algorithm works; I only wish I could understand how it works!

[AlexTP]

I never seen this algo, but we can easily understand how example works. It tests all ODD (n mod 2=1) numbers (variable p) from 3 to square_root(n). On each step, if n is divided by p, n is NOT prime.

This is school level.

[Thaddy]

@AlextTP
My last entry is not miller-rabin, but a naive test to see if miller-rabin is correct.
You did not read the whole thread. Miller-rabin is about probabilities, either by determinstic values or by random probing.
The deterministic values come from verification runs and in the wiki entry you will find them.
It is much faster than simple sieves or other naive prime algorithms on big primes.
Note that the current implementation I gave  is hampered by the limits of native unsigned integer types.(Uint64, Uint32)
Ideally you would implement it using a biginteger library. (As I wrote before)
But the algorithm is correct and the example is useable for primes in a range < ~4.7 billion.
With the chosen values it covers the primes in that range 100%

That is "not quite" highschool math....

[abouchez]

You have a very fast Miler-rabin algorithm implemented in mORMot RSA.
It has been validated with 2048-bit RSA keys on several platforms.
With a pre-defined list of small known primes, it has very good performance in practice.

Perhaps a bit overkill for your initial request, but worth considering as a working solution.

See https://github.com/synopse/mORMot2/blob/master/src/crypt/mormot.crypt.rsa.pas#L1400

[Thaddy]

@abouchez

More undiscovered beauty in mORMot. That is an implementation that does exactly what is needed and I suggested.
KUDOS to the mORMot team.

It also follows the same - correct miller-rabin - algorithm but on a production level for high primes.
I suggest to use my example to study HOW it works in somewhat simple code (it is not by any means "simple") and use the mORMot code in production.
for <~4,7 billion the results of both implementations are the same of course.
I suppose, but did not test, that the mORMot implementation can handle more than just ~4.7 billion.
E.g. it might do primes in the range of > 2 trillion with test a = 2, 3, 5, 7, and 11

[ad1mt]

Thanks for your help with this algorithm.
Found a candidate prime 801 digits long:  1000...0001537

[Thaddy]

But do you now know how it....works? 

[ad1mt]

But do you now know how it....works? 

No... I trusted your code!   

Thanks for your help with this... getting this code working with my big integer library found some more bugs, which had escaped my testing until now.

[srvaldez]

hello ad1mt
just in case you want to look at a different implementation of Miller-Rabin prime test, here's one in Basic https://raw.githubusercontent.com/stephanbrunker/big_integer/refs/heads/master/example_bigint_isprime.bas
I changed r = randomstring(128) to r = randomstring(416) and the for loop to 4000
it produced a pseudo prime of 1002 digits in less than a minute
please take no offense at this post, Basic is not Pascal but it's easy to understand and it's only about 270 lines of code
@Thaddy
I think it's basically the same as your code so please don't feel offended

[ad1mt]

@srvaldez... thanks for your suggestion.

please take no offense at this post, Basic is not Pascal

I don't have a problem with Basic. When I first started hobby programming at home in 1981, I used a Basic compiler on a Sinclair Spectrum to write a Draughts playing program!
That program has evolved over the decades and is now written in Free Pascal. I still sometimes work on trying to improve it, but I think I've reached the limits of ability on that project.

[srvaldez]

ad1mt
the reason that I posted that link was to present to you a different way of generating possible random primes
how are you progressing with your code?