Problem 66

Problem 66 - Project Euler
ちゃんと連分数の循環節を掃きだす関数を書いた。あとでそこだけ詳しくメモる。
ペル方程式は理解しきれなかったのでかなり答えに近いものを調べてしまった…
だというのに実装するのにかなり苦労した。そしてまたgmpに頼ってしまった。
ううむ。ダメダメだなぁ。
ちなみに調べたのはここ。

• http://www004.upp.so-net.ne.jp/s_honma/pell/pell.htm
• ペル方程式の解の列挙方法 - まめめも

かなり参考になりました。

## 連分数の循環節を掃きだす
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 <- l-p*a
  p1 <- p <- (n-l^2)/p
  om <- (om0 - l)/p
  i <- 2
  repeat{
    i <- i+1
    a <- rep.set[i] <- floor(om)
    l <- - l - a*p
    p <- (n - l^2)/p
    om <- (om0 - l)/p
    if(l1==l  && p1==p) break
  }
  return(rep.set[-length(rep.set)])
}

## ペル方程式の最小解
pel.min <- function(n){
  if(sqrt(n)%%1==0) return(0)
  qn <- rep.frac(n)
  if(length(qn)>2) qn <- qn[-length(qn)]
  x_1 <- as.bigz(1)
  x <- as.bigz(qn[1])
  y_1 <- as.bigz(0)
  y <- as.bigz(1)
  for(i in 2:length(qn)){
    x_2 <- x_1
    x_1 <- x
    x <- qn[i]*x_1 + x_2
    y_2 <- y_1
    y_1 <- y
    y <- qn[i]*y_1 + y_2
  }
  if(x^2 - n*y^2 == 1){
    return(c(x, y))
  }else{
    return(c(pow.bigz(x,2)+n*pow.bigz(y,2), 2*x*y))
  }
}

max.min <- 0
ans <- 0
candi <- as.bigz(1)
for(i in 1:1000){
  candi <- pel.min(i)[1]
  if(candi > max.min){
    max.min <- candi
    ans <- i
  }
}
ans