Related, there are a few cryptocurrencies that used things related to finding large primes as part of their proof of work functions. It turns out that ~8 years ago, a really fast primality test implementation could make you a lot of money. (For some period of time I was the author and maintainer of mining software for riecoin. Why, I have no idea, except that I like prime numbers.)

This article omits the number one optimization for fast primality testing: Montgomery multiplication

https://en.m.wikipedia.org/wiki/Montgomery_modular_multiplic...

It forms the basis of fast practical modular exponentiation implementations.

Niall Emmart, then in academia and now I believe at Nvidia, released a really, spectacularly fast GPU bigint library, CGBN: https://github.com/NVlabs/CGBN

It's still the fastest batch modexp implentation that I'm aware of. It's breathtaking, if I can gush in geek for a moment.

(Someday I should write up the story of how that helped me dominate production of a small cryptocurrency for about 5 years. Thanks, Niall - owe you a beer if we ever cross paths!)

Also worth noting that python includes an entirely fine modexp in the three-argument form of pow(x, y, m) --> compute x^y % m

With that, you can very easily implement a fermat or miller-rabin primality test if you want to roll your own, which is quite fun. If you don't, the gmp library provides a very nice mpz_probab_prime() function. Gmp is obviously faster, but it's hard to beat the fun of a two liner fermat test for playing with big primes.