计算伯努利数
在这个例子中，我们实现并分析了一种用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]:
奇伯努利数的推导
功能 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]:
调用一个函数来显示伯努利数字
In [ ]:
display(60)
更有效地计算伯努利数字列表
功能 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
# 在B中存储第n个伯努利数（A的第一个元素）
B[n + 1] = A[1]
end
# 我们用伯努利数返回向量
return B
end
Out[0]:
我们通过数字列表并仅输出具有奇数分子的数字。
In [ ]:
for (n, b) in enumerate(BernoulliList(60))
# 检查奇数分子
isodd(numerator(b)) && println("B($(n-1)) = $b")
end
结论
在这个例子中，我们研究并实现了一种用Julia语言计算伯努利数的算法。 我们考虑了两种方法：一次计算一个数字，一次计算整个列表，第二种选择结果更有效。 结果以人类可读的格式显示，只有具有奇数编号的数字，因为它们是最感兴趣的。 这是解决数学分析，数论和组合学问题的重要元素。
该示例是使用[罗塞塔代码]的材料开发的（https://rosettacode.org/wiki/Bernoulli_numbers ）