円周率の計算（Machin の式）

Ruby で遊んでいたら、高校生時代に友達と円周率を計算して遊んでいたことを思い出した。当時はまだコンピュータを持っておらず、関数電卓で円に内接する多角形の面積を計算したり、解析概論から Machin の式を見つけ出して計算したりしていた。友達が PC-9801 を持っていたが、どうやって倍精度計算をするか、その限界を越えるにはどうしたらいいか、といった話をグダグタとしていた。

Machin の式はこんな感じ。

この式が良いのは、冪乗があるものの、加減乗除だけで構成されていて単純なことだ。
ruby で書くと以下のようになる。
ruby は整数の計算は桁数の制限がないことと、この式はひたすら項を加減していくだけなので、全部の項に factor をかけておくだけで、小数の計算精度を考えずに計算できる。

def Machin()
  prec = 100
  matched_index = 0
  loop = 30

  factor = 1
  for i in 1..prec do
    factor *= 10
  end
  factor1 = 16 * factor
  factor2 =  4 * factor
  sum = factor1 / 5 - factor2 / 239

  div = 3
  for i in 1..loop
    sum += -factor1 / (div * 5**div)   + factor1 / ((div+2) * 5**(div+2))
    sum -= -factor2 / (div * 239**div) + factor2 / ((div+2) * 239**(div+2))
    div += 4
    puts "#{sprintf("%8d",i)} #{sum}"
  end
  return sum
end

$start = Time.now
puts "Machin の式による円周率"
Machin()
puts Time.now-$start

古い MacBook Pro (Core i7-4578U 3GHz dual core) で実行してみると意外に速い。1 msec 程度で小数点以下 30桁近く計算できているようにみえる。

$ ruby machin.rb
Machin の式による円周率
       1 31416210293250344250468325171164080697062446182134305548606396697574415082581156863085774715303620079
       2 31415926824043995172402598360735758604897579967201355063017547035914903487748109950086128003694245533
       3 31415926536235547619955045938203118459271322712991732260010405689815976684675012175940708988985388463
       4 31415926535898358474857006722516843955090730354281179593931838398495000890439868309465181373830789786
       5 31415926535897932947473748577153435433787472210278399039415887536036600624560331366617000933453000146
       6 31415926535897932385393245926718652825091820036365355561499505447671429038153835703541996763189559804
       7 31415926535897932384627502076754923578603952858843005625356465856445914191569978625053179643878350850
       8 31415926535897932384626435346265610108418224608598626348719124702975728463343742259229373035461723957
       9 31415926535897932384626433834967774246425946616320634070726846710697736185351464266957897915043585053
      10 31415926535897932384626433832798181320873204152280609055092074978365028242887424241942263143313044890
      11 31415926535897932384626433832795033456017380588539006987908612446065286640820240779409963401710977707
      12 31415926535897932384626433832795028848774340169568717799892980661436541525752937262258421934060087564
      13 31415926535897932384626433832795028841981785615050020762231123465728069457680425238581026447563602177
      14 31415926535897932384626433832795028841971709044323951350064488696988037559753790469840994549636967408
      15 31415926535897932384626433832795028841971694016301271684489234984290065900959680994986868403485258600
      16 31415926535897932384626433832795028841971693993784982166044833812760297891191671226977100393717248832
      17 31415926535897932384626433832795028841971693993751109427006988812698398829976010400673834116408142191
      18 31415926535897932384626433832795028841971693993751058287322259793318396065343944029846128154609955813
      19 31415926535897932384626433832795028841971693993751058209867272216277219504082451389153487461969263172
      20 31415926535897932384626433832795028841971693993751058209749625351706126852156677062114887680753448201
      21 31415926535897932384626433832795028841971693993751058209749446196953453517888686699842356993297743027
      22 31415926535897932384626433832795028841971693993751058209749445923497087355714539445875926316609121046
      23 31415926535897932384626433832795028841971693993751058209749445923078806095979536420509089060197389981
      24 31415926535897932384626433832795028841971693993751058209749445923078165048578878411934919267800531363
      25 31415926535897932384626433832795028841971693993751058209749445923078164064377944349349611767129831896
      26 31415926535897932384626433832795028841971693993751058209749445923078164062864424622439667791268307192
      27 31415926535897932384626433832795028841971693993751058209749445923078164062862093586998418682337244371
      28 31415926535897932384626433832795028841971693993751058209749445923078164062862089991840979712133067678
      29 31415926535897932384626433832795028841971693993751058209749445923078164062862089986288946017364411264
      30 31415926535897932384626433832795028841971693993751058209749445923078164062862089986280361562208912209
0.001099

思ったより早かったので、調子に乗って、桁数を増やしてみた。また、項を加えても変わらない桁までを計算できたとみなして、その桁だけを表示するように変えてみた。

## digit 桁（前の計算でマッチしたところ）から文字列の長さまでの s1, s2 を比較する。
## 計算が一致しているところまでを正しく計算できているとみなして、その桁数を返す
def check_match(s1,s2,digit)   
  matched_index = -1
  for i in digit..s1.length do
    if s1[i] == s2[i]
      matched_index = i
    else
      return matched_index
    end
  end
  return -1
end

def Machin()
  prec = 100000
  matched_index = 0
  loop = 10000

  str = "314"
  
  factor = 1
  for i in 1..prec do
    factor *= 10
  end
  factor1 = 16 * factor
  factor2 =  4 * factor
  sum = factor1 / 5 - factor2 / 239

  div = 3
  for i in 1..loop
    sum += -factor1 / (div * 5**div)   + factor1 / ((div+2) * 5**(div+2))
    sum -= -factor2 / (div * 239**div) + factor2 / ((div+2) * 239**(div+2))
    div += 4
    #puts "#{sprintf("%8d",i)} #{sum}"

    str_prev = str
    str = sum.to_s
    matched_index = check_match(str_prev, str, matched_index)
    puts "#{i} 回目の計算 (#{matched_index}桁 / #{(Time.now - $start)} sec)"
    puts "3.#{str[1..matched_index]}"  #計算結果の見た目を整える
    if matched_index < 0
      break
    end
  end
  return sum
end

$start = Time.now
puts "Machin の式による円周率"
Machin()
puts Time.now-$start

実行してみると、見た目で桁数が増えていくので楽しい。

$ ruby pi-machin.rb
Machin の式による円周率
1 回目の計算 (2桁 / 1.757707 sec)
3.14
2 回目の計算 (3桁 / 1.844319 sec)
3.141
3 回目の計算 (7桁 / 1.924601 sec)
3.1415926
4 回目の計算 (9桁 / 1.999731 sec)
3.141592653
5 回目の計算 (12桁 / 2.077495 sec)
3.141592653589
6 回目の計算 (16桁 / 2.154938 sec)
3.1415926535897932
7 回目の計算 (18桁 / 2.233316 sec)
3.141592653589793238
8 回目の計算 (21桁 / 2.31235 sec)
3.141592653589793238462
9 回目の計算 (24桁 / 2.404382 sec)
3.141592653589793238462643
10 回目の計算 (27桁 / 2.483269 sec)
3.141592653589793238462643383

......

273 回目の計算 (764桁 / 24.741014 sec)
3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721135000

25秒ほどで 760桁ほど計算できた。ただ、結果が合っているのかは確認していません．．．