Rubyでラマヌジャンに挑戦

Ruby メモ

先日行われたLL行事http://ll.jus.or.jp/2009のライトニングトーク、「HP50gでラマヌジャンに挑戦 by 大野典宏」という演目がありました。

大野さんの日記にも出演予告(2009-08-26 - 大野の日常)があります。
ラマヌジャンの公式にしたがって円周率πの数値計算をしてみるという話し、演目では Maximaで確かめながらヒューレット・パッカードの電卓HP50g(http://h20426.www2.hp.com/product/calculators/jp/ja/graphing/50g/index.aspx)で挑戦し、大変面白く結果が出ていた模様。
で、まあ、折角なので Ruby で計算してみるわけです。級数展開はスライド(PDF)3枚目に出てくる奴、上の方、それを割算して計算結果としてπが得られるようにしたもの。

何も考えずに Float で

1.9のRuby(ruby 1.9.1p243 (2009-07-16 revision 24175) [i386-mswin32])の irb。取り敢えず階乗の計算を定義する。

irb(main):001:0> class Integer; def fact; self>0 ?(1..self).inject(:*) : 1; end
; end
=> nil

何気に 1.9 の #inject 記法良いものです。

irb(main):004:0> y = 1 / ((0..41).inject(0){|r,n| r + (4*n).fact * (1103 + 26390*n) / ((4.0**n) * (99**n) * n.fact)**4} * 2 * Math::sqrt(2) / (99 ** 2))
=> 3.14159265358979

計算出来てるというか有効桁数なに、こういう計算で Float は無いよね。ちなみに項数41にしてるんだけど、42項で NaN。

irb(main):005:0> y = 1 / ((0..42).inject(0){|r,n| r + (4*n).fact * (1103 + 26390*n) / ((4.0**n) * (99**n) * n.fact)**4} * 2 * Math::sqrt(2) / (99 ** 2))
=> NaN

BigDecimal で何も考えない

というわけで BigDecimalライブラリを使います。

irb(main):011:0> require 'bigdecimal'
=> true
irb(main):012:0> require 'bigdecimal/util'
=> true

util を require しておくと #to_d メソッドが使える。

irb(main):017:0> y = (99.0**2).to_d / (0..200).map{|n| ( (4*n).fact * (1103+26390*n) ).to_s.to_d / (( (4**n) * (99**n) * n.fact )**4).to_s.to_d }.inject(:+) / (2.0.to_d.sqrt(1000)*2)
=> #<BigDecimal:b7e204,'0.3141592653 5897932384 6264338327 9502855789 3248380888
＜云々＞

ずっと続きますが、一行目35桁くらいしか正しくないです、「 95028」まで。

  • 取り敢えず 200項ばかり足してみた、ここでも「#inject(:+)」としてみる
  • 割とやみくもに「#to_d」してます
  • 「#to_s.to_d」注意、大きい整数計算は BigNum になるのだけど Integer型は直接「#to_d」出来ない、「#to_f」経由にすると Floatオーバーフローするので「#to_s」経由にします
  • 前は「Math::sqrt」使ってたけど、BigDecimal には「#sqrt(精度)」メソッドがあるのでそれを使う、今回 1000桁にしたのは余り意味ないです、1000なら一瞬、10000だと一拍待つも5秒は掛からない感じ(coresolo 1.2GHz/1.5G)

BigDecimal でちょっと考える

35桁しか精度が出ないのはなんなのでもう少し考えてみます。
その前に、階乗の計算、1.9ならメソッド名に「!」が使えるのでそうしましょう、なんかちょっとそれっぽくなります「.」一つはさまるけど、そしてその計算結果は Integerのクラス変数に残す様にします(参考: フラナガンasin:4873113946 P.363 ハッシュで階乗) 気休めですね。

irb(main):023:0> class Integer; @@factorial=Hash.new{|h,k| h[k] = k>1 ? k*h[k-1] : 1}; def !; @@factorial[self]; end; end
=> nil
irb(main):024:0> class Integer; def self.factorial; @@factorial; end; end
=> nil

「Integer.factorial」はなんというか、まあ、確認用かな。
それとπの桁数確認、目の子で見るのは大変です。π計算結果 1.6万桁辺りを参考に確認メソッドを作っておきます。

module BD
  PI = <<EOp.gsub(/\s/, '')
3.

1415926535 8979323846 2643383279 5028841971 6939937510
＜云々＞
EOp

  def PI.compare(y)
    y_to_s = y.to_s
    pi = '0.3'
    each_char.with_index do |v, i|
      next if i<2; (puts "#{pi}:#{i-1}"; break) unless v == y_to_s[i+1]; pi << v;
    end # each_char.with_index do |v, i|
  end # def PI.compare(y)
end

モジュールを作って「BD::PI.compare(y)」としました、y がレシーバになる方が分かり易いかもしれないけどそこまではこだわらない。最初 2項飛ばしてるのは「3.14」とか「0.314」とかの表記ずれ対策。
で、

irb(main):003:0> t=Time.now; BD::PI.compare (99.0**2).to_d / (0..200).map{|n| ((4*n).! * (1103+26390*n) ).to_s.to_d.div( (( (4**n) * (99**n) * n.! )**4).to_s.to_d, 1000 ) }.inject(:+) / (2.0.to_d.sqrt(1000)*2); puts Time.now-t
0.314159265358979323846264338327950288419716939937510582097494459230781640628620
89986280348253421170679821480865132823066470938446095505822317253594081284811174
50284102701938521105559644622948954930381964428810975665933446128475648233786783
16527120190914564856692346034861045432664821339360726024914127372458700660631558
81748815209209628292540917153643678925903600113305305488204665213841469519415116
09433057270365759591953092186117381932611793105118548074462379962749567351885752
72489122793818301194912983367336244065664308602139494639522473719070217986094370
27705392171762931767523846748184676694051320005681271452635608277857713427577896
09173637178721468440901224953430146549585371050792279689258923542019956112129021
96086403441815981362977477130996051870721134999999837297804995105973173281609631
85950244594553469083026425223082533446850352619311881710100031378387528865875332
08381420617177669147303598253490428755468731159562863882353787593751957781857780
5321712268066130019278766111959092164201989380952:1007
3.6875
=> nil
  • 200項で 1000桁、4秒。まあ良いところじゃないでしょうか
  • 級数各項の割算をちゃんと桁数指定するのがポイントでした
    • √2 も一応それに合わせる

まあでも手軽にやるにはこの辺が限界かな。項目数と桁数指定のバランスを取るのに何度か実行しないといけないので計算時間が掛かるようになってくるとつらい。

irb(main):004:0> t=Time.now; BD::PI.compare (99.0**2).to_d / (0..520).map{|n| ((4*n).! * (1103+26390*n) ).to_s.to_d.div( (( (4**n) * (99**n) * n.! )**4).to_s.to_d, 4000 ) }.inject(:+) / (2.0.to_d.sqrt(4000)*2); puts Time.now-t
0.314159265358979323846264338327950288419716939937510582097494459230781640628620
＜云々＞
3966655730:4008
90.828125
=> nil

4000桁、520項、1分半。

BigDecimal でもう少し

1万桁くらいまで挑戦しようか。ポイントだった級数各項の割算の桁数指定なんだけど、一律に大きくする必要ないような気もする。ちょっとやってみる。

irb(main):030:0> t=Time.now; BD::PI.compare (99.0**2).to_d / (0..1255).map{|n| ( (4*n).! * (1103+26390*n) ).to_s.to_d.div( (( (4**n) * (99**n) * n.! )**4).to_s.to_d, 10000-n*10000/1255+10 ) }.inject(:+) / (2.0.to_d.sqrt(10000)*2); puts (Time.now-t)/60
＜云々＞
785667227966:10011
11.8442708333333
=> nil

「10000-n*10000/1255+10」線型に調整してみた、1万桁1255項、12分。時間表示60で割ってるの注意
ちなみにこの調整しないと26分。

irb(main):111:0> t=Time.now; BD::PI.compare (99.0**2).to_d / (0..1255).map{|n| ( (4*n).! * (1103+26390*n) ).to_s.to_d.div( (( (4**n) * (99**n) * n.! )**4).to_s.to_d, 10000 ) }.inject(:+) / (2.0.to_d.sqrt(10000)*2); puts (Time.now-t)/60
＜云々＞
78566722796:10010
25.9786458333333
=> nil

実際に必要になる桁数は指数的に減少すると想定されるので、ちょっといろいろやってみたんだけど、計算自体やメソッド呼び出しにコストが掛かるのかあんまり違いが出てこない。
さっきの線型な調整でいいや、別に。

ええと

だいたい桁数見るのに本物と比べるのはずるいです。ほんとなら本物はないわけで、計算するときの有効桁数と級数の収束性等からこの桁までは大丈夫な筈だ、と云うべきものでしょう。
まあ、やってみた、と言うだけの計算でした。

そういえば、Ruby会議2009の時に言ってた Decimal (次世代数値演算ライブラリDecimalという再発明の意義 - 日本Ruby会議2009) って何処にあるんだろう。