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 a trivial task. 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 bernoulli(n) is made. ============================================================ 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: "I would like to learn more efficient ways to compute irregular 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) is made. 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. External links: [1] http://primes.utm.edu/top20/page.php?id=26 [2] http://primes.utm.edu/top20/page.php?id=25 [3] http://www.ams.org/mcom/2007-76-257/S0025-5718-06-01887-4/home.html
[Last edited: 2007-01-24]