The List of Irregular Primes
Peter Luschny
[ http://primes.utm.edu/top20/page.php?id=26 ]
[ http://www.research.att.com/~njas/sequences/A000928 ]
To keep things simple, the routine should do fairly well
in a range up to, say, n=1000 or n=2000.
However, I do not want to use any computer algebra
except using a precomputed list of primes up to n.
Below I try to give such a routine, written for Maple.
I would like to learn more efficient ways. Any improvements,
different implementations or algorithms are highly welcome.
============================================================
IrrRegPrimes := proc(len)
local a,t,p,m,m1,j,R,F,P,r,i,q,n;
P := select(isprime, [$2..len]);
t := array(0..len): a := []; F := {};
t[0] := 1; p := 1;
for m from 1 to len do
m1 := m + 1; t[m] := 1/m1;
for j from m by - 1 to 1 do
t[j - 1] := j * (t[j - 1] - t[j ]);
od;
if irem(m1, 2) = 1 then
if irem(p, m1) = 0 and irem(denom(t[0]),m1) = 0
then a := [op(a),m1] fi;
R := {}; n := numer(t[0]);
for q in P do
if irem(n,q) = 0 then R := R union {q} fi
od;
R := R minus F;
if R <> {} then
F := F union R;
p := p * mul(i,i=R);
fi;
fi;
od;
a end:
============================================================
t:=time(); IrrRegPrimes(n); time()-t;
n=1000 -> 282 seconds,
n=2000 -> 4373 seconds.
===============================================================
[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, 1105, 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, 1729, 1733, 1753, 1759, 1777, 1787, 1789,
1811, 1831, 1847, 1871, 1877, 1879, 1889, 1901, 1933, 1951, 1979,
1987, 1993, 1997]
================================================================
Peter
Peter Luschny
> I would like to compute the list of irregular primes.
> To keep things simple, the routine should do fairly well
> in a range up to, say, n=1000 or n=2000.
> Below I try to give such a routine, written for Maple.
There is a subtle error in the routine given. A corrected
and simplified version can now be found at:
http://www.luschny.de/math/primes/irregular.html
Peter
Gottfried Helms
Do you know the site bernoulli.org of Bernd Kellner? He has
an efficient algorithm to compute bernoulli-numbers (I think
he has done up to 100 000). Two important ideas are (if I
got things right) first to cancel the systematic primefactors
of a bernoulli-number, and second to use some approximation
using zeta and the fact, that the result must be integer...
Bernd has an article about this at his page.
Regards -
Gottfried Helms
Peter Luschny
>>> Peter Luschny wrote:
>>> I would like to compute the list of irregular primes.
>>> To keep things simple, the routine should do fairly well
>>> in a range up to, say, n=1000 or n=2000.
>>> Below I try to give such a routine, written for Maple.
>>> http://www.luschny.de/math/primes/irregular.html
>> Do you know the site bernoulli.org of Bernd Kellner? He has
>> an efficient algorithm to compute bernoulli-numbers (I think
>> he has done up to 100 000).
Yes, I know. This is the kind of stuff to use when your
aim is to break some number crunching records. Whereas
my approach is minimalistic.
As I said, my aim is the range n=1..2000. I like to have
just a little program to demonstrate some basic
mathematical and/or algorithmic ideas.
>> Two important ideas are (if I
>> got things right) first to cancel the systematic primefactors
>> of a bernoulli-number, and second to use some approximation
>> using zeta and the fact, that the result must be integer...
Bernd's idea can be summed up in 6 lines Maple code.
Kellner := proc(n) local p,k,T,K,N,Z,A;
T := mul(p,p=select(isprime,{seq(k+1,k=divisors(n))})):
K := 2*n!/(2*Pi)^n;
Z := Zeta(n);
A := (-1)^(iquo(n,2)+1)*ceil(T*K*Z);
A/T end:
[Note: The above implementation is valid for even n >= 10.]
You see the difficulties. First of all the computation
of n! is involved. Not nice if n gets large. Secondly you
have to compute Zeta to a very high precision.
Kellner uses the prime factorization of n! and some
tightly bound precisions to compute the quantities.
Recently I found these slides on the Internet:
William Stein, Computing Bernoulli Numbers.
http://modular.math.washington.edu/talks/bernoulli/current.pdf
There it is suggested to compute the Zeta function as
(K and T as above)
N := ceil((K*T)^(1/(n-1)));
Z := mul(1/(1-p^(-n)),p=select(isprime,seq(i,i=2..N)));
Observing that PARI is (very much) faster than Maple
or Maxima or Mathematica 5.1 Stein posed the question:
"Do you know who came up with or implemented the idea
in PARI for computing Bernoulli numbers quickly by
approximating the zeta function and using Classen
and von Staudt’s identification of the denominator
of the Bernoulli number?"
In my opinion the answer is: This is Satz 2.7.13 in
the diploma thesis of Bernd C. Kellner 2002.
"Uber irregul"are Paare h"oherer Ordnungen."
http://www.bernoulli.org/~bk/irrpairord.pdf
In the slides someone (Kevin?? Karim??) answers:
"Henri [Cohen?] did, and wrote the initial implementation.
I wrote the current one (same idea, faster details).
The idea independently came up (Bill Daly) on pari-dev
as a speed up to Euler-Mac Laurin formulae for zeta or
gamma/loggamma."
Gruesse Peter
Thomas Mautsch
> Peter Luschny wrote:
>
>> I would like to compute the list of irregular primes.
>> To keep things simple, the routine should do fairly well
>> in a range up to, say, n=1000 or n=2000.
Maybe you would have gotten more suggestions
if you had given the definition of (ir)regular primes
instead of the links:
A prime number p is called 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;
and the Bernoulli numbers are defined by the series expension
x / ( exp(x)-1 ) = sum( B_n/n! * x^n, n=0..infinity ).
>> Below I try to give such a routine, written for Maple.
>
> There is a subtle error in the routine given. A corrected
> and simplified version can now be found at:
> http://www.luschny.de/math/primes/irregular.html
I think the basic problem about both your programs
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 afterwords.
And, of course, you only need to calculate
Bernoulli numbers with even index.
>>> I would like to learn more efficient ways. Any improvements,
>>> different implementations or algorithms are highly welcome.
Here is a suggestion:
isIrregular := proc(p) global a,b; local i,n;
for i from 2 to p-1 do a[i] := a[i-1]/i mod p end do:
for n from 1 to (p-3)/2 do
a[n+n] := a[n+n-1]/(n+n) mod p;
b[n] := a[n+n]/2 - add( b[i]*a[n+n-i-i+1],i=0..n-1) mod p;
if % = 0 then return true end if;
a[n+n+1] := a[n+n]/(n+n+1) mod p;
end do:
false
end proc:
IrregPrimes := proc(len) global a,b; local p,F;
a := array(1..len-1); a[1] := 1;
b := array(0..iquo(len-3,2)); b[0] := 1;
p := 3; F := NULL;
while p <= len do
if isIrregular(p) then F := F,p end if;
p := nextprime(p)
end do:
F
end proc:
The test
t:=time(): IrregPrimes(2000): time()-t; %%;
takes a couple of minutes on my computer,
but not nearly half an hour.
Have fun!