For some additional excitement, I've been searching for primes $p \gg q = 104729$, where $q$ is of course the ten-thousandth prime. It seems that the best way to search for prime candidates $p$ is to do trial division $p \pmod {\ell }$ for all $\ell $ belonging to the set of all primes between $2$ and $q$, and if $p \not\equiv 0 \pmod{\ell}$, check that $2^{p-1} \equiv 1 \pmod{p}$. If so, declare that $p$ might be prime. There is obviously no point in computing $2^{p-1} \pmod{p}$ if $3 \mid p$. Competing with the time it takes to run Mathematica's PrimeQ[ ], in Python 2.73 I wrote something like:
from numpy import * primes=array([2,3,...,104729]) def trialdivandfermat(p): z=0 for q in range(0,9999): if p%int(primes[z])==0: break else: z=z+1 if z==9999: if pow(2,p-1,p)==1: return p break for p in range(1+10**1000,10**10000,2): trialdivandfermat(p)
Specifically searching for primes of the form $p = m^n \cdot h \pm 1$, $m$ odd, $h < m^n$ even, I have written that into my function trialdivandfermat(p).
Question: Are there faster ways to search for prime candidates of this form? How much trial division should I do if this is even the correct approach?