「今週のお題：円周率に近似できる分数」問題解答 ( CodeIQ )

CodeIQ 数学

はじめに

CodeIQという所で週次で出題される「今週のお題」シリーズに解答しましたので、ネタバレです。

codeiq.jp

問題

今回は次のような問題でした。

  • 入力で与えられる n ( $n\le 11$ ) に対して、「円周率$\pi$と小数点n桁までが一致する分数 ( 有理数 )」の内、最小の分母を出力する。

おそらく「n桁までが一致」ではなく、「n桁以上が一致」としないと無理が出るような気がしますので、そのように解釈して解いています。

解説

さて。「最小の分母」ですから、尺取り法で分母・分子を相互に増やしていって、初めてn桁一致になるのがいつになるか見ても良いのですが、ちょうど「これは連分数でイけるな」という感覚がありました。

円周率と連分数

連分数の話については、例えば円周率.jpの連分数のページに説明があるのですが、 $$ \pi = 3+\dfrac{1}{7+\dfrac{1}{15+\dfrac{1}{1+\dfrac{1}{292+\cdots}}}} $$ というように無限に続く分数で$\pi$を表すことができること、これを途中で打ち切ることで、$\pi$に極めて近い有理数を得ることができること ( 先まで延ばすほど$\pi$に近づく ) ということが特徴です。

連分数に現れる数字は、数列$a_1=\pi,\,a_{n+1}=\frac{1}{a_n-\lfloor a_n \rfloor}$に対して、$\lfloor a_n \rfloor$ によって求めることができます。

さて、これを途中で打ち切る訳ですが、「これと別系統で$\pi$に近づく分数があるのでは?」という疑問が湧いてくるかもしれません。これについてはちょっと説明が上手くまとまらないので割愛しますが、連分数の系統以外で近づくものはない、と言えます。

後は精度の問題です。単純に打ち切る段数を変えると、 $$ \begin{eqnarray} 3+\dfrac{1}{7+\dfrac{1}{15+\dfrac{1}{1}}} &=& \frac{355}{113} &=& 3.141592\cdots\\ 3+\dfrac{1}{7+\dfrac{1}{15+\dfrac{1}{1+\dfrac{1}{292}}}} &=& \frac{103993}{33102}&=& 3.141592653\cdots \\ \end{eqnarray} $$ と、一気に一致する桁数が跳ね上がるところがあります。

しかし、それでは指定された桁数に対して過剰な精度である、言い換えると、もっと小さな分母で十分なのに、大きな分母の数値を計算してしまう可能性があります。この例では292を素直に使ったのが過剰な精度を引き起こした原因であり、最適な精度に調整するには、292の代わりにもうちょっと小さい値を探すことになります。

実装

ということで、提出版のRubyの実装です。

連分数を途中で打ち切った分数の分子・分母は、連分数の各係数$q_n$に対して、$x_n=q_n\cdot x_{n-1}+x_{n-2}$という漸化式の数列$x_n$ として計算できますので、そのように漸化式を回しています。連分数自体を求めるのは標準の浮動小数点演算でやってますが、もっと精度を要求されるなら、任意精度の数値を使う必要が出て来るでしょう。

最後、より小さい分母で最適な精度となるよう調整を行うために bsearch での探索を行います。採用する数字が十分大きければ精度は確保できるため、find-minimumモードでの探索になります。

p ->n{
  uc,up,dc,dp,r,b=3,1,1,0,Math::PI-3,10**n
  t=(Math::PI*b).floor
  loop{
    q,r=(1/r).divmod(1.0)
    ut,dt=up+uc*q,dp+dc*q
    return dp+dc*(1..q).bsearch{|e|(up+uc*e)*b/(dp+dc*e)==t} if ut*b/dt==t
    uc,up,dc,dp=ut,uc,dt,dc
  }
}[gets.to_i]

終わりに

今回は、あじさんからこのようなtweetが出ていて、果たして同じ武器を使っているのか…というのが気になっていたのですが。どうなんでしょうか。

さて。随分と間が空いてしまいました。が、なるべく余裕があれば続けていきたいと思います。今回は再開ということで。