平方根を連分数展開したときの循環節を掃きだす

Project EulerのProblem 66を解く際に平方根の連分数を求める関数を書きました。これについてちょっとメモを。
それと当然ながらProject EulerのProblem 64〜66あたりのネタバレになっちゃいます。注意して下さい。

## 連分数の循環節を掃きだす
rep.frac <- function(n){
  if(sqrt(n)%%1==0) return(0)
  om0 <- om <- sqrt(n)
  rep.set <- numeric(0)
  ## a_0
  a <- rep.set[1] <- floor(om)
  p <- n-a^2
  l <- a
  om <- (om0 + l)/p
  ## 循環節のひとつめ(ここに戻ったら終了)
  a <- rep.set[2] <- floor(om)
  l1 <- l <- p*a - l
  p1 <- p <- (n-l^2)/p
  om <- (om0 + l)/p
  i <- 2
  repeat{
    i <- i+1
    a <- rep.set[i] <- floor(om)
    l <- a*p - l
    p <- (n - l^2)/p
    om <- (om0 + l)/p
    if(l1==l  && p1==p) break
  }
  return(rep.set[-length(rep.set)])
}

連分数展開

ある数 $\omega$ について、次の形の分数にすることを連分数展開といいます。
$\omega = a_0 + \frac{1}{a_1 + \frac{1}{a_2 + \frac{1}{a_3 + \frac{1}{a_4 + \cdots}}}}$
それで、平方根も含めて二次無理数と呼ばれる無理数を連分数展開した場合、 $a_0, a_1, a_2,...$ からなる数列が循環節を持つということが知られています(参考：http://aozoragakuen.sakura.ne.jp/suuronN/node78.html)。
$a_0, a_1, a_2...$ を求めること自体はそれほど難しいことではありません。
ある数 $\omega$ が与えられたとしましょう。まず、 $\omega$ を越えない最大の整数が $a_0$ です。次に、 $\omega$ から $a_0$ を引いた残り、この逆数を $\omega_1$ とし、これを越えない最大の整数が $a_1$ になります。つまり、

$\omega_n$ を越えない最大の整数を $a_n$ とする
$\omega_n - a_n$ の逆数を $\omega_{n+1}$ とする

の繰り返しによって $a_n$ を順次決定していくわけです(参考：連分数 - Wikipedia)。
循環節がそれほど長くならない場合はこの方法で問題ないのですが、ある程度以上操作を繰り返すとあっというまに誤差が拡大して答が不正確なものとなってしまいます。

平方根の連分数展開

精度を維持するためには極力計算を整数のままで行なうことがポイントになります。ある数 $N$ についてその平方根を連分数展開し、 $a_0,a_1,a_2...$ を求めていく方法を考えてみます。
まず、
$\omega_0 = \sqrt{N}$
とおきます。次に、上記の手順に従って、
$a_0 = \lfloor\omega_0\rfloor$
$\omega_1 = (\omega_0 - a_0)^{-1}$
ここで $\omega_1$ を有理化し、整数にできる部分はまとめて別の文字に置き換えてしまいます。
$(\omega_0 - a_0)^{-1} = \frac{\omega_0+a_0}{N - a_0^2} = \frac{\omega_0 + l_1}{p_1}$
この式の最後の形式で今後の $\omega_n$ を表現していくことにしましょう。
もう一段階計算を進めます。
$a_1 = \lfloor\omega_1\rfloor$
$\array{\omega_1 &=& (\omega_1 - a_1)^{-1}\\ &=& \left{\frac{\omega_0 + l_1}{p_1} - a_1\right}^{-1}\\ &=& \left{\frac{\omega_0 - (p_1a_1-l_1)}{p_1}\right}^{-1}\\ &=& \left{\frac{\omega_0 - l_2}{p_1}\right}^{-1}\\ &=& \frac{\omega_0 + l_2}{(N-l_2^2)/p_1}\\ &=& \frac{\omega + l_2}{p_2}$
これより、次の漸化式が得られます。
$a_n = \lfloor\omega_n\rfloor$
$l_{n+1} = p_na_n - l_n$
$p_{n+1} = \frac{N-l_{n+1}^2}{p_n}$
$\omega_{n+1} = \frac{\omega_0 + l_{n+1}}{p_{n+1}}$
ただしここで
$l_0 = 0$
$p_0 = 1$
とします。
ここで導入した $l$ 、 $p$ はいずれも整数で、漸化式の計算に必要なパラメータも全て整数ですから、誤差の伝播を気にすることなく計算を進めることができます。これをRで実装したのが最初に示した関数というわけです。