数論よくわからん…
いや説明聞いたらたぶんその場では分かった気になるんだろうけどすぐに忘れてしまう.
よくわからんものは結果だけ暗記しちゃった方が実用的には良いのかもしれないなーとかそんなことを思うこの頃.
いや,「実用」ってなんだ.
まあそんなことよりも
ペル方程式って
という方程式で,が平方数じゃない自然数(2, 3, 5, 6, 7, 8, 10,...)のときに整数解を求めましょうという問題.
ペル方程式は自明な解(x=1, y=0)以外に必ず解を持っていて,ある解(x, y)について
として得られるともペル方程式の解になる.それどころか(x, y)を最小解とすれば,全ての解を上の式から求めることができる.そして,
とすると,
として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
おー(理屈はワカラン)