💾 Archived View for tilde.team › ~steve › blog › prime.gmi captured on 2023-07-22 at 16:51:03. Gemini links have been rewritten to link to archived content
⬅️ Previous capture (2021-11-30)
-=-=-=-=-=-=-
Recently I wrote about the RSA algorithm and showed some python code to experiment with that. The starting point of RSA is having 2 large prime numbers, which I took from a website. I was interested in how one can generate these big numbers. I came across the Miller-Rabin primality test.
It's a probabilistic test - it can have false positives but not false negatives. There is a parameter k such that an inner loop does k iterations. Complexity is O(k log^3 n) where n is the input number and error rate is at most 4^(-k).
Here is a simple python implementation. The idea is pick a random number with d bits using the python function random.getrandbits and then test it using the Miller-Rabin test. Choosing k > 10 should be good enough.
import random def prime_test(n,k): r = 0 n2 = n - 1 while n2 & 1 == 0: r += 1 n2 = n2 >> 1 d = n2 << 1 for _ in range(k): a = random.randint(2, n - 2) x = pow(a, d, n) if x == 1 or x == n - 1: continue for _ in range(r-1): x = pow(x,2,n) if x == n - 1: break return "Composite" return "Probably Prime"
You run the test like this:
bit=500 while True: num = random.getrandbits(bit) if prime_test(num, 20) == 'Probably Prime': print(num) break
On my laptop it takes about 5ms for 100bit number, 200ms for 500bit and a few seconds for 1000bit. I think that a C implementation could be at least 10x faster or more. Maybe I'll try implement it in Julia, for fun.
-------------------------------------
As I thought, naive implementation in Julia runs about 10x faster, when comparing the same amount of "while True" iterations. I never coded in Julia before; it's really like Python, but without the ":" at the end of lines and with "end" at the end of blocks. pow is powermod, and I had to improvise when randomly choosing large integers. In Fact Python is very nice in that regard: integers are actually arbitrary-precision; in Julia they have BigInt but in Python it's transparent. In C, C++ I would have had to write it myself.