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.
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.
There is a subtle error in the routine given. A corrected and simplified version can now be found at:
>> 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:
>>> I would like to learn more efficient ways. Any improvements, >>> different implementations or algorithms are highly welcome.
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.
>>> 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.
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."
>> 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
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.