Engee 文档
Notebook

计算伯努利数

在这个例子中,我们实现并分析了一种用Julia编程语言查找伯努利数的算法。

导言

伯努利数是一系列有理数,在数论,数学分析和数学的其他分支中起着重要作用。 例如,它们出现在将函数分解为泰勒级数、欧拉-麦克劳林公式以及计算自然数的幂和中。 Akiyama-Tanigawa算法使得使用复发关系有效地计算伯努利数成为可能。

主要部分

计算单个伯努利数的函数的实现

功能 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代码]的材料开发的(https://rosettacode.org/wiki/Bernoulli_numbers