finding irregular primes below a limit with PARI/gp
Irregular primes were defined/known by Kummer, around 1850. Kummer studied cyclotomic extensions of Q and, I think, their ring of integers (in modern parlance).
The series of commands below has PARI/gp find the irregular primes (the number of them) below 125,000; this and more was done in a paper of Wagstaff in 1978. He found:
3559 irregular primes of index 1,
875 irregular primes of index 2,
153 irregular primes of index 3,
16 irregular primes of index 4 and
2 irregular primes of index 5.
= 4605 irregular primes up to 125,000 (Wagstaff).
Using PARI/gp, Bernoulli numbers (as provided by berfrac()),
%126 = 4605 // PARI-gp output
This agrees with Wagstaff’s computation in 1978. The method above is not the best: there are congruences known for the Bernoulli numbers. It took a few days on PARI/gp.