Download code

Jump to: navigation, search

Back to Bernoulli_numbers_(Erlang)

Download for Windows: zip

Download for UNIX: zip, tar.gz, tar.bz2

rational.erl

 1 % The authors of this work have released all rights to it and placed it
 2 % in the public domain under the Creative Commons CC0 1.0 waiver
 3 % (http://creativecommons.org/publicdomain/zero/1.0/).
 4 % 
 5 % THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 6 % EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 7 % MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
 8 % IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
 9 % CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
10 % TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
11 % SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
12 % 
13 % Retrieved from: http://en.literateprograms.org/Bernoulli_numbers_(Erlang)?oldid=10283
14 
15 -module(rational).
16 -export([is_rational/1, make/2, make/1, add/2,
17          sub/2, mult/2, divide/2, inv/1, neg/1]).
18 
19 -record(rational, {top, bot}).
20 
21 gcd(A, 0) -> A;
22 gcd(A, B) -> gcd(B, A rem B).
23 
24 lcm(A, B) -> (A*B) div gcd(A, B).
25 
26 is_rational(X) when record(X, rational) -> true;
27 is_rational(_) -> false.
28 
29 make(Top, Bot) ->
30     Gcd = gcd(Top, Bot),
31     #rational{top = Top div Gcd, bot = Bot div Gcd}.
32 
33 make(Num) when is_integer(Num) -> make(Num, 1);
34 make(Num) when record(Num, rational) -> Num.
35 
36 inv(X) ->
37     A = make(X),
38     make(A#rational.bot, A#rational.top).
39 
40 neg(X) ->
41     A = make(X),
42     make(-A#rational.top, A#rational.bot).
43 
44 mult(X, Y) ->
45     A = make(X),
46     B = make(Y),
47     Top = A#rational.top * B#rational.top,
48     Bot = A#rational.bot * B#rational.bot,
49     Gcd = gcd(Top, Bot),
50     make(Top div Gcd, Bot div Gcd).
51 
52 divide(X, Y) -> mult(X, inv(Y)).
53 
54 add(X, Y) ->
55     A = make(X),
56     B = make(Y),
57     Lcm = lcm(A#rational.bot, B#rational.bot),
58     Top = A#rational.top * (Lcm div A#rational.bot) +
59           B#rational.top * (Lcm div B#rational.bot),
60     make(Top, Lcm).
61 
62 sub(X, Y) -> add(X, neg(Y)).
63 


hijacker
hijacker
hijacker
hijacker

bernoulli.erl

 1 % The authors of this work have released all rights to it and placed it
 2 % in the public domain under the Creative Commons CC0 1.0 waiver
 3 % (http://creativecommons.org/publicdomain/zero/1.0/).
 4 % 
 5 % THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 6 % EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 7 % MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
 8 % IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
 9 % CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
10 % TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
11 % SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
12 % 
13 % Retrieved from: http://en.literateprograms.org/Bernoulli_numbers_(Erlang)?oldid=10283
14 
15 -module(bernoulli).
16 -export([bernoulli/1]).
17 
18 bernoulli(N) ->
19     Name = bernoulli,
20     case ets:info(Name) of
21     undefined ->
22         ets:new(Name, [public, named_table]);
23     _ -> true
24     end,
25     case ets:lookup(Name, N) of
26     [] ->
27         Val = bernoulli_i(N),
28         ets:insert(Name, {N, Val}),
29         Val;
30     [{N, Val}] -> Val;
31     Else -> Else
32     end.
33 bernoulli_i(0) -> 1;
34 bernoulli_i(1) -> rational:make(-1,2);
35 bernoulli_i(N) when N rem 2 == 1 -> 0;
36 bernoulli_i(N) when N rem 6 == 0 ->
37     Val = rational:sub(rational:make(N+3, 3), bernoulli_a(N, N div 6)),
38     rational:divide(Val, choose(N+3, N));
39 bernoulli_i(N) when N rem 6 == 2 ->
40     Val = rational:sub(rational:make(N+3, 3), bernoulli_a(N, (N-2) div 6)),
41     rational:divide(Val, choose(N+3, N));
42 bernoulli_i(N) when N rem 6 == 4 ->
43     Val = rational:sub(rational:make(-(N+3), 6), bernoulli_a(N, (N-4) div 6)),
44     rational:divide(Val, choose(N+3, N)).
45 %This can be sped up by not re-calculating the binary coefficient every time or
46 %memoizing that function.
47 bernoulli_a(N, M) -> bernoulli_a(N, M, 0).
48 
49 bernoulli_a(_N, 0, Sum) -> Sum;
50 bernoulli_a(N, M, Sum) ->
51     Val = rational:mult(choose(N+3, N-6*M), bernoulli(N-6*M)),
52     bernoulli_a(N, M-1, rational:add(Sum, Val)).
53 min(X, Y) when X<Y -> X;
54 min(_X, Y) -> Y.
55 
56 choose(N, R) ->
57     choose(N, min(R, N-R), 1, 1).
58 choose(_N, 0, Top, Bottom) -> Top div Bottom;
59 choose(N, R, Top, Bottom) -> choose(N-1, R-1, Top*N, Bottom*R).


hijacker
hijacker
hijacker
hijacker