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:
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.
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)).