ペル方程式

数論よくわからん…
いや説明聞いたらたぶんその場では分かった気になるんだろうけどすぐに忘れてしまう.
よくわからんものは結果だけ暗記しちゃった方が実用的には良いのかもしれないなーとかそんなことを思うこの頃.
いや,「実用」ってなんだ.
まあそんなことよりも

ペル方程式って

x^2-ny^2=1
という方程式で,nが平方数じゃない自然数(2, 3, 5, 6, 7, 8, 10,...)のときに整数解を求めましょうという問題.
ペル方程式は自明な解(x=1, y=0)以外に必ず解を持っていて,ある解(x, y)について
x_k+y_k\sqrt{n}=(x+y\sqrt{n})^k
として得られるx_ky_kもペル方程式の解になる.それどころか(x, y)を最小解とすれば,全ての解を上の式から求めることができる.そして,
\alpha=x+y\sqrt{n}
\beta=x-y\sqrt{n}
とすると,
x_k=(\alpha^k+\beta^k)/2
y_k=(\alpha^k-\beta^k)/(2\sqrt{n})
としてk番目の解を求められる.以上Wikipediaの焼き直し.

よく分からないのでやってみよう

n=5の場合の最小解は(x, y) = (9, 4)らしいです.

## n = 5の場合のペル方程式の解判定
Pelln5 <- function(x, y){
  x^2 - 5*y^2 == 1
}
## > Pelln5(9, 4)
## [1] TRUE

それで,最小解が分かってれば全ての解を生成できるはずです.

## n = 5の場合のペル方程式の解を生成
## n = 5のときの最小解は(x, y) = (9, 4)
x <- 9
y <- 4
n <- 5
a <- x + y * sqrt(n)
b <- x - y * sqrt(n)
## k番目の解
answer <- function(k){
  xk <- (a^k + b^k)/2
  yk <- (a^k - b^k)/(2*sqrt(n))
  return(c(xk, yk))
}
## > answer(1)
## [1] 9 4
## > answer(2)
## [1] 161  72

しかもこれで全ての解を網羅できるということなので,いくつかの解を力技でもとめて比較してみる.

## Brute force
limit <- 10000
ans.x <- ans.y <- numeric(0)
for(i in 1:limit){
  if(sum(hit <- Pelln5(i, 1:limit)) != 0){
    ans.x <- c(ans.x, i)
    ans.y <- c(ans.y, (1:limit)[hit])
  }
}

## チカラワザで求めた解
data.frame(ans.x, ans.y)
## > data.frame(ans.x, ans.y)
##   ans.x ans.y
## 1     9     4
## 2   161    72
## 3  2889  1292

## 公式から求めた解
matrix(answer(1:3), ncol=2)
## > matrix(answer(1:3), ncol=2)
##      [,1] [,2]
## [1,]    9    4
## [2,]  161   72
## [3,] 2889 1292

おー(理屈はワカラン)