#!/usr/bin/perl -w
# http://www.cyut.edu.tw/~ckhung/
# 從命令列上讀入 prec, 計算 pi=3.1416 ... 精準至小數以下 prec 位.

# 選項:
# -p prec	指定精準度
# -v		verbose, 過程中多印一些訊息

# 演算法: 利用 pi ~= 2^k * sin (pi / 2^k) 所導出的遞迴關係式:
#     pi(k+1) = pi(k) / cos(pi / 2^(k+1))
#     cos(th/2) = sqrt((1 + cos(th) / 2)
# 來計算 pi. (想像單位圓內接正 2^k * 6 邊形)
# 每次運算約增加 2 個 bit 的精準度, 故迴圈要執行 prec * 10 / 3 / 2 次

use Getopt::Std;
use Math::BigFloat;

my (
    %opts,		# 命令列選項
    $prec,		# 精確到小數下幾位?
    $p,			# pi 的近似值
    $c,			# 正多邊形的 "一頁" (很窄) 的內角 cos 值
    $i,
);

getopts('vp:', \%opts);

$prec = $opts{p} || 60;
Math::BigFloat::div_scale($prec);
# 也可以直接去修改一個全域變數: (不建議這麼做)
# $Math::BigFloat::div_scale = $prec;

$p = Math::BigFloat->new("3");
$c = Math::BigFloat->new("0.5");
$c = (($c + 1) / 2)->bsqrt();
## 上句其實是簡寫. 原來比較麻煩的寫法如下, 反而還不能執行):
# $c = Math::BigFloat::bsqrt(
#    Math::BigFloat::fdiv( Math::BigFloat::fadd($c, 1), 2)
# );
for ($i=1; $i<=$prec/3*5; ++$i) {
    printf "%3d: %s\n", $i, substr($p, 0, $prec+2) if $opts{v};
    $c = (($c + 1) / 2)->bsqrt();
#    $c = Math::BigFloat::bsqrt(($c + 1) / 2);	# 一樣的意思, 但無 inheritance
    $p = $p / $c;
}
printf "%s\n", substr($p, 0, $prec+2);