Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Ruby script to compute Pi with arctan's formula.(old ver.)
#! /usr/local/bin/ruby
#*********************************************
# 円周率計算 by Arctan 系の公式
#
# 1: Machin
# 2: Klingenstierna
# 3: Euler
# 4: Euler(2)
#*********************************************
#
class CalcPiArctan
# 公式
FORMULA = ["Machin", "Klingenstierna", "Euler", "Euler2"] # 公式名
KK = [ # 公式内係数
[16, 1, 5, -4, 1, 239 ], # 1: Machin
[32, 1, 10, -4, 1, 239, -16, 1, 515], # 2: Klingenstierna
[20, 1, 7, 8, 3, 79 ], # 3: Euler
[16, 1, 5, -4, 1, 70, 4, 1, 99] # 4: Euler(2)
]
# 保存ファイル名
FILE_PRE = "PI_" # プリフィックス
FILE_EXT = ".txt" # 拡張子
# 出力文字列 ( コンソール、テキストファイル共通 )
STR_TITLE = "** Pi Computation by Arctan method **"
STR_FORMULA = " Formula = [ %s ]"
STR_DIGITS = " Digits = %d"
STR_TIME = " Time = %f seconds"
# 多桁計算
MAX_DIGITS = 100000000 # 1つの int で8桁扱う
def initialize(x, y)
# 結果出力ファイル名生成
@fname = "#{FILE_PRE}#{FORMULA[x - 1]}#{FILE_EXT}"
puts
puts "[ Output file : #{@fname} ]"
puts
@formula = FORMULA[x-1] # 公式名取得
@ks = KK[x-1].size / 3 # 計算対象公式の項数
@kk = KK[x-1] # 計算対象公式の係数
@l = y # 計算桁数
@l1 = (@l / 8) + 1 # 配列サイズ
n = @l / Math::log10(@kk[2]) + 1
@n = (n / 2).truncate + 1 # 計算項数
end
# 計算・結果出力
def calc
# ==== 計算開始時刻
t0 = Time.now
# ==== 配列宣言・初期化
s = Array.new(@l1 + 2, 0) # 総和
a = Array.new(@ks, nil) # 公式の各項
a.each_index {|i| a[i] = Array.new(@l1 + 2, 0)}
q = Array.new(@l1 + 2, 0) # 各項の和
# ==== 計算初期値
0.upto(@ks - 1) do |i|
a[i][0] = @kk[i * 3] * @kk[i * 3 + 2]
a[i][0] *= -1 if @kk[i * 3] < 0
# 分子が 1 で無ければ、さらに分子の値で除算しておく
a[i] = long_div(a[i], @kk[i * 3 + 1]) unless @kk[i * 3 + 1] == 1
end
# ==== 計算
1.upto(@n) do |k|
# 各項の分数計算 ( 1つ前の計算結果に追加で除算・乗算 )
0.upto(@ks - 1) do |i|
# 本来 x * x で除算したい(した方が高速だ)が、
# x が 8桁以下でも x * x が8桁超になるケースがあるため2回に分けている。
a[i] = long_div(a[i], @kk[i * 3 + 2])
a[i] = long_div(a[i], @kk[i * 3 + 2])
# 分子が 1 でない時
# 本来 x * x を乗算したい(した方が高速だ)が、
# x が 8桁以下でも x * x が8桁超になるケースがあるため2回に分けている。
unless @kk[i * 3 + 1] == 1
a[i] = long_mul(a[i], @kk[i * 3 + 1])
a[i] = long_mul(a[i], @kk[i * 3 + 1])
end
end
# 各項の加減算
1.upto(@ks - 1) do |i|
if @kk[i * 3] >= 0 # 加算
q = i == 1 ? long_add(a[i - 1], a[i]) : long_add(q, a[i])
else # 減算
q = i == 1 ? long_sub(a[i - 1], a[i]) : long_sub(q, a[i])
end
end
# 1 / (2 * k - 1) の部分
q = long_div(q, 2 * k - 1)
# 総和計算
s = (k % 2) == 0 ? long_sub(s, q) : long_add(s, q)
end
# ==== 計算終了時刻
t1 = Time.now
# ==== 計算時間
tt = t1 - t0
# ==== 結果出力
display(tt, s)
rescue => e
raise
end
# ロング + ロング
def long_add(a, b)
z = Array.new(@l1 + 2, 0)
cr = 0
(@l1 + 1).downto(0) do |i|
z[i] = a[i] + b[i] + cr
if z[i] < MAX_DIGITS
cr = 0
else
z[i] -= MAX_DIGITS
cr = 1
end
end
return z
rescue => e
raise
end
# ロング - ロング
def long_sub(a, b)
z = Array.new(@l1 + 2, 0)
br = 0
(@l1 + 1).downto(0) do |i|
z[i] = a[i] - b[i] - br
if z[i] >= 0
br = 0
else
z[i] += MAX_DIGITS
br = 1
end
end
return z
rescue => e
raise
end
# ロング * ショート
def long_mul(a, b)
z = Array.new(@l1 + 2, 0)
cr = 0
(@l1 + 1).downto(0) do |i|
w = a[i]
z[i] = (w * b + cr) % MAX_DIGITS
cr = (w * b + cr) / MAX_DIGITS
end
return z
rescue => e
raise
end
# ロング / ショート
def long_div(a, b)
z = Array.new(@l1 + 2, 0)
r = 0
0.upto(@l1 + 1) do |i|
w = a[i]
z[i] = (w + r) / b
r = ((w + r) % b) * MAX_DIGITS
end
return z
rescue => e
raise
end
# 結果出力
def display(tt, s)
puts STR_TITLE
puts STR_FORMULA % @formula
puts STR_DIGITS % @l
puts STR_TIME % tt
# ファイル出力
out_file = File.open(@fname, "w")
out_file.puts STR_TITLE
out_file.puts STR_FORMULA % @formula
out_file.puts STR_DIGITS % @l
out_file.puts STR_TIME % tt
out_file.puts "\n %d." % s[0]
1.upto(@l1 - 1) do |i|
out_file.printf("%08d:", (i - 1) * 8 + 1) if i % 10 == 1
out_file.printf(" %08d", s[i])
out_file.puts if i % 10 == 0
end
rescue => e
raise
end
end
if __FILE__ == $0
begin
# ==== 使用公式番号入力
print "1:Machin, 2:Klingenstierna, 3:Euler, 4:Euler(2) : "
f = gets.to_i
f = 1 if f < 1 || f > 8 # 範囲外なら 1:Machin と判断
# ==== 計算桁数入力
print "Please input number of Pi Decimal-Digits : "
n = gets.to_i
# 計算クラスインスタンス化
obj = CalcPiArctan.new(f, n)
# 円周率計算
obj.calc
rescue => e
$stderr.puts "[#{e.class}] #{e.message}\n"
e.backtrace.each{ |tr| $stderr.puts "\t#{tr}" }
end
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment