Julia 中关于 Π 的一些趣味

图片由 Cormullion 提供，代码 在此。
这篇文章以 Jupyter notebook 的形式提供 在此。

Julia 中的 π

作者：Simon Byrne

像大多数技术语言一样，Julia 为 π 提供了一个变量常量。但是 Julia 的处理方式有点特殊。

pi

π = 3.1415926535897...

它也可以通过 Unicode 符号访问（您可以在 REPL 或笔记本中通过 TeX 自动完成功能获得它，方法是输入 \pi 然后按 Tab 键）。

π

π = 3.1415926535897...

您会注意到它打印的格式与普通浮点数不同：这是因为它不是浮点数。

typeof(pi)

Irrational{:π}

π 和其他一些无理常数以特殊的 Irrational 值存储，而不是四舍五入为 Float64。它们的行为类似于普通的数值，除了它们可以自动转换为任何浮点类型，而无需任何中间舍入。

1 + pi # integers are promoted to Float64 by default

4.141592653589793

Float32(1) + pi # Float32

4.141593f0

这对于与任意精度 BigFloat 一起使用特别有用，因为 π 可以计算到完整的精度（而不是被截断为 Float64 然后转换回来）。

BigFloat(1) + pi # 256 bits by default

4.141592653589793238462643383279502884197169399375105820974944592307816406286198

如果 π 存储为 Float64，我们将得到

BigFloat(1) + Float64(pi)

4.141592653589793115997963468544185161590576171875000000000000000000000000000000

实际上，BigFloat（使用 MPFR 库）将根据需要计算到当前精度下的 π，该精度通过 setprecision 设置。这提供了一种简单的方法来获取它的数字。

# to 1024 bits
setprecision(BigFloat, 1024) do
    BigFloat(pi)
end

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724586997

最后几个打印的数字可能不正确，因为从 BigFloat 的内部二进制格式转换为用于打印的十进制表示形式。但这只是一个显示问题，内部二进制表示是正确舍入到最后一位的。

Irrational 的另一个巧妙特性是不等式是正确的。

Float64(pi) < pi < nextfloat(Float64(pi))

true

通过内联汇编指令获取 π

作者：Simon Byrne

Julia 提供了一个非常底层的 llvmcall 接口，允许用户直接编写 LLVM 中间表示，包括使用内联汇编。以下代码片段调用了 fldpi 指令（“floating point load pi”），该指令将常量 π 加载到浮点寄存器堆栈上（这仅适用于 x86 和 x86_64 架构）。

function asm_pi()
    Base.llvmcall(
    """ %pi = call double asm "fldpi", "={st}"()
        ret double %pi""",
    Float64, Tuple{})
end

asm_pi (generic function with 1 method)

asm_pi()

3.141592653589793

我们可以查看生成的实际代码。

@code_native asm_pi()

    	.section	__TEXT,__text,regular,pure_instructions
    Filename: In[10]
    	pushq	%rbp
    	movq	%rsp, %rbp
    Source line: 2
    	fldpi
    	fstpl	-8(%rbp)
    	movsd	-8(%rbp), %xmm0         ## xmm0 = mem[0],zero
    	popq	%rbp
    	retq

如果您想知道这些指令的其他部分在做什么

1. pushq 和 movq 添加到 调用堆栈帧 中。

2. fldpi 将 π 推送到 x87 浮点寄存器堆栈。

• x87 是较旧的传统浮点指令集，可以追溯到原始的 Intel 8087 协处理器。

3. fstpl 和 movsd 将值移动到 SSE 浮点寄存器 xmm0。

• Julia 与大多数现代软件一样，使用更新的 SSE 指令集进行浮点运算。这也使我们能够利用诸如 SIMD 操作 之类的东西。

4. popq 和 retq 弹出调用堆栈帧。

使用泰勒级数展开式计算 π

作者：Luis Benet，墨西哥国立自治大学（UNAM）物理科学研究所。

这将演示如何使用 TaylorSeries.jl 包通过各种泰勒级数展开式来计算 π。

using TaylorSeries

马德哈瓦公式

标准三角恒等式之一是

$\tan\left( \frac{\pi}{6} \right) = \frac{1}{\sqrt{3}}.$

因此，通过取 $6 \arctan(x)$ 在 0 处的泰勒展开式，我们可以通过在 $1/\sqrt{3}$处对其进行评估来获得 $\pi$ 的值，该值在收敛半径内。

我们使用 BigFloat 获取 37 阶泰勒级数。

series1 = 6atan( Taylor1(BigFloat, 37) )
convert(Taylor1{Rational{BigInt}},series1)

6//1 t - 2//1 t³ + 6//5 t⁵ - 6//7 t⁷ + 2//3 t⁹ - 6//11 t¹¹ + 6//13 t¹³ - 2//5 t¹⁵ + 6//17 t¹⁷ - 6//19 t¹⁹ + 2//7 t²¹ - 6//23 t²³ + 6//25 t²⁵ - 2//9 t²⁷ + 6//29 t²⁹ - 6//31 t³¹ + 2//11 t³³ - 6//35 t³⁵ + 6//37 t³⁷ + 𝒪(t³⁸)

请注意，上述级数仅包含奇数幂，因此在这种情况下我们将使用 18 个系数。

在 $1/\sqrt{3}$处评估该表达式，我们得到

pi_approx1 = evaluate(series1, 1/sqrt(big(3)))

3.141592653647826046431202390582141253830948237428790668441592864548346569098516

然后，37 阶泰勒展开式得到的值与 $\pi$ 的差异为

abs(pi - pi_approx1)

5.803280796855900730263836963377883805368484746664827224053016281231814650118929e-11

为了获得更精确的结果，我们可以简单地增加展开的阶数。

series2 = 6atan( Taylor1(BigFloat,99) ) # 49 coefficients of the series
pi_approx2 = evaluate(series2, 1/sqrt(BigInt(3)))

3.141592653589793238462643347272152237127662423839333289949470742535834074912581

abs(pi - pi_approx2)

3.600735064706950697553577253102547384977198233137361734413175534929622111373249e-26

此公式是 马德哈瓦 或 格雷戈里-莱布尼茨级数 之一。

$\pi = 6 \sum_{n=0}^{\infty} (-1)^n \frac{(1/\sqrt{3})^{2n+1}}{2n+1}.$