使用 Viete 公式计算 Pi 近似值时卡在 ocaml 的无限循环中

Question

我在计算给定精度的 Pi 近似值时遇到问题。我已经得出结论，无限循环问题是由我的循环退出条件引起的，但是我不知道确切的问题是什么。在我看来退出条件应该是这样的

abs(current_aproximation - previous_approximation) < precision

代码如下：

let pi(prec) = 
let rec loop(curr, prev) = 
    if(abs_float( (2. /. curr) -. (2. /. prev) ) < prec) then  // problematic line
        (2. /. curr)
    else
        loop(curr *. sqrt(2. +. curr) /. 2., curr)
    in loop(sqrt(0.5), 1.);;

感谢您提供解决问题的任何提示。

Answer 1

如果您修改代码以打印出 curr 的值，您会发现您很快到达了一个固定点：

9.88131291682e-324
9.88131291682e-324
9.88131291682e-324
9.88131291682e-324
9.88131291682e-324
9.88131291682e-324
9.88131291682e-324
9.88131291682e-324
9.88131291682e-324
(...)

根据 Wikipedia 你的初始值应该是 (sqrt 2.) /. 2. 而不是 sqrt (0.5)。但即使进行了此修改，也不可能在不达到前面描述的固定点的情况下请求小于 0.61 的精度。

我的猜测是 float 不够精确，无法以这种方式表达该算法。

Answer 2

您正在计算不同的产品：在一次迭代后，您应该计算 sqrt(2)/2 * sqrt(2+sqrt(2))/2 但您计算的是 sqrt(2)/2 * sqrt( 2+平方(2)/2)/2.

这个算法怎么样？

let pi prec =
   let rec p2 xn root =
     let nroot = sqrt(2. +. root) in
     let xm = xn *. (nroot /. 2.) in
     if (abs_float (( 2. /. xm ) -. (2. /. xn))) < prec
     then xm
     else p2 xm nroot
   in 2. /. (p2 1.0 0.0)

它的计算结果为：

# pi 0.1 ;;
- : float = 3.12144515225805197
# pi 0.01 ;;
- : float = 3.14033115695475251
# pi 0.001 ;;
- : float = 3.14127725093277288
# pi 0.0001 ;;
- : float = 3.14157294036709134

使用 Viete 公式计算 Pi 近似值时卡在 ocaml 的无限循环中

Stuck in infinite loop in ocaml while calculating Pi approximation using Viete's formula

ocaml

pi

approximation