[math-fun] Bell number question
It looks as if B(p) == 2 (mod p) for prime p. Is this true, and are there "Bell pseudoprimes" in this respect?
On 7/16/07, David Wilson <davidwwilson@comcast.net> wrote:
It looks as if B(p) == 2 (mod p) for prime p.
Is this true, and are there "Bell pseudoprimes" in this respect?
This is the special case n = 0 of the "Touchard recurrence" B(n+p) = B(n+1) + B(n) mod p for all natural n and prime p. It is easily proved from the "umbral" identity B^n = (B+1)^n usually used to define the Bell numbers.
From this incidentally it follows that they are periodic modulo a prime: B(n + (p^p-1)/(p-1)) = B(n) mod p; it is conjectured that the above equation gives the minimum period.
This has been verified for p <= 101 by by S.S.Wagstaff --- I don't have the reference, but it should be easy to locate on the web. There is an earlier paper in Math. Comp. (circa 1987) explaining the background, which I can dig out for anybody who's interested. Fred Lunnon
On 7/16/07, David Wilson <davidwwilson@comcast.net> wrote:
It looks as if B(p) == 2 (mod p) for prime p.
Is this true, and are there "Bell pseudoprimes" in this respect?
There are none < 2500, according to a Q&D Maple search I did. It would not be too difficult to widen this range substantially. If it turns out that p is prime just when B(p) mod p = 2, you will have discovered the Wilson Theorem Jr. !! [the Sr. version being that p is prime just when (p-1)! mod p = 1] WFL
On 7/16/07, David Wilson <davidwwilson@comcast.net> wrote:
It looks as if B(p) == 2 (mod p) for prime p.
Besides the Touchard recurrence argument, it is well-known that B(n) = \sum_i S(n, i) where S(n, i) denotes second-kind Stirling numbers; and it can be shown (using a refinement of the Dobinski sum mentioned below) that for prime p, S(p, i) mod p = 0 unless i = 1,p, when S(p, i) = 1.
Is this true, and are there "Bell pseudoprimes" in this respect?
Yes indeed: in order of smallest factor, the first example is m = 21361 = 41*521 Verifying that B(m) mod m = 2 for this value of m is not entirely trivial. However, by the Chinese remainder theorem the result is equivalent to B(m) mod p = 2 for p | m; furthermore, it is a theorem that for prime p and natural q, in umbral notation B(p q) mod p = (B + 1)^q = B(q+1). [See for example Lunnon,~W.~F.~&~Pleasants,~P.~A.~B.~&~Stephens,~N.~M. Arithmetic Properties of Bell Numbers to Composite Modulus I, Acta Arith 35, pp~1-16 (1979).] Putting these together, we need only verify instead that B(522) mod 41 = B(42) mod 521 = 2 which would easily be evaluated directly under Maple 9, using bell(n) from the combinat package. But you can't stop progress: under Maple 10, bell(n) takes a stupendously long time to evaluate anything over n = 170 or so, despite its predecessor being easily good for n = 2500. Fortunately, a simple workaround via sums of Stirling numbers with(combinat); p := 522; q := 41; add(stirling2(p, i), i = 0..p) mod q; add(stirling2(q, i), i = 0..q) mod p; grinds to a conclusion in a minute (G4 Macintosh). A much faster alternative --- unfortunately requiring rather delicate estimation of precision and convergence --- would be the Dobinski infinite sum B(n) = (1/e) \sum_{i >= 0} n^i / i! --- but even this will not easily cope with B(m) above directly. Fred Lunnon
On 7/21/07, Fred lunnon <fred.lunnon@gmail.com> wrote:
Fortunately, a simple workaround via sums of Stirling numbers
with(combinat); p := 522; q := 41; add(stirling2(p, i), i = 0..p) mod q; add(stirling2(q, i), i = 0..q) mod p;
grinds to a conclusion in a minute (G4 Macintosh).
Heigh-ho --- now for the inevitable correction --- should have read with(combinat); p := 521; q := 41; add(stirling2(p+1, i), i = 0..p+1) mod q; add(stirling2(q+1, i), i = 0..q+1) mod p; WFL
Well, Wilson's Theorem Jr. is out the window. And of course, Bell pseudoprime is infinitely preferable to Wilson pseudoprime. So I remain in obscurity. At least I've discovered a primality test that is at once computationally difficult and unreliable. Would you be the discoverer of the first Bell pseudoprime, [MD]r. Lunnon? ----- Original Message ----- From: "Fred lunnon" <fred.lunnon@gmail.com> To: "math-fun" <math-fun@mailman.xmission.com> Sent: Saturday, July 21, 2007 12:32 AM Subject: Re: [math-fun] Bell number question
On 7/16/07, David Wilson <davidwwilson@comcast.net> wrote:
It looks as if B(p) == 2 (mod p) for prime p.
Besides the Touchard recurrence argument, it is well-known that
B(n) = \sum_i S(n, i)
where S(n, i) denotes second-kind Stirling numbers; and it can be shown (using a refinement of the Dobinski sum mentioned below) that for prime p,
S(p, i) mod p = 0 unless i = 1,p, when S(p, i) = 1.
Is this true, and are there "Bell pseudoprimes" in this respect?
Yes indeed: in order of smallest factor, the first example is
m = 21361 = 41*521
Verifying that B(m) mod m = 2 for this value of m is not entirely trivial. However, by the Chinese remainder theorem the result is equivalent to
B(m) mod p = 2 for p | m;
furthermore, it is a theorem that for prime p and natural q, in umbral notation
B(p q) mod p = (B + 1)^q = B(q+1).
[See for example Lunnon,~W.~F.~&~Pleasants,~P.~A.~B.~&~Stephens,~N.~M. Arithmetic Properties of Bell Numbers to Composite Modulus I, Acta Arith 35, pp~1-16 (1979).]
Putting these together, we need only verify instead that
B(522) mod 41 = B(42) mod 521 = 2
which would easily be evaluated directly under Maple 9, using bell(n) from the combinat package.
But you can't stop progress: under Maple 10, bell(n) takes a stupendously long time to evaluate anything over n = 170 or so, despite its predecessor being easily good for n = 2500. Fortunately, a simple workaround via sums of Stirling numbers
with(combinat); p := 522; q := 41; add(stirling2(p, i), i = 0..p) mod q; add(stirling2(q, i), i = 0..q) mod p;
grinds to a conclusion in a minute (G4 Macintosh).
A much faster alternative --- unfortunately requiring rather delicate estimation of precision and convergence --- would be the Dobinski infinite sum
B(n) = (1/e) \sum_{i >= 0} n^i / i!
--- but even this will not easily cope with B(m) above directly.
Fred Lunnon
_______________________________________________ math-fun mailing list math-fun@mailman.xmission.com http://mailman.xmission.com/cgi-bin/mailman/listinfo/math-fun
-- No virus found in this incoming message. Checked by AVG Free Edition. Version: 7.5.476 / Virus Database: 269.10.11/909 - Release Date: 7/20/2007 4:39 PM
On 7/22/07, David Wilson <davidwwilson@comcast.net> wrote:
Well, Wilson's Theorem Jr. is out the window.
And of course, Bell pseudoprime is infinitely preferable to Wilson pseudoprime.
So I remain in obscurity.
At least I've discovered a primality test that is at once computationally difficult and unreliable.
But not as unreliable as my typing: some gremlin introduces invisible typos into text as I copy and paste it, which only become visible after the send button is pressed ... Dobinski's sum should actually have read: B(n) = (1/e) \sum_{i >= 0} i^n / i! A Maple script implementing this method is (hopefully) appended below: it finds B(n) for n = 173 in 1/4 sec, for n = 2999 in 20 sec on a Mac G4. [The Maple 10 library function took 102 sec for n = 170, was aborted after 480 sec for n = 173; Maple 9 aborted after similar time for n = 2999.]
Would you be the discoverer of the first Bell pseudoprime, [MD]r. Lunnon?
As far as I know --- tho' I don't plan to have it engraved on my tombstone --- nor the 66,000 or so digits of the corresponding Bell number! Note that there may be a smaller one --- my program searched in order of smallest prime factor p, and found only that there are no more pseudoprimes for p < 60. I can prove that B(p^2)mod p = 3 for (odd) prime p; so the next candidate would be n = 61*67 = 4087 ... Fred Lunnon # Bell number B(n) via Dobinsky sum bell := proc (n) local i,i_0,i_1,digo,digs,beln; digo := Digits; Digits := 10; i_0 := fsolve(i*log(i) = n, i); # index max term i_1 := n + i_0; # number of terms required digs := n^2/i_0/log(10.0); # precision > log_10 of max term Digits := max(10, round(digs)); beln := add(evalf(i)^n/factorial(evalf(i)), i = 0..i_1)/exp(1.0); # floating-point evaluation! Digits := digo; # restore precision round(beln) end; tim := time(); # test n := 173; bell(n)mod n; # 2 for n prime tim := time() - tim;
I pushed the direct search for Bell/Wilson pseudoprimes up to n <= 4096, taking 4 hours --- there are none in this range. On 7/22/07, Fred lunnon <fred.lunnon@gmail.com> wrote:
Note that there may be a smaller one --- my program searched in order of smallest prime factor p, and found only that there are no more pseudoprimes for p < 60.
An over-simplification: this cofactor search program actually enumerates by quotient n/q of largest prime factor q, finding all possible such q. Anyhow, a 10 hour run took n/q up to 73, finding just one more pseudoprime: n = 69 * 16218646893090134590535390526854205539989357 (46 digits) Doesn't look as though there are too many of these things around! These searches employ a Dobinski-based Bell(n) function, also a recurrence- based Bell(n) mod m, which is faster for large n and small m. Maple and Magma implementations available on request [you interested, JJA?] Fred Lunnon
On 7/21/07, Fred lunnon <fred.lunnon@gmail.com> wrote: [...snip...]
A Maple script implementing this method is (hopefully) appended below: it finds B(n) for n = 173 in 1/4 sec, for n = 2999 in 20 sec on a Mac G4. [The Maple 10 library function took 102 sec for n = 170, was aborted after 480 sec for n = 173; Maple 9 aborted after similar time for n = 2999.]
I ran your maple code using Maple 7 (the only version I have) on an XP PC with Celeron 2.80 GHz, 504 MB Ram, and n = 173 took 0.21 seconds, and n = 2999 took 59 seconds. Probably I should have mentioned before (and you could have saved a couple of hours of search time) that (also with Maple 7) I searched up to n = 4000 using the built in Bell in the combinat package finding no Bell-pseudo-primes, and although I didn't time it, it took definitely less than 2 hours. Jim
hello, Maybe you should give a try with mathematica, the function is BellB[n] It is much faster than maple. bell(1000) takes forever to compute and almost instantly with mathematica. simon plouffe
From: "Simon Plouffe" <simon.plouffe@gmail.com>
Maybe you should give a try with mathematica, the function is BellB[n] It is much faster than maple.
bell(1000) takes forever to compute and almost instantly with mathematica.
Jack Schmidt just wrote on the sage-devel mailing list that GAP takes under a second to compute Bell(1000) - that means that SAGE would also compute it under a second because it uses GAP's Bell, http://modular.math.washington.edu/sage/doc/html/ref/module-sage.combinat.co... Alec
On 7/24/07, Alec Mihailovs <alec@mihailovs.com> wrote:
From: "Simon Plouffe" <simon.plouffe@gmail.com>
Maybe you should give a try with mathematica, the function is BellB[n] It is much faster than maple.
bell(1000) takes forever to compute and almost instantly with mathematica.
Jack Schmidt just wrote on the sage-devel mailing list that GAP takes under a second to compute Bell(1000) - that means that SAGE would also compute it under a second because it uses GAP's Bell,
http://modular.math.washington.edu/sage/doc/html/ref/module-sage.combinat.co...
Alec
Hmmm. John Cannon reported that his 2.4GHz Opteron took 0.9 sec to find B(2999) using the Magma transcription of my implementation. At a rough guess, this would scale down at most 0.05 sec for B(1000) ... WFL
From: "Fred lunnon" <fred.lunnon@gmail.com>
Hmmm. John Cannon reported that his 2.4GHz Opteron took 0.9 sec to find B(2999) using the Magma transcription of my implementation. At a rough guess, this would scale down at most 0.05 sec for B(1000) ...
Here is the related part of William Stein's post in sage-devel: -------------------Start of the (edited) quote------------------ On 7/24/07, Jack Schmidt <Jack.Schmidt.SciMath@gmail.com> wrote:
GAP takes under a second to compute Bell(1000), compared to over a minute (and going) for maple on the same computer.
Gap takes about 1 second for this on one of my test machines. Nick Alexander wrote an optimized native SAGE function that computes *all* bell numbers up to n very quickly, much more quickly than Gap computes even one! sage: time v=expnums(1001,1) # this computes *all* Bell numbers up to Bell(1000) CPU times: user 0.12 s, sys: 0.00 s, total: 0.12 s Wall time: 0.13 sage: v[1000] == bell_number(1000) # uses GAP; about a second True sage: v[990] == bell_number(990) True sage: time v=expnums(2000,1) CPU times: user 1.36 s, sys: 0.04 s, total: 1.40 s Wall time: 2.35 How do you compute BellB in Mathematica? I get In[2]:= BellB[1000] Out[2]= BellB[1000] -- William -- William Stein Associate Professor of Mathematics University of Washington http://www.williamstein.org -----------------------End of (edited) quote------------------- Alec
On 7/25/07, Alec Mihailovs <alec@mihailovs.com> wrote:
Here is the related part of William Stein's post in sage-devel:
-------------------Start of the (edited) quote------------------ ... Nick Alexander wrote an optimized native SAGE function that computes *all* bell numbers up to n very quickly, much more quickly than Gap computes even one!
Presumably then the algorithm used here is the usual recursion B^(n+1) = (B+1)^n. [In this range, recursion seems to run consistently 4x slower than Dobinski --- the apparent consistency of this factor is unexpected!].
sage: time v=expnums(1001,1) # this computes *all* Bell numbers up to Bell(1000) CPU times: user 0.12 s, sys: 0.00 s, total: 0.12 s Wall time: 0.13 ...
Since most of the time should be taken up by multi-precision arithmetic, it's hard to understand why using "native SAGE" (Python) makes any difference --- surely only the MP package driving the hardware can manage that.
How do you compute BellB in Mathematica?
I get
In[2]:= BellB[1000] Out[2]= BellB[1000]
-- William
Something missing here? It would be interesting to know which algorithms other programmers use! Fred Lunnon
On 7/24/07, James Buddenhagen <jbuddenh@gmail.com> wrote:
I ran your maple code using Maple 7 (the only version I have) on an XP PC with Celeron 2.80 GHz, 504 MB Ram, and n = 173 took 0.21 seconds, and n = 2999 took 59 seconds.
Say 3x slower than on my Mac Powerbook ...
Probably I should have mentioned before (and you could have saved a couple of hours of search time) that (also with Maple 7) I searched up to n = 4000 using the built in Bell in the combinat package finding no Bell-pseudo-primes, and although I didn't time it, it took definitely less than 2 hours.
Now that's a bit strange --- under Maple 9, the library version is noticeably slower than mine, although not badly so. Under Maple 10, it is hopelessly outclassed. Your timings suggest that under Maple 7 it is actually 6x faster! What _is_ going on here, I wonder? Why should anybody have bothered to modify combinat[bell]() at all between versions, let alone substantially degrade it? Or is there some other explanation? By the way --- I am told that there are difficulties with implementations of CAS systems on the Mac G4 and G5 processors, involving the operating system failing fully to exploit the 64-bit hardware --- as well, I suspect, software companies' understandable reluctance to invest in development for a minority architecture. The Magma version of my bell() function --- slightly revised in order to improve efficiency for small n, and portability --- apparently runs 90x faster on a state-of-the-art PC than on my Mac. Sickening, isn't it! This prompts the thought that my cofactor search program (also in Maple) would probably run faster on a PC: currently it has eliminated all n = r q, where q is any prime, and r is any natural < 74. I shan't take up bandwidth with it here; but if anybody requests the code, I'm happy to email it personally. Fred Lunnon
participants (5)
-
Alec Mihailovs -
David Wilson -
Fred lunnon -
James Buddenhagen -
Simon Plouffe