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)
# Создаём вектор для промежуточных вычислений размером n+1, элементы типа Rational{BigInt}
A = Vector{Rational{BigInt}}(undef, n + 1)
# Запускаем цикл по m от 0 до n
for m = 0 : n
# Инициализируем первый элемент согласно формуле: A[m+1] = 1 / (m+1)
A[m + 1] = 1 // (m + 1)
# Обратный цикл по j от m до 1 для обновления элементов A
for j = m : -1 : 1
# Обновляем A[j] по рекуррентной формуле
A[j] = j * (A[j] - A[j + 1])
end
end
# Возвращаем первое значение вектора A, которое соответствует числу Бернулли 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)
# Вычисляем все числа Бернулли от 0 до n
B = map(bernoulli, 0 : n)
# Определяем ширину для выравнивания вывода: максимальное количество символов в числителе
pad = mapreduce(x -> ndigits(numerator(x)) + Int(x < 0), max, B)
# Вычисляем ширину заголовка (например, номера B(60) будет 2 символа)
argdigits = ndigits(n)
# Итерация по индексам от 0 до n
for i = 0 : n
# Проверяем, является ли числитель числа Бернулли нечётным
if numerator(B[i + 1]) & 1 == 1
# Выводим результат в виде "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)
# Вектор для промежуточных значений
A = Vector{Rational{BigInt}}(undef, len + 1)
# Вектор для хранения результата
B = similar(A)
# Цикл по индексу n от 0 до len
for n in 0 : len
# Инициализируем A[n+1] = 1/(n+1)
A[n + 1] = 1 // (n + 1)
# Внутренний цикл обновляет A[j]
for j = n : -1 : 1
A[j] = j * (A[j] - A[j + 1])
end
# Сохраняем n-е число Бернулли (первый элемент A) в B
B[n + 1] = A[1]
end
# Возвращаем вектор с числами Бернулли
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))
# Проверка на нечётность числителя
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