> Miller-Rabin primality test (40 rounds, ~10^-24 false positive rate)

It's way better than that. You are using the simplest upper bound for the false positive rate, which is 1/t^4 where t is the number of rounds. More sophisticated analysis can give better bounds.

See the paper "Average Case Error Estimates for the Strong Probable Prime Test" by Ivan Damgård, Peter Landrock and Carl Pomerance, in Mathematics of Computation Vol. 61, No. 203, Special Issue Dedicated to Derrick Henry Lehmer (Jul., 1993), pp. 177-194. Here's a PDF: https://math.dartmouth.edu/~carlp/PDF/paper88.pdf

Here are the bounds given there for t rounds testing a candidate of k bits. I'll give them as Mathematica function definitions because I happen to have them in a Mathematica notebook.

1. This one is valid for k >= 2.

  p1[k_, t_] := k^2 4^(2 - Sqrt[k])

Note this one does not depend on t, and for small k does not give a very useful bound. For 160 the bound is 0.00992742. For large k the story is different. Testing an 8192 bit number this gives a bound of 3.45661 x 10^-46. That's good enough for almost all applications so in most cases if you want an 8192 bit prime one round is good enough.

2. This one is for t = 2, k >= 88 or 3 <= t <= k/9, k >= 21.

  p2[k_, t_] := k^(3/2) 2^t t^(-1/2)  4^(2 - Sqrt[t k])

For k = 160 this is valid for 2 <= t <= 17. For t = 17 it gives 4.1 x 10^-23.

3. This one is for t >= k/9, k >= 21.

  p3[k_, t_] := 
 7/20 2^(-5 t) + 1/7 k^(15/4) 2^(-k/2 - 2 t) + 12 k 2^(-k/4 - 3 t)

For k = 160 this is valid for t >= 18. At 18 it gives 9.75 x 10^-26. At 40 it gives 1.80 x 10^-41.

4. This one is for t >= k/4, k >= 21.

  p4[k_, t_] := 1/7 k^(15/4) 2^(-k/2 - 2 t)

For k = 160 this is valid for k >= 40. At 40 it gives the same bound as p3.

So bottom line is that with your current 40 rounds your false positive rate is under 1.80 x 10^-41, considerably better than 10^-24.

If 10^-24 is an acceptable rate for this application, 18 rounds is sufficient giving a rate under 9.7 x 10^-25.

BTW, the larger the k the lower the rate. I've often seen people looking for 1024+ bit primes doing 64 or more rounds. The simplest 1/4^t bound gives 2.9 x 10^-39. OpenSSL for example does 64 for k up to 2048, and 128 for larger k.

For k = 1024 a mere 6 rounds beats that with a bound of 8.8 x 10^-41.

For k = 2048 it only takes 3 rounds to get 4.4 x 10^-41.

For k = 4096 a mere 2 sounds gives 3.8 x 10^-48.

If we had a population of 1 trillion people, each using 1000 things that needed a 4096 bit prime, and that frequently rekeyed so they needed 1000 new primes per second, and every star in the observable universe also had such a civilization consuming 4096 bit primes at that rate, and they were all using 2 rounds of Miller-Rabin, there would be around 24 false positives a year across the whole universe.

If everyone upped it to 3 rounds there would, across the whole universe, be a false positive approximately every 44 billion years.