Pell方程式を解くための実装メモ

Project Euler をやってると、Pell方程式 X²-DY²=1 に関連する問題をたまに見かける。これを解くとき、連分数を使うやり方がよく使われるようだけど、半分位の手数で解が得られるやり方がここで紹介されていたので、Haskell で実装してみた。

■ やり方

だいたいこんな感じ

g-1=0	g0=0	gn+1 = -gn + knhn
h-1=0	h0=1	hn+1 = (D - gn+12)/hn
k-1=0	k0=[√D]	kn+1 = [(k0 + gn)/hn]
x-1=0	x0=1	xn+1 = gnyn + hnyn-1
y-1=0	y0=0	yn+1 = yn-1 + knyn

で、g_n とg_n+1 で同じ値が連続したら、

X = (xn2 + Dyn2)/hn

Y = 2xnyn/hn

h_n とh_n+1 で同じ値が連続したら、

X' = (xnxn+1 + Dynyn+1)/hn

Y' = (xnyn+1 + xn+1yn)/hn

X = X'2 + DY'2

Y = 2X'Y'

■ 例

D = 23の場合

n	g	h	k	x	y
-1	0	0	0	0	1
0	0	1	4	1	0
1	4	7	1	4	1
2	3	2	3	5	1
3	3	7	1	19	4

24 = (52 + 23*12)/2

5 = 2*5*1/2

D = 29の場合

n	g	h	k	x	y
-1	0	0	0	0	1
0	0	1	5	1	0
1	5	4	2	5	1
2	3	5	1	11	2
3	2	5	1	16	13

70 = (11*16 + 29*2*3)/5

13 = (11*3 + 16*2)/5

9801 = 702 + 29*132

1820 = 2*70*13

■ コード

pell d=pell' d (0,0,0,0,1) (0,1,floor(sqrt(fromIntegral d)),1,0) 
  where 
   pell' d (g,h,k,x,y) n@(gn,hn,kn,xn,yn)
    | g'== gn    = let xr = div (xn^2 + d*yn^2) hn
                       yr = div (2*xn*yn) hn
                   in (xr, yr)
    | h'== hn    = let xr = div (x'*xn + d*y'*yn) hn
                       yr = div (x'*yn + xn*y') hn
                   in (xr^2 + d*yr^2, 2*xr*yr)
    | otherwise  = pell' d n (g',h',k',x',y')
    where k0 = floor$sqrt$fromIntegral d
          g' = (-gn) + kn*hn
          h' = div (d - g'^2) hn
          k' = div (k0 + g') h'
          y' = y + kn*yn
          x' = g'*y' + h'*yn

■ 結果

ghci> let isNotSquare n= (round . sqrt $ fromIntegral n) ^ 2 /= n
ghci> mapM_ (\n->print$(n,pell n)) $filter isNotSquare [2..15]
(2,(3,2))
(3,(2,1))
(5,(9,4))
(6,(5,2))
(7,(8,3))
(8,(3,1))
(10,(19,6))
(11,(10,3))
(12,(7,2))
(13,(649,180))
(14,(15,4))
(15,(4,1))
ghci> maximumBy (\a b->(on compare (snd.snd))a b )$ map (\n->(n,pell n)) $filter isNotSquare [2..2000]
(1621,(6298101812493732343034974500091457815529942308667051412857352310169665125001,156429324369979112128445583345098338627552043874824108399177922442751050500))

Problem 66 で、D≦1000でXが最大になるものを求める問題があるけど、D≦10000でもすぐに答えが出てくる。