mk-mode BLOG

このブログは自作の自宅サーバに構築した Debian GNU/Linux で運用しています。
PC・サーバ構築等の話題を中心に公開しております。(クローンサイト: GitHub Pages

ブログ開設日2009-01-05
サーバ連続稼働時間
Reading...
Page View 合計
Reading...
今日
Reading...
昨日
Reading...

Ruby - 円周率計算(by マチンの公式)!

[ プログラミング, 数学 ] [ Ruby, 円周率 ]

こんばんは。

前回は、C++ による「マチンの公式による円周率計算」のアルゴリズムを紹介しました。

今日は、同じアルゴリズムを Ruby で実現してみました。
アルゴリズム等については、上記リンクの記事を参照してください。

実際、大体同じです。

以下、Ruby によるサンプルスクリプトです。

記録

0. 前提条件

  • Linux Mint 14 Nadia (64bit) での作業を想定。
  • Ruby 2.0.0-p0 を使用。

1. Ruby スクリプト作成

今回作成した Ruby ソースは以下の通りです。

calc_pi_machin.rb
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
#*********************************************
# 円周率計算 by マチンの公式
#*********************************************
#
L  = 1000                 # 算出桁数
L1 = ((L / 8) + 1)        # 配列サイズ
L2 = (L1 + 1)             # 配列サイズ + 1
N  = ((L / 1.39794) + 1)  # 計算項数

class CalcPiMachin
  def initialize
  end

  # 計算・結果出力
  def calc
    # 配列宣言・初期化
    s = Array.new(L2 + 1, 0)  # 総和
    a = Array.new(L2 + 1, 0)  # マチンの公式の前の項
    b = Array.new(L2 + 1, 0)  # マチンの公式の後の項
    q = Array.new(L2 + 1, 0)  # マチンの公式の前の項+後の項

    # マチンの公式
    a[0] = 16 * 5
    b[0] = 4 * 239
    1.upto(N) do |k|
      a = long_div(a, 5 * 5)
      b = long_div(b, 239 * 239)
      q = long_sub(a, b)
      q = long_div(q, 2 * k - 1)
      unless k % 2 == 0
        s = long_add(s, q)
      else
        s = long_sub(s, q)
      end
    end

    # 結果出力
    display(s)
  end

  # ロング + ロング
  def long_add(a, b)
    z = Array.new(N, 0)

    carry = 0
    L2.downto(0) do |i|
      z[i] = a[i] + b[i] + carry
      if z[i] < 100000000
        carry = 0
      else
        z[i] -= 100000000
        carry = 1
      end
    end

    return z
  end

  # ロング - ロング
  def long_sub(a, b)
    z = Array.new(N, 0)

    borrow = 0
    L2.downto(0) do |i|
      z[i] = a[i] - b[i] - borrow
      if z[i] >= 0
        borrow = 0
      else
        z[i] += 100000000
        borrow = 1
      end
    end

    return z
  end

  # ロング / ショート
  def long_div(a, b)
    z = Array.new(N, 0)

    r = 0
    0.upto(L2) do |i|
      w = a[i]
      z[i] = (w + r) / b
      r = ((w + r) % b) * 100000000
    end

    return z
  end

  # 結果出力
  def display(s)
    printf("%7d. ", s[0])
    1.upto(L1 - 1) do |i|
      printf("%08d ", s[i])
    end
    printf("\n")
  end
end

# メイン処理
begin
  # 計算クラスインスタンス化
  obj_calc = CalcPiMachin.new

  # 円周率計算
  obj_calc.calc
rescue => e
  # エラーメッセージ
  puts "[例外発生] #{e}"
end

ちなみに、Gist でもスクリプトを公開しているが、こちらのスクリプトは汎用性を高めたものとなっている。

2. 実行

実際に実行して検証してみる。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
$ ruby calc_pi_machin.rb
      3. 14159265 35897932 38462643 38327950 28841971 69399375 10582097
49445923 07816406 28620899 86280348 25342117 06798214 80865132 82306647
09384460 95505822 31725359 40812848 11174502 84102701 93852110 55596446
22948954 93038196 44288109 75665933 44612847 56482337 86783165 27120190
91456485 66923460 34861045 43266482 13393607 26024914 12737245 87006606
31558817 48815209 20962829 25409171 53643678 92590360 01133053 05488204
66521384 14695194 15116094 33057270 36575959 19530921 86117381 93261179
31051185 48074462 37996274 95673518 85752724 89122793 81830119 49129833
67336244 06566430 86021394 94639522 47371907 02179860 94370277 05392171
76293176 75238467 48184676 69405132 00056812 71452635 60827785 77134275
77896091 73637178 72146844 09012249 53430146 54958537 10507922 79689258
92354201 99561121 29021960 86403441 81598136 29774771 30996051 87072113
49999998 37297804 99510597 31732816 09631859 50244594 55346908 30264252
23082533 44685035 26193118 81710100 03137838 75288658 75332083 81420617
17766914 73035982 53490428 75546873 11595628 63882353 78759375 19577818
57780532 17122680 66130019 27876611 19590921 64201989

C++ 版と同じ結果が得られました。

3. 参考サイト


前回の C++ と今回の Ruby でのコーディングで、マチンの公式での円周率計算についての計算機科学的な考え方は理解できました。
(実際、もう少し綺麗にコーディングできるでしょうが)

以上。

Comments