Документация Engee
Notebook

Вычисление чисел Бернулли

В этом примере мы реализуем и анализируем алгоритм нахождения чисел Бернулли на языке программирования Julia.

Введение

Числа Бернулли — это последовательность рациональных чисел, которые играют важную роль в теории чисел, математическом анализе и других разделах математики. Они появляются, например, в разложении функций в ряд Тейлора, в формуле Эйлера-Маклорена и при вычислении сумм степеней натуральных чисел. Алгоритм Акиямы–Танигавы позволяет эффективно вычислять числа Бернулли с помощью рекуррентных соотношений.

Основная часть

Реализация функции для вычисления отдельного числа Бернулли

Функция bernoulli(n) вычисляет n-е число Бернулли по алгоритму Акиямы–Танигавы. Параметр n - индекс числа Бернулли. Функция возвращает рациональное число типа Rational{BigInt}.

In [ ]:
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
Out[0]:
bernoulli (generic function with 1 method)

Вывод нечётных чисел Бернулли

Функция display(n) выводит на экран только нечётные числа Бернулли от до .
Нечётные числа Бернулли имеют нечётные числители.

In [ ]:
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
Out[0]:
display (generic function with 1 method)

Вызов функции для отображения чисел Бернулли до

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

Более эффективное вычисление списка чисел Бернулли

Функция BernoulliList(len) вычисляет сразу весь список чисел Бернулли от до , что более эффективно, чем вычислять их по одному.

In [ ]:
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
Out[0]:
BernoulliList (generic function with 1 method)

Перебираем список чисел и выводим только те, у которых числитель нечётный

In [ ]:
for (n, b) in enumerate(BernoulliList(60))
    # Проверка на нечётность числителя
    isodd(numerator(b)) && println("B($(n-1)) = $b")
end
B(0) = 1//1
B(1) = 1//2
B(2) = 1//6
B(4) = -1//30
B(6) = 1//42
B(8) = -1//30
B(10) = 5//66
B(12) = -691//2730
B(14) = 7//6
B(16) = -3617//510
B(18) = 43867//798
B(20) = -174611//330
B(22) = 854513//138
B(24) = -236364091//2730
B(26) = 8553103//6
B(28) = -23749461029//870
B(30) = 8615841276005//14322
B(32) = -7709321041217//510
B(34) = 2577687858367//6
B(36) = -26315271553053477373//1919190
B(38) = 2929993913841559//6
B(40) = -261082718496449122051//13530
B(42) = 1520097643918070802691//1806
B(44) = -27833269579301024235023//690
B(46) = 596451111593912163277961//282
B(48) = -5609403368997817686249127547//46410
B(50) = 495057205241079648212477525//66
B(52) = -801165718135489957347924991853//1590
B(54) = 29149963634884862421418123812691//798
B(56) = -2479392929313226753685415739663229//870
B(58) = 84483613348880041862046775994036021//354
B(60) = -1215233140483755572040304994079820246041491//56786730

Заключение

В этом примере мы изучили и реализовали алгоритм вычисления чисел Бернулли на языке Julia. Мы рассмотрели два подхода: вычисление по одному числу и вычисление всего списка сразу, второй вариант оказался более эффективным. Результат выводится в удобочитаемом формате, причем только числа с нечётными числителями, так как они представляют наибольший интерес. Это важный элемент при решении задач математического анализа, теории чисел и комбинаторики.

Пример разработан с использованием материалов Rosetta Code