﻿ Irregular Bernoulli and Euler Primes
```KEYWORDS: Irregular Primes, Bernoulli numbers,
Euler numbers, computing modulo p, Akiyama-Tanigawa
algorithm, Ernst Eduard Kummer, Thomas Mautsch,
OEIS A000928, OEIS A120337, OEIS A128197.
```

```"The pioneering mathematician Kummer, over the period 1847-1850,
used his profound theory of cyclotomic fields to establish
a certain class of primes called 'regular' primes. ...
It is known that there exist an infinity of irregular primes;
in fact it is a plausible conjecture that only an asymptotic
fraction 1/Sqrt(e) ~ 0.6 of all primes are regular."

[Ribenboim, cited after OEIS A000928.]

The Computation of Irregular Primes.
(1) Bernoulli Irregular Primes

Definition:  A prime number p is called B-regular if none
of the Bernoulli numbers with even indices up to (p-3),
B_2, B_4, B_6, ..., B_(p-3), is divisible by p.

Computing the list of B-irregular primes efficiently is not

To keep things simple, I wanted to have a routine which does
fairly well in a range up to, say, n = 1000 or n = 2000.

Below I give such a routine, written in Maple. However, I did
not use any computer algebra except using a precomputed list
of primes up to n. No call to Maple's built-in function

============================================================

BIrregularPrimes := proc(len)
local t,m,j,F,P,k,lenP;

t := array(0..len): t[0] := 1;
F := {}; k := 1;
P := select(isprime, [\$2..len]); lenP := nops(P);

for m from 2 to len do

t[m - 1] := 1 / m;
for j from m-1 by - 1 to 1 do
t[j - 1] := j * (t[j - 1] - t[j]);
od;

if irem(m, 2) = 0 then next fi;
if iquo(m, 2) = P[k] then k := k+1 fi;

for j from k to lenP do
if irem(numer(t[0]),P[j]) = 0 then F := F union {P[j]} fi;
od;

od;
F end:

-----------------------------------------------------------
t:=time(): BIrregularPrimes(n): time()-t; %%; nops(%);
n=1000 ->  282 seconds,  (64 IPs),
n=2000 -> 4331 seconds, (121 IPs).
============================================================

I published this routine in the newsgroup sci.math.symbolic
and wrote:

primes. Any improvements, different implementations or algorithms
are highly appreciated."

And I got a very fine answer from Thomas Mautsch:

"I think the basic problem about your program is that you first
calculate the Bernoulli numbers directly, and Bernoulli numbers
become very quickly very large.

It might be better to do this calculation modulo p separately
for each prime number, instead of calculating huge Bernoulli
numbers and only reducing the result modulo primes afterwards."

And here is the suggestion of Thomas:

===============================================================

isBIrregular := proc(p)
local a,b,n,ah,half,quarter;

a := [1]; b:=[1];
half := iquo(p,2);
quarter := -half/2 mod p;

for n from 1 to half-1 do
ah := a[1]*quarter/n mod p;
a := [ah / (n-half) mod p, a[]];
ah - add(b[i]*a[i], i=1..n) mod p;
if % = 0 then return true else b := [b[],%] end if;
end do:

false end proc:

BIrregPrimes := proc(len) local p,F;

p := 3; F := NULL;

while p <= len do
if isBIrregular(p) then F := F,p end if;
p := nextprime(p)
end do:

F end proc:

-----------------------------------------------------------
t:=time(): BIrregPrimes(n): time()-t; %%; nops([%]);
n=1000 ->  8 seconds,  (64 IPs),
n=2000 -> 59 seconds, (121 IPs).
===============================================================

I really like the solution of Thomas, a programming pearl!
And it is fast! What a summer for good old Kummer! It took him
almost 20 years to find the first 10 or so. Thanks a lot, Thomas!

Next is a table of Bernoulli irregular primes up to len = 2000.

===============================================================
[ 37,   59,   67,  101,  103,  131, 149, 157, 233, 257, 263, 271,
283,  293,  307,  311,  347,  353, 379, 389, 401, 409, 421, 433,
461,  463,  467,  491,  523,  541, 547, 557, 577, 587, 593, 607,
613,  617,  619,  631,  647,  653, 659, 673, 677, 683, 691, 727,
751,  757,  761,  773,  797,  809, 811, 821, 827, 839, 877, 881,
887,  929,  953,  971, 1061, 1091, 1117, 1129, 1151, 1153, 1193,
1201, 1217, 1229, 1237, 1279, 1283, 1291, 1297, 1301, 1307,
1319, 1327, 1367, 1381, 1409, 1429, 1439, 1483, 1499, 1523,
1559, 1597, 1609, 1613, 1619, 1621, 1637, 1663, 1669, 1721,
1733, 1753, 1759, 1777, 1787, 1789, 1811, 1831, 1847, 1871,
1877, 1879, 1889, 1901, 1933, 1951, 1979, 1987, 1993, 1997]
================================================================

(2) Euler Irregular Primes

Definition:  A prime p is Euler-irregular (E-irregular)
iff it divides an Euler number E(2n) with 0< 2n <p-1.

The following Maple procedure computes Bernoulli irregular primes
as well as Euler irregular primes.

IrregularPrimes(100,B); # Computes Bernoulli irregular primes
IrregularPrimes(100,E); # Computes Euler irregular primes

------------------------------------------------------------
IrregularPrimes := proc(len,typ) local t,m,j,F,eg,bg,p,maxp;

eg := n -> (-1)^iquo(n-1,4)*2^iquo(1-n,2)*ceil((n mod 4)/3);
bg := n -> 1/n;

t  := array(0..len): t[0] := 1; F := {};

for m from 2 to len do
t[m - 1] := `if`(typ = B, bg(m), eg(m));
for j from m - 1 by - 1 to 1 do
t[j - 1] := j * (t[j - 1] - t[j]);
od;
if irem(m, 2) = 0 then next fi;

p := nextprime(m);
maxp := min(abs(t[0]), len);
while p <= maxp do
if abs(t[0]) mod p = 0
then print(typ[m-1], p); F := F union {p} fi;
p := nextprime(p);
od;
od;
F end:
------------------------------------------------------------

The only difference in the computation depending on the 'typ'
is the call of bg(m) and eg(m) respectively in the sixth line.
No call to Maple's built-in function bernoulli(n) or euler(n)

Next is a table of Euler irregular primes up to len = 2000.

============================================================
[19,  31,  43,  47,  61,  67,  71,  79, 101, 137, 139, 149,
193, 223, 241, 251, 263, 277, 307, 311, 349, 353, 359, 373,
379, 419, 433, 461, 463, 491, 509, 541, 563, 571, 577, 587,
619, 677, 691, 709, 739, 751, 761, 769, 773, 811, 821, 877,
887, 907, 929, 941, 967, 971, 983,1013,1019,1031,1039,1049,
1051,1069,1151,1163,1187,1223,1229,1231,1277,1279,1283,1291,
1307,1319,1361,1381,1399,1409,1423,1427,1429,1439,1447,1453,
1523,1531,1559,1583,1601,1621,1637,1663,1693,1697,1723,1733,
1759,1787,1801,1831,1867,1873,1877,1879,1889,1901,1907,1931,
1933,1951,1987,1993,1997]
============================================================

(3) Weak and Strong Irregular Primes

We say a prime is weak irregular iff
it is Bernoulli irregular or Euler irregular.

19,   31,   37,   43,   47,   59,   61,   67,   71,   79,
101,  103,  131,  137,  139,  149,  157,  193,  223,  233,
241,  251,  257,  263,  271,  277,  283,  293,  307,  311,
347,  349,  353,  359,  373,  379,  389,  401,  409,  419,
421,  433,  461,  463,  467,  491,  509,  523,  541,  547,
557,  563,  571,  577,  587,  593,  607,  613,  617,  619,
631,  647,  653,  659,  673,  677,  683,  691,  709,  727,
739,  751,  757,  761,  769,  773,  797,  809,  811,  821,
827,  839,  877,  881,  887,  907,  929,  941,  953,  967,
971,  983, 1013, 1019, 1031, 1039, 1049, 1051, 1061, 1069,
1091, 1105, 1117, 1129, 1151, 1153, 1163, 1187, 1193, 1201,
1217, 1223, 1229, 1231, 1237, 1277, 1279, 1283, 1291, 1297,
1301, 1307, 1319, 1327, 1361, 1367, 1381, 1399, 1409, 1423,
1427, 1429, 1439, 1447, 1453, 1483, 1499, 1523, 1531, 1559,
1583, 1597, 1601, 1609, 1613, 1619, 1621, 1637, 1663, 1669,
1693, 1697, 1721, 1723, 1729, 1733, 1753, 1759, 1777, 1787,
1789, 1801, 1811, 1831, 1847, 1867, 1871, 1873, 1877, 1879,
1889, 1901, 1907, 1931, 1933, 1951, 1979, 1987, 1993, 1997.

We say a prime is strong irregular iff
it is Bernoulli irregular and Euler irregular.

67,  101,  149,  263,  307,  311,  353,  379,  433,  461,
463,  491,  541,  577,  587,  619,  677,  691,  751,  761,
773,  811,  821,  877,  887,  929,  971, 1151, 1229, 1279,
1283, 1291, 1307, 1319, 1381, 1409, 1429, 1439, 1523, 1559,
1621, 1637, 1663, 1733, 1759, 1787, 1831, 1877, 1879, 1889,
1901, 1933, 1951, 1987, 1993, 1997.