Bernoulli numbers (Erlang)

From LiteratePrograms

Other implementations: Erlang | Haskell | Python

Contents 1 Rational Numbers 2 Memoize 3 Implementation 4 Running

Rational Numbers

For generating Bernoulli numbers we first need a way to represent rational numbers. We'll make a module for handling regular number functions: add, subtract, multiply, divide, inverse, and negation.

<<rational.erl>>=
-module(rational).
-export([is_rational/1, make/2, make/1, add/2,
         sub/2, mult/2, divide/2, inv/1, neg/1]).
-record(rational, {top, bot}).
gcd(A, 0) -> A;
gcd(A, B) -> gcd(B, A rem B).
lcm(A, B) -> (A*B) div gcd(A, B).
is_rational(X) when record(X, rational) -> true;
is_rational(_) -> false.
make(Top, Bot) ->
    Gcd = gcd(Top, Bot),
    #rational{top = Top div Gcd, bot = Bot div Gcd}.
make(Num) when is_integer(Num) -> make(Num, 1);
make(Num) when record(Num, rational) -> Num.
inv(X) ->
    A = make(X),
    make(A#rational.bot, A#rational.top).
neg(X) ->
    A = make(X),
    make(-A#rational.top, A#rational.bot).
mult(X, Y) ->
    A = make(X),
    B = make(Y),
    Top = A#rational.top * B#rational.top,
    Bot = A#rational.bot * B#rational.bot,
    Gcd = gcd(Top, Bot),
    make(Top div Gcd, Bot div Gcd).
divide(X, Y) -> mult(X, inv(Y)).
add(X, Y) ->
    A = make(X),
    B = make(Y),
    Lcm = lcm(A#rational.bot, B#rational.bot),
    Top = A#rational.top * (Lcm div A#rational.bot) +
          B#rational.top * (Lcm div B#rational.bot),
    make(Top, Lcm).
sub(X, Y) -> add(X, neg(Y)).

Memoize

The values of the bernoulli numbers need to be memoized to avoid recalculation over and over of the same number. This should be transparent to the caller of this function. Ets can handle this. We check if the ets table exists and create one if it does not. We then check if the value has already been calculated and return it if it has. Otherwise we calculate it, store it in the ets table, and then return it.

<<memoize>>=
bernoulli(N) ->
    Name = bernoulli,
    case ets:info(Name) of
    undefined ->
        ets:new(Name, [public, named_table]);
    _ -> true
    end,
    case ets:lookup(Name, N) of
    [] ->
        Val = bernoulli_i(N),
        ets:insert(Name, {N, Val}),
        Val;
    [{N, Val}] -> Val;
    Else -> Else
    end.

Implementation

Now to implement the actual function. The computation can be performed more efficiently according to the following identities due to Srinivasa Ramanujan:

where

<<impl>>=
bernoulli_i(0) -> 1;
bernoulli_i(1) -> rational:make(-1,2);
bernoulli_i(N) when N rem 2 == 1 -> 0;
bernoulli_i(N) when N rem 6 == 0 ->
    Val = rational:sub(rational:make(N+3, 3), bernoulli_a(N, N div 6)),
    rational:divide(Val, choose(N+3, N));
bernoulli_i(N) when N rem 6 == 2 ->
    Val = rational:sub(rational:make(N+3, 3), bernoulli_a(N, (N-2) div 6)),
    rational:divide(Val, choose(N+3, N));
bernoulli_i(N) when N rem 6 == 4 ->
    Val = rational:sub(rational:make(-(N+3), 6), bernoulli_a(N, (N-4) div 6)),
    rational:divide(Val, choose(N+3, N)).

Now we implement A:

<<a>>=
%This can be sped up by not re-calculating the binary coefficient every time or
%memoizing that function.
bernoulli_a(N, M) -> bernoulli_a(N, M, 0).
bernoulli_a(_N, 0, Sum) -> Sum;
bernoulli_a(N, M, Sum) ->
    Val = rational:mult(choose(N+3, N-6*M), bernoulli(N-6*M)),
    bernoulli_a(N, M-1, rational:add(Sum, Val)).

The last piece we need is to calculate the binary coefficient:

<<bincoef>>=
min(X, Y) when X<Y -> X;
min(_X, Y) -> Y.
choose(N, R) ->
    choose(N, min(R, N-R), 1, 1).
choose(_N, 0, Top, Bottom) -> Top div Bottom;
choose(N, R, Top, Bottom) -> choose(N-1, R-1, Top*N, Bottom*R).

Now we put it all together in a module:

<<bernoulli.erl>>=
-module(bernoulli).
-export([bernoulli/1]).
memoize
impl
a
bincoef

Running

Now the code can be run from the erlang shell:

1> c(rational).
{ok,rational}
2> c(bernoulli).
{ok,bernoulli}
3> bernoulli:bernoulli(6).
{rational,1,42}