## 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.

3559+875+153+16+2

= 4605 irregular primes up to 125,000 (Wagstaff).

Using PARI/gp, Bernoulli numbers (as provided by berfrac()),

k=125000;

g2=primes(k);

np=primepi(k);

w3=vector(np);

for(Z=6,k/2,a=2*Z;b=numerator(abs(bernfrac(a)));for(X=1,np,if(g2[X]>a+2,if(lift(Mod(b,g2[X]))<1,w3[X]=1))));

sum(X=1,np,w3[X])

%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.