Calculating Bernoulli numbers
In this example, we implement and analyze an algorithm for finding Bernoulli numbers in the Julia programming language.
Introduction
Bernoulli numbers are a sequence of rational numbers that play an important role in number theory, mathematical analysis, and other branches of mathematics. They appear, for example, in the decomposition of functions into a Taylor series, in the Euler-Maclaurin formula, and in calculating sums of powers of natural numbers. The Akiyama–Tanigawa algorithm makes it possible to efficiently calculate Bernoulli numbers using recurrence relations.
The main part
Implementation of a function for calculating a single Bernoulli number
Function bernoulli(n) calculates the nth Bernoulli number using the Akiyama–Tanigawa algorithm. Parameter n is the index of the Bernoulli number. The function returns a rational number of the type Rational{BigInt}.
function bernoulli(n)
# Creating a vector for intermediate calculations of size n+1, elements of type Rational{BigInt}
A = Vector{Rational{BigInt}}(undef, n + 1)
# We run the cycle on m from 0 to n
for m = 0 : n
# Initialize the first element according to the formula: A[m+1] = 1 / (m+1)
A[m + 1] = 1 // (m + 1)
# Reverse cycle on j from m to 1 to update the elements of A
for j = m : -1 : 1
# We update A[j] using the recurrent formula
A[j] = j * (A[j] - A[j + 1])
end
end
# We return the first value of vector A, which corresponds to the Bernoulli number B_n
return A[1]
end
Derivation of odd Bernoulli numbers
Function display(n) displays only odd Bernoulli numbers from before .
Odd Bernoulli numbers have odd numerators.
function display(n)
# We calculate all Bernoulli numbers from 0 to n
B = map(bernoulli, 0 : n)
# Defining the width for output alignment: the maximum number of characters in the numerator
pad = mapreduce(x -> ndigits(numerator(x)) + Int(x < 0), max, B)
# Calculating the width of the header (for example, the number B(60) will be 2 characters)
argdigits = ndigits(n)
# Iterating through the indexes from 0 to n
for i = 0 : n
# We check whether the numerator of the Bernoulli number is odd
if numerator(B[i + 1]) & 1 == 1
# We output the result as "B(0) = 1/1"
println(
"B(", lpad(i, argdigits), ") = ",
lpad(numerator(B[i + 1]), pad), " / ", denominator(B[i + 1])
)
end
end
end
Calling a function to display Bernoulli numbers up to
display(60)
More efficient calculation of a list of Bernoulli numbers
Function BernoulliList(len) calculates the entire list of Bernoulli numbers from before which is more efficient than calculating them one at a time.
function BernoulliList(len)
# Vector for intermediate values
A = Vector{Rational{BigInt}}(undef, len + 1)
# Vector for storing the result
B = similar(A)
# Cycle by index n from 0 to len
for n in 0 : len
# Initialize A[n+1] = 1/(n+1)
A[n + 1] = 1 // (n + 1)
# The inner loop updates A[j]
for j = n : -1 : 1
A[j] = j * (A[j] - A[j + 1])
end
# Storing the nth Bernoulli number (the first element of A) in B
B[n + 1] = A[1]
end
# We return the vector with Bernoulli numbers
return B
end
We go through the list of numbers and output only those with an odd numerator.
for (n, b) in enumerate(BernoulliList(60))
# Checking for the odd numerator
isodd(numerator(b)) && println("B($(n-1)) = $b")
end
Conclusion
In this example, we have studied and implemented an algorithm for calculating Bernoulli numbers in the Julia language. We considered two approaches: calculating one number at a time and calculating the entire list at once, the second option turned out to be more effective. The result is displayed in a human-readable format, with only numbers with odd numerators, as they are of the greatest interest. This is an important element in solving problems of mathematical analysis, number theory and combinatorics.
The example was developed using materials from Rosetta Code