■プロでは無いので、もっと簡単に計算する方法は無いものか。
以下の3つの公式を使う。これなら単純に、
組み込み関数でどこまで出来るかで充分。
atan 1 = π/4、acos(-1) = π、asin(1) = π/2
上記を変形して以下のようにします。acos(-1)はそのままです。
4 * acos = π、asin(1) * 2 = π
30万の精度を出そうとすると以下の通り、オーバーフローになる。
*** - Program stack overflow. RESET
どこが限界かはまたの機会に。。。
■アークタンジェントを使う方法
というか、こっちの方がBBPよりもコードが短く精度が高い。。。
BBPの方がはるかに優れたアルゴリズムですけれど。。。
3千桁までは直ぐに答えが出る。
※atan：Arctan tanの逆数 arctan 1 = π/4
$ echo "(* 4 (atan 1))" | clisp -q | tail -1
3.1415925
$ echo "(* 4.0L0 (atan 1.0L0))" | clisp -q | tail -1
3.141592653589793238L0
$ echo "(setf (EXT:LONG-FLOAT-DIGITS) 100) (* 4.00L0 (atan 1.00L0))" | clisp -q | tail -1
3.14159265358979323846264338327950288418L0
■小数点以下30桁
$ echo "(setf (EXT:LONG-FLOAT-DIGITS) 1000) (* 4.00L0 (atan 1.00L0))" | clisp -q | tail -1 | \
awk -F\. '{print $2}' | awk -F L '{print $1}' | \
sed s/"."/"&\n"/g | awk '{if (NR % 10 == 0) print $1 ":"; else print $1}' | \
xargs echo -n | sed s/" "//g | sed s/":"/"\n"/g | head -30 > atan.log
$ head -30 check.txt > atan.check; \
diff atan.log atan.check
■小数点以下300桁
$ echo "(setf (EXT:LONG-FLOAT-DIGITS) 10000) (* 4.00L0 (atan 1.00L0))" | \
clisp -q | tail -1 | \
awk -F\. '{print $2}' | awk -F L '{print $1}' | \
sed s/"."/"&\n"/g | awk '{if (NR % 10 == 0) print $1 ":"; else print $1}' | \
xargs echo -n | sed s/" "//g | sed s/":"/"\n"/g | head -300 >atan1.log
$ head -300 check.txt > atan1.check; \
diff atan1.log atan1.check
■小数点以下3000桁
$ echo "(setf (EXT:LONG-FLOAT-DIGITS) 100000) (* 00L0 (atan 1.00L0))" | \
clisp -q | tail -1 | \
awk -F\. '{print $2}' | awk -F L '{print $1}' | \
sed s/"."/"&\n"/g | awk '{if (NR % 10 == 0) print $1 ":"; else print $1}' | \
xargs echo -n | sed s/" "//g | sed s/":"/"\n"/g | head -3000 >atan2.log
$ head -3000 check.txt > atan2.check; \
diff atan2.log atan2.check
■アークコサインを使う方法
アークタンジェントよりも3万桁以降が速いが、
直ぐに答えを出したいならば、3千桁までだろう。
cos(π)=-1の逆関数、acos(-1)=π
■小数点以下30桁
$ echo "(setf (EXT:LONG-FLOAT-DIGITS) 1000)(acos -1.0L0)" | \
clisp -q | tail -1 | \
awk -F\. '{print $2}' | awk -F L '{print $1}' | \
sed s/"."/"&\n"/g | awk '{if (NR % 10 == 0) print $1 ":"; else print $1}' | \
xargs echo -n | sed s/" "//g | sed s/":"/"\n"/g | head -30 > acos.log
$ diff acos.log atan.check
■小数点以下300桁
$ echo "(setf (EXT:LONG-FLOAT-DIGITS) 10000)(acos -1.0L0)" | \
clisp -q | tail -1 | \
awk -F\. '{print $2}' | awk -F L '{print $1}' | \
sed s/"."/"&\n"/g | awk '{if (NR % 10 == 0) print $1 ":"; else print $1}' | \
xargs echo -n | sed s/" "//g | sed s/":"/"\n"/g | head -300 > acos1.log
$ diff acos1.log atan1.check
■小数点以下3000桁
$ echo "(setf (EXT:LONG-FLOAT-DIGITS) 100000)(acos -1.0L0)" | \
clisp -q | tail -1 | \
awk -F\. '{print $2}' | awk -F L '{print $1}' | \
sed s/"."/"&\n"/g | awk '{if (NR % 10 == 0) print $1 ":"; else print $1}' | \
xargs echo -n | sed s/" "//g | sed s/":"/"\n"/g | head -3000 > acos2.log
$ diff acos2.log atan2.check
■小数点以下3万桁
$ echo "(setf (EXT:LONG-FLOAT-DIGITS) 1000000)(acos -1.0L0)" | \
clisp -q | tail -1 | \
awk -F\. '{print $2}' | awk -F L '{print $1}' | \
sed s/"."/"&\n"/g | awk '{if (NR % 10 == 0) print $1 ":"; else print $1}' | \
xargs echo -n | sed s/" "//g | sed s/":"/"\n"/g | head -30000 > acos3.log
$ diff acos3.log atan3.check
■アークサインを使う方法。
こちらも3千桁までは直ぐに結果が出る。
asin(-1)=-π/2、asin(1)=π/2
$ echo "(* 2 (asin 1))" | clisp -q | tail -1
3.1415927
$ echo "(* 2 (asin 1.0L0))" | clisp -q | tail -1
3.1415926535897932385L0
$ echo "(setf (EXT:LONG-FLOAT-DIGITS) 100000)(* 2 (asin 1.0L0))" | \
clisp -q | tail -1 | \
awk -F\. '{print $2}' | awk -F L '{print $1}' | \
sed s/"."/"&\n"/g | awk '{if (NR % 10 == 0) print $1 ":"; else print $1}' | \
xargs echo -n | sed s/" "//g | sed s/":"/"\n"/g | head -3000 > asin2.log
$ diff asin2.log atan2.check
□おまけ---
■debianパッケージのpiコマンドの中身
$ dpkg -L pi | grep 'bin/'
/usr/bin/pi
$ sudo apt-get install dpkg-dev
$ sudo apt-get source pi
■「cln-1.2.2」
$ grep = cln-1.2.2/examples/pi.cc
long digits = 100;
if (argc == 2 && !strcmp(argv[1],"--help")) {
if (argc == 2 && !strcmp(argv[1],"--version")) {
if (argc == 2 && !strcmp(argv[1],"--bibliography")) {
if (argc == 2 && isdigit(argv[1][0])) {
digits = atol(argv[1]);
cl_F p = pi(float_format(digits));
cpf.default_float_format = float_format(p);
istreambuf_iterator<char> i = buf.rdbuf();
■１行目：値を指定しなかった場合は、小数点以下100桁返します。
これはマニュアルにもあります。
$ pi | awk -F\. '{print $2}' | wc -c
100
$ LANG=C man pi | grep -A 2 100
The pi command prints 100 decimal digits of Archimedes' constant pi or
a number of digits specified by an integer parameter on the command
line.
■２〜４行目：以下の通り説明を返します。
$ LANG=C man pi | grep "\-\-"
--help Show summary of options.
--version
--bibliography
■５〜６行目：数値に変換します。
「atol」は、「L0」等の文字が含まれていた場合の対処でしょうか。
■７行目がメインです。
$ find cln-1.2.2/src/float/transcendental | \
for list in `xargs`;do \
grep cl_F "$list" > /dev/null 2>&1 && echo "$list"; \
done | sed s%"cln-1.2.2/src/float/transcendental/"%%g | \
grep pi | column
cl_F_pi_f.cc cl_F_pi.cc cl_F_roundpi2.cc
cl_F_pi_def.cc cl_F_pi_var.cc
cl_F_roundpi.cc cl_LF_pi.cc
■以下のように振り分け
$ grep return cln-1.2.2/src/float/transcendental/cl_F_pi.cc
, return cl_SF_pi;
, return cl_FF_pi;
, return cl_DF_pi;
, return pi(TheLfloat(y)->len);
■ヘッダファイルを探します。
$ find cln-1.2.2/src/float/transcendental/ | \
for list in `xargs`;do \
grep cl_SF_pi "$list" > /dev/null 2>&1 && echo "$list"; \
done | sed s%"cln-1.2.2/src/float/transcendental/"%%g | grep .h
cl_F_tran.h
■「cl_round_pi」と「cl_round_pi2」が本体のようです。
$ grep ^extern cln-1.2.2/src/float/transcendental/cl_F_tran.h | grep pi
extern const cl_SF cl_SF_pi;
extern const cl_FF cl_FF_pi;
extern const cl_DF cl_DF_pi;
extern cl_LF cl_LF_pi; // as long as it has ever been computed
extern const cl_LF pi (uintC len); // computes it even further
extern const cl_F_div_t cl_round_pi (const cl_F& x);
extern const cl_F_div_t cl_round_pi2 (const cl_F& x);
■アルキメデスのアルゴリズムです。
$ grep -A 8 const cln-1.2.2/src/float/transcendental/cl_F_roundpi.cc
const cl_F_div_t cl_round_pi (const cl_F& x)
{
if (float_exponent(x) <= 0)
// Exponent <=0 -> |x|<1 -> |x/pi| < 1/2, also Division unnötig
return cl_F_div_t(0,x); // Quotient 0, Rest x
else
// x durch pi (mit hinreichender Genauigkeit) dividieren
return round2(x,pi(x));
}
$ grep -A 8 ^const cln-1.2.2/src/float/transcendental/cl_F_roundpi2.cc
const cl_F_div_t cl_round_pi2 (const cl_F& x)
{
if (float_exponent(x) < 0)
// Exponent <0 -> |x|<1/2 -> |x/(pi/2)| < 1/2, also Division unnötig
return cl_F_div_t(0,x); // Quotient 0, Rest x
else
// x durch pi/2 (mit hinreichender Genauigkeit) dividieren
return round2(x,scale_float(pi(x),-1));
}
■アルキメデスの手順は歴史が長いので、色々解りやすい資料があります。
参考：円周率を求めて
http://www.chikyo.co.jp/maths/pdf/free02.pdf
■ところで、「cl_LF_pi.cc」が気になります。
ラマヌジャンのアルゴリズムですか。
$ grep -v "\/\/" cln-1.2.2/src/float/transcendental/cl_LF_pi.cc | grep ^const
const cl_LF compute_pi_brent_salamin (uintC len)
const cl_LF compute_pi_brent_salamin_quartic (uintC len)
const cl_LF compute_pi_ramanujan_163 (uintC len)
const cl_LF compute_pi_ramanujan_163_fast (uintC len)
const cl_LF pi (uintC len)
参考：cln 1.3.2
http://fossies.org/dox/cln-1.3.2/cl__LF__pi_8cc.html
参考：cl_LF_pi.cc
http://www.koders.com/cpp/fid738153EC24692DE267F68A7BC6E2B54799243682.aspx?s=mdef%3Acompute
■その他の方法
参考：計算プログラム紹介
http://xn--w6q13e505b.jp/program/introduction.html
参考：第8回勉強会
http://nseg.jp/?%E7%AC%AC8%E5%9B%9E%E5%8B%89%E5%BC%B7%E4%BC%9A
