Bernoulli-number method wrong for input > 1

Question

I'm trying to implement a method that returns the n:th bernoulli number, like so:

Object subclass: #Bernoulli.
Bernoulli class extend [


    "******************************************************
    *   Psuedo code for bernoulli method I'm working from: 
    *                                                 
    *   function B(n)                                       
    *       B[1] <-- 1                                      
    *       for m <-- 2 to n+1 do                           
    *           B[m] <-- 0                                  
    *           for k <-- 1 to m - 1 do                     
    *               B[m] <-- B[m] − BINOM (m, k-1) * B[k]   
    *           B[m] <-- B[m]/m                               
    *       return B[n+1]                                   
    *                                                       
    *   function BINOM (n, k)                              
    *       r <-- 1                                         
    *       for i <-- 1 to k do                             
    *           r <-- r · (n − i + 1)/i                     
    *       return r                                        
    ******************************************************"
    bernoulli: n [
        "Initialize variables"
        | B size innerLoop temp|
        size := n + 1.
        B := Array new: size.

        B at: 1 put: 1.                                                             "B[1] <-- 1"
        2 to: size do: [:m |                                                        "for m <-- 2 to (n+1) do"
            B at: m put: 0.                                                             "B[m] <-- 0"

            innerLoop := m - 1.
            1 to: innerLoop do: [:k |                                                   "for k <-- 1 to (m-1) do"
                B at: m put: (B at: m) - (Bernoulli binom: m k: (k-1)) * (B at: k).         "B[m] <-- B[m] − BINOM(m, k-1) * B[k]"
            ].    

            B at: m put: (B at: m) / m.                                                 "B[m] <-- B[m] / m"
        ].
        ^(B at: size)                                                               "return B[n+1]"
    ]

    binom: n k:k [
        | r i |                    
        r := 1.                         "r <-- 1"           
        1 to: k do: [:i |               "for i <-- 1 to k do"
            r := r * (n - i + 1) / i.       "r <-- r * (n - i + 1)/i"
        ].
        ^r                              "return r"
    ]
]



z := Bernoulli bernoulli: 3.
z printNl.

(I've done my best to comment the code).

However, for inputs n > 1, I get the wrong bernoulli number:

n = 0 --> 1 (correct).
n = 1 --> -1/2 (correct)
n = 2 --> 2/3 (should be 1/6)
n = 3 --> -7/12 (should be 0)
n = 4 --> 77/45 (should be -1/30)
n = 5 --> 3157/9720 (should be 0)

My guess is that I've implemented the inner loop or inner-inner loop wrong somehow, since input n < 2 works correctly (and n < 2 skips the inner-inner loop entirely). The psuedo code I'm working from could also be incorrect, but I doubt it since I made it work in COBOL just yesterday. The binom method works correctly, I've tested this myself.

Even so, I can't figure out why this isn't working as it should. Any help is appreciated.

Leandro Caniglia · Answer 1 · 2016-05-05T15:04:19.620

3

Here is a translation of the original algorithm. The receiver is your parameter n and the answer is the array of all Bernoulli numbers from 0 to n.

Integer >> bernoulliNumbers
  | bernoulli |
  bernoulli := Array new: self + 1.
  bernoulli at: 1 put: 1.
  2 to: self + 1 do: [:m | 
    bernoulli at: m put: 0.
    1 to: m - 1 do: [:k | 
      bernoulli
        at: m
        put: (bernoulli at: m) - ((m binom: k - 1) * (bernoulli at: k))].
    bernoulli at: m put: (bernoulli at: m) / m].
  ^bernoulli

Example:

5 bernoulliNumbers -> #(1 -1/2 1/6 0 -1/30 0)

EDIT

Here is the binom: method that I created in the Integer class

binom: k
  | r |
  r := 1.
  1 to: k do: [:i | r := r * (self - i + 1) / i].
  ^r

edited May 05 '16 at 15:04

answered May 04 '16 at 14:43

Leandro Caniglia

14,495
4
29
51

This code gives me a compliation error. I'm running GNU Smalltalk 3.2.4. At least it shows me that the psuedo code is correct, but my implementation isn't working. – Mossmyr May 05 '16 at 10:00
GNU Smalltalk is not oriented towards Smalltalk development. For that you want a Smalltalk with the full tool set, especially a debugger that lets you do fix and continue – Stephan Eggermont May 05 '16 at 12:10
@Mossmyr could you please be more specific? What compilation error do you get? – Leandro Caniglia May 05 '16 at 12:25
@LeandroCaniglia `Object: Integer error: Invalid argument nil: key not found`. However, I found out why my code doesn't work. I needed parentheses around `(Bernoulli binom: m k: (k-1)) * (B at: k)` for some reason. I worked it out after looking at your implementation and comparing it to mine, but I can't mark your answer as "accepted" as it doesn't quite answer my question. – Mossmyr May 05 '16 at 14:59
@Mossmyr take a look at my *EDIT*. It should help. – Leandro Caniglia May 05 '16 at 15:07
@LeandroCaniglia Same compilation error. Please note that I just copied your code, pasted it into a file 'test.st' and then ran `gst test.st` in the terminal. – Mossmyr May 05 '16 at 15:15
1

@Mossmyr Please give it a try in Pharo first. There you will be able to understand better what goes wrong and that will make it easier for you to figure out how to improve the code for GNU-Smalltalk. – Leandro Caniglia May 05 '16 at 15:31
1

@Mossmyr you're probably not using the GNU syntax for class extension. The code in this answer works if you extend your class via `Integer extend [ ...` rather than `Integer >> ...` – lurker May 06 '16 at 16:18
@lurker You are right, tried it just now. But I am running GNU SmallTalk 3.2.4, so why can't I use the "GNU syntax"? – Mossmyr May 06 '16 at 16:58
@Mossmyr not sure what you mean. The GNU syntax works for me on 3.2.4. The GNU syntax is `Integer extend [ ...`. – lurker May 06 '16 at 16:59
@lurker I misunderstood you, I thought you said that `Integer >> ...` was the GNU syntax, but you actually said that `Integer extend [ ...` is the GNU syntax. That's why I was confused. – Mossmyr May 06 '16 at 17:17
@Mossmyr ok, that make sense. I thought that may be the case. You should accept this answer if it answers your question now. – lurker May 06 '16 at 17:19
@lurker Well, I mentioned why I havn't accepted this as my answer yet a few comments up, but I should probably elaborate. Leandro's code works well and is correct, but it doesn't answer why my own code doesn't work, which was my question. _Through examining his code_, I found out I was missing two parentheses that was causing my error, but this clarity wasn't provided by the answer itself (and took a lot of time to find, in fact). I'm new on stackoverflow, so I'm asking for clarity: If an answer isn't proper, but I found a solution through examining the answer, should I mark it as accepted? – Mossmyr May 06 '16 at 17:23
@Mossmyr I understand that thought process. If you haven't already, you could upvote this answer as being "useful". On SO you are actually allowed to answer your own question as well. You could post an answer explaining what you found was wrong and accept your own answer. – lurker May 06 '16 at 17:29

score 2 · Accepted Answer · answered May 06 '16 at 17:43

2

You (me) are missing parentheses around the (Bernoulli binom: m k: (k-1)) * (B at: k) clause. It should look like this:

B at: m put: (B at: m) - ((Bernoulli binom: m k: (k-1)) * (B at: k)).

answered May 06 '16 at 17:43

Mossmyr

909
2
10
26

2

I'm a little surprised that `*` didn't have precedence over `-` and obviate the need for those extra parentheses. But it's true! In Smalltalk, the arithmetic expression is taken left-to-right unless dictated by the parentheses. – lurker May 06 '16 at 18:21

Bernoulli-number method wrong for input > 1

2 Answers2