Bernoulli numbers: Difference between revisions

Using a GMP thick binding available at http://www.codeforge.com/article/422541

USE GMP.Rationals, GMP.Integers, Ada.Text_IO, Ada.Strings.Fixed, Ada.Strings;

END IF;

END LOOP;

{{out}}

<pre>

{{works with|ALGOL 68G|Any - tested with release 2.8.3.win32}}

Uses the LONG LONG INT mode of Algol 68G which allows large precision integers.

# Show Bernoulli numbers B0 to B60 as rational numbers #

OD

END

{{out}}

<pre>

=={{header|AppleScript}}==

To be able to handle numbers up to B(60) and beyond, this script represents the numbers with lists whose items are the digit values — which of course requires custom math routines.

set listMathScript to getListMathScript(10) -- Script object providing custom list math routines.

end task

{{output}}

B(0) = 1 / 1

B(1) = 1 / 2

B(56) = -2479392929313226753685415739663229 / 870

B(58) = 84483613348880041862046775994036021 / 354

=={{header|Bracmat}}==

= B Bs answer indLn indexLen indexPadding

, n numberPadding p solPos solidusPos sp

& str$!answer

)

<pre>B(0)= 1/1

B(1)= 1/2

=={{header|C}}==

{{libheader|GMP}}

#include <stdlib.h>

#include <gmp.h>

return 0;

}

{{out}}

<pre>

{{libheader|Mpir.NET}}

Translation of the C implementation

using Mpir.NET;

using System;

}

}

{{out}}

<pre>

{{libheader|MathNet.Numerics}}

using System;

using System.Console;

}

}

{{out}}

<pre>

{{libheader|System.Numerics}}

Algo based on the example provided in the header of this RC page (the one from Wikipedia). <br/> Extra feature - one can override the default of 60 by supplying a suitable number on the command line. The column widths are not hard-coded, but will adapt to the widths of the items listed.

using System.Numerics;

using System.Collections.Generic;

}

}

{{out}}

Default (nothing entered on command line):

{{Works with|C++11}}

{{libheader|boost}}

* Configured with: --prefix=/Library/Developer/CommandLineTools/usr --with-gxx-include-dir=/usr/include/c++/4.2.1

* Apple LLVM version 9.1.0 (clang-902.0.39.1)

return 0;

{{out}}

<pre>

=={{header|Clojure}}==

ns test-project-intellij.core

(println q ":" (format-ans ans)))

{{out}}

<pre>

Be advised that the pseudocode algorithm specifies (j * (a[j-1] - a[j])) in the inner loop; implementing that as-is gives the wrong value (1/2) where n = 1, whereas subtracting a[j]-a[j-1] yields the correct value (B[1]=-1/2). See [http://oeis.org/A027641 the numerator list].

(loop with a = (make-array (list (1+ n)))

for m from 0 to n do

n

(numerator r)

{{out}}

{{Trans|Ruby}}

class Bernoulli

puts "B(%2i) = %*i/%i" % [i, max_width, v.numerator, v.denominator]

end

{{Trans|Python}}

Version 1: compute each number separately.

def bernoulli(n)

width = b_nums.map{ |b| b.numerator.to_s.size }.max

b_nums.each_with_index { |b,i| puts "B(%2i) = %*i/%i" % [i, width, b.numerator, b.denominator] unless b.zero? }

{{Trans|Python}}

Version 2: create faster generator to compute array of numbers once.

def bernoulli2(limit)

width = b_nums.map{ |b| b.numerator.to_s.size }.max

b_nums.each_with_index { |b,i| puts "B(%2i) = %*i/%i" % [i, width, b.numerator, b.denominator] unless b.zero? }

{{out}}

<pre>

This uses the D module from the Arithmetic/Rational task.

{{trans|Python}}

auto bernoulli(in uint n) pure nothrow /*@safe*/ {

foreach (immutable b; berns)

writefln("B(%2d) = %*d/%d", b[0], width, b[1].tupleof);

The output is exactly the same as the Python entry.

Thanks Rudy Velthuis for the [https://github.com/rvelthuis/DelphiBigNumbers Velthuis.BigRationals] library.<br>

program Bernoulli_numbers;

end;

readln;

=={{header|EchoLisp}}==

Only 'small' rationals are supported in EchoLisp, i.e numerator and demominator < 2^31. So, we create a class of 'large' rationals, supported by the bigint library, and then apply the magic formula.

(lib 'bigint) ;; lerge numbers

(lib 'gloops) ;; classes

(define-method div (Rational Rational) (lambda (r q)

(normalize (Rational (* r.a q.b) (* r.b q.a)))))

{{Output}}

;; Bernoulli numbers

;; http://rosettacode.org/wiki/Bernoulli_numbers

(B 1) → 1 / 2

=={{header|Elixir}}==

defmodule Rational do

import Kernel, except: [div: 2]

end

{{out}}

=={{header|F sharp|F#}}==

{{libheader|MathNet.Numerics.FSharp}}

open MathNet.Numerics

open System

printf ""

0

{{out}}

<pre>

=={{header|Factor}}==

One could use the "bernoulli" word from the math.extras vocabulary as follows:

[

0 1 1 "%2d : %d / %d\n" printf

58 : 84483613348880041862046775994036021 / 354

60 : -1215233140483755572040304994079820246041491 / 56786730

Alternatively a method described by Brent and Harvey (2011) in "Fast computation of Bernoulli, Tangent and Secant numbers" https://arxiv.org/pdf/1108.0286.pdf is shown.

n 1 + 0 <array> :> tab

1 1 tab set-nth

"%2d : %d / %d\n" printf

] each

It gives the same result as the native implementation, but is slightly faster.

...

=={{header|Fermat}}==

{{out}}<pre>

0 1

=={{header|FreeBASIC}}==

{{libheader|GMP}}

' compile with: fbc -s console

' uses gmp

Print :Print "hit any key to end program"

Sleep

{{out}}

<pre>B( 0) = 1/1

=={{header|Frink}}==

{

a = new array

}

{{out}}

=={{header|FunL}}==

FunL has pre-defined function <code>B</code> in module <code>integers</code>, which is defined as:

def B( n ) = sum( 1/(k + 1)*sum((if 2|r then 1 else -1)*choose(k, r)*(r^n) | r <- 0..k) | k <- 0..n )

for i <- 0..60 if i == 1 or 2|i

{{out}}

=={{header|GAP}}==

Print(a, "\n");

od;

[ 56, -2479392929313226753685415739663229/870 ]

[ 58, 84483613348880041862046775994036021/354 ]

=={{header|Go}}==

import (

}

}

{{out}}

<pre>

The implementation of the algorithm is in the function bernoullis. The rest is for printing the results.

import System.Environment

berno i = 1 % (i + 1)

ulli _ [_] = []

{{Out}}

<pre>B(0) = 1/1

====Derivation from Faulhaber's triangle====

import Data.Ratio (Ratio, denominator, numerator, (%))

rjust :: Int -> a -> [a] -> [a]

{{Out}}

<pre>Bernouillis from Faulhaber triangle:

The following works in both languages:

procedure main(args)

procedure align(r,n)

return repl(" ",n-find("/",r))||r

Sample run:

See [[j:Essays/Bernoulli_Numbers|Bernoulli Numbers Essay]] on the J wiki.

'''Task:'''

B 0 1

B 1 -1/2

B56 -2479392929313226753685415739663229/870

B58 84483613348880041862046775994036021/354

=={{header|Java}}==

public class BernoulliNumbers {

return A[0];

}

<pre>B(0 ) = 1

B(1 ) = 1 / 2

'''BigInt Stubs''':

# def lessOrEqual(x; y): # x <= y

# def long_add(x;y): # x+y

def long_mod(x;y):

# A fraction is represented by [numerator, denominator] in reduced form, with the sign on top

| if $a == $b then ["0", "1"]

else add($a; [ ($b[0]|negate), $b[1] ] )

'''Bernoulli Numbers''':

def bernoulli(n):

reduce range(0; n+1) as $m

.[$j-1] = multiply( [($j|tostring), "1"]; minus( .[$j-1] ; .[$j]) ) ))

| .[0] # (which is Bn)

'''The task''':

{{out}}

The following output was obtained using the previously mentioned BigInt library.

0: ["1","1"]

1: ["1","2"]

56: ["-2479392929313226753685415739663229","870"]

58: ["84483613348880041862046775994036021","354"]

=={{header|Julia}}==

A = Vector{Rational{BigInt}}(undef, n + 1)

for m = 0 : n

for (n, b) in enumerate(BernoulliList(60))

isodd(numerator(b)) && println("B($(n-1)) = $b")

Produces virtually the same output as the Python version.

{{works with|Commons Math|3.3.5}}

object Bernoulli {

if (n % 2 == 0 || n == 1)

System.out.printf("B(%-2d) = %-1s%n", n, Bernoulli(n))

{{out}}

Produces virtually the same output as the Java version.

{{libheader|luagmp}}

{{works with|LuaJIT|2.0-2.1}}

local gmp = require 'gmp' ('libgmp')

local ffi = require'ffi'

gmp.z_clears(n,d)

gmp.q_clear(rop)

{{out}}

<pre>> time ./bernoulli_gmp.lua

=={{header|Maple}}==

{{out}}

<pre>[[0, 1], [1, 1/2], [2, 1/6], [4, -1/30], [6, 1/42], [8, -1/30], [10, 5/66], [12, -691/2730], [14, 7/6], [16, -3617/510], [18, 43867/798], [20, -174611/330], [22, 854513/138], [24, -236364091/2730], [26, 8553103/6], [28, -23749461029/870], [30, 8615841276005/14322], [32, -7709321041217/510], [34, 2577687858367/6], [36, -26315271553053477373/1919190], [38, 2929993913841559/6], [40, -261082718496449122051/13530], [42, 1520097643918070802691/1806], [44, -27833269579301024235023/690], [46, 596451111593912163277961/282], [48, -5609403368997817686249127547/46410], [50, 495057205241079648212477525/66], [52, -801165718135489957347924991853/1590], [54, 29149963634884862421418123812691/798], [56, -2479392929313226753685415739663229/870], [58, 84483613348880041862046775994036021/354], [60, -1215233140483755572040304994079820246041491/56786730]]</pre>

=={{header|Mathematica}} / {{header|Wolfram Language}}==

Mathematica has no native way for starting an array at index 0. I therefore had to build the array from 1 to n+1 instead of from 0 to n, adjusting the formula accordingly.

Do[

a[[m]] = 1/m;

, {m, 1, n + 1}];

]

{{out}}

<pre>{0,1}

(Note from task's author: nobody is forced to use any specific algorithm, the one shown is just a suggestion.)

{{out}}

<pre>0 1

=={{header|Nim}}==

import strformat

for (n, b) in values:

let s = fmt"{($b.num).alignString(maxLen, '>')} / {b.denom}"

{{out}}

=={{header|PARI/GP}}==

{{out}}

<pre>0 1

{{libheader|BigDecimalMath}}

Tested with fpc 3.0.4

(* Taken from the 'Ada 99' project, https://marquisdegeek.com/code_ada99 *)

end;

end.

{{out}}

Instead of doing the same calculations over and over again, I retain the A array until the final Bernoulli number is produced.

use strict;

use warnings;

bernoulli_print();

The output is exactly the same as the Python entry.

We can also use modules for faster results. E.g.

{{libheader|ntheory}}

for my $n (0 .. 60) {

my($num,$den) = bernfrac($n);

printf "B(%2d) = %44s/%s\n", $n, $num, $den if $num != 0;

with identical output. Or:

for my $n (0 .. 60) {

my($num,$den) = split "/", bernfrac($n);

printf("B(%2d) = %44s/%s\n", $n, $num, $den||1) if $num != 0;

with the difference being that Pari chooses <math>B_1</math> = -&frac12;.

{{libheader|Phix/mpfr}}

{{trans|C}}

<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>

<span style="color: #008080;">include</span> <span style="color: #000000;">builtins</span><span style="color: #0000FF;">/</span><span style="color: #004080;">mpfr</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span>

<span style="color: #0000FF;">{</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">d</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_free</span><span style="color: #0000FF;">({</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">d</span><span style="color: #0000FF;">})</span>

<span style="color: #000000;">rop</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpq_free</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rop</span><span style="color: #0000FF;">)</span>

{{out}}

<pre>

=={{header|PicoLisp}}==

Brute force and method by Srinivasa Ramanujan.

(de fact (N)

(test (berno N) (berno-brute N)) )

=={{header|PL/I}}==

declare i fixed binary;

declare B complex fixed (31);

(3 A, column(10), F(32), 2 A);

end;

The above uses GCD (see Rosetta Code) extended for 31-digit working.

=={{header|Python}}==

===Python: Using task algorithm===

def bernoulli(n):

width = max(len(str(b.numerator)) for i,b in bn)

for i,b in bn:

{{out}}

===Python: Optimised task algorithm===

Using the optimization mentioned in the Perl entry to reduce intermediate calculations we create and use the generator bernoulli2():

A, m = [], 0

while True:

width = max(len(str(b.numerator)) for i,b in bn2)

for i,b in bn2:

Output is exactly the same as before.

=={{header|Quackery}}==

[ 1+

char / over find

44 swap - times sp

{{out}}

{{incorrect|rsplus|This example is incorrect: It is not executable and if made executable (with 'library(gmp)') it returns completely different and wrong results -- not the ones shown here. The R code needs complete rewrite and the 'pracma' library will not be of any help.}}

library(pracma)

cat("B(",idx,") = ",n,"/",d,"\n", sep = "")

}

{{out}}

<pre>

use the same emmitter... it's just a matter of how long to wait for the emission.

;; For: http://rosettacode.org/wiki/Bernoulli_numbers

(list 1/1 (app abs 1/2) 1/6 -1/30 1/42 -1/30 _ ...))

; timing only ...

{{out}}

First, a straighforward implementation of the naïve algorithm in the task description.

{{works with|Rakudo|2015.12}}

my @a;

for 0..$n -> $m {

my $form = "B(%2d) = \%{$width}d/%d\n";

{{out}}

<pre>B( 0) = 1/1

{{works with|Rakudo|2015.12}}

my @a;

for 0..* -> $m {

my $form = "B(%d)\t= \%{$width}d/%d\n";

{{out}}

<pre>B(0) = 1/1

{{works with|Rakudo|2016.12}}

this.key => this.key * (this.value - prev.value)

}

my $form = "B(%d)\t= \%{$width}d/%d\n";

Same output as memoization example

<br>

:::::::::::: where &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; <big><big> <math> \binom kr</math> </big></big> &nbsp; &nbsp; &nbsp; is a binomial coefficient. <br>

parse arg N .; if N=='' | N=="," then N= 60 /*Not specified? Then use the default.*/

numeric digits max(9, n*2) /*increase the decimal digits if needed*/

/*──────────────────────────────────────────────────────────────────────────────────────*/

perm: procedure expose !.; parse arg x,y; if !.P.x.y\==. then return !.P.x.y

{{out|output|text=&nbsp; when using the default input:}}

<pre>

=={{header|Ruby}}==

{{trans|Python}}

ar = []

0.step do |m|

b_nums.each_with_index {|b,i| puts "B(%2i) = %*i/%i" % [i, width, b.numerator, b.denominator] unless b.zero? }

{{out}}

<pre>

=={{header|Rust}}==

// the optimized function. The speeds vary significantly: relative

// speeds of optimized:iterator:naive implementations is 625:25:1.

}

}

{{out}}

<pre>

'''With Custom Rational Number Class'''<br/>

(code will run in Scala REPL with a cut-and-paste without need for a third-party library)

case class BFraction( numerator:BigInt, denominator:BigInt ) {

require( denominator != BigInt(0), "Denominator cannot be zero" )

println( f"$label%-6s $num / ${b.den}" )

}}

{{out}}

<pre>b(0) 1 / 1

=={{header|Scheme}}==

{{works with|Chez Scheme}}

(define bernoulli

(else

(printf "B(~2@a) = ~a~%" index (rational-padded (car numbers) max-numerator-length))

{{out}}

<pre>$ scheme --script bernoulli.scm

function automatically writes repeating decimals in parentheses, when necessary.

include "bigrat.s7i";

end if;

end for;

{{out}}

=={{header|Sidef}}==

Built-in:

Recursive solution (with auto-memoization):

n.is_one && return 1/2

var Bn = bernoulli_number(n) || next

printf("B(%2d) = %44s / %s\n", n, Bn.nude)

Using Ramanujan's congruences (pretty fast):

return 1/2 if n.is_one

binomial(n+3, n - 6*k) * __FUNC__(n - 6*k)

})) / binomial(n+3, n)

Using Euler's product formula for the Riemann zeta function and the Von Staudt–Clausen theorem (very fast):

n.is_zero && return 1

(-1)**(n/2 + 1) * int(ceil(d*K / z)) / d

The Akiyama–Tanigawa algorithm:

var a = []

for m in (0..60) {

}

 

{{out}}

=={{header|SPAD}}==

{{works with|FriCAS, OpenAxiom, Axiom}}

for n in 0..60 | (b:=bernoulli(n)$INTHEORY; b~=0) repeat print [n,b]

Package:[http://fricas.github.io/api/IntegerNumberTheoryFunctions.html?highlight=bernoulli IntegerNumberTheoryFunctions]

Uses the Frac type defined in the [http://rosettacode.org/wiki/Arithmetic/Rational#Swift Rational] task.

public func bernoulli<T: BinaryInteger & SignedNumeric>(n: Int) -> Frac<T> {

print("B(\(n)) = \(b)")

{{out}}

=={{header|Tcl}}==

for {set m 0} {$m <= $n} {incr m} {

lappend A [list 1 [expr {$m + 1}]]

foreach {n num denom} $result {

puts [format {B_%-2d = %*lld/%lld} $n $len $num $denom]

{{out}}

<pre>

{{works with|Visual Basic .NET|2013}}

{{libheader|System.Numerics}}

Imports System.Numerics 'BigInteger

End Sub 'bernoulli_BigInt

{{out}}

<pre>

{{libheader|Wren-fmt}}

{{libheader|Wren-big}}

import "/big" for BigRat

var b = bernoulli.call(n)

if (b != BigRat.zero) Fmt.print("B($2d) = $44i / $i", n, b.num, b.den)

{{out}}

{{trans|EchoLisp}}

Uses lib GMP (GNU MP Bignum Library).

fcn init(_a,_b){ var a=_a, b=_b; normalize(); }

fcn toString{ "%50d / %d".fmt(a,b) }

}

fcn __opDiv(n){ self(a*n.b,b*n.a) } // Rat / Rat

fcn B(N){ // calculate Bernoulli(n)

var A=List.createLong(100,0); // aka static aka not thread safe

}

A[0]

{{out}}

<pre>