2012年5月2日水曜日

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

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

■ やり方
だいたいこんな感じ
g-1=0g0=0gn+1 = -gn + knhn
h-1=0h0=1hn+1 = (D - gn+12)/hn
k-1=0k0=[√D]kn+1 = [(k0 + gn)/hn]
x-1=0x0=1xn+1 = gnyn + hnyn-1
y-1=0y0=0yn+1 = yn-1 + knyn
で、gn とgn+1 で同じ値が連続したら、
X = (xn2 + Dyn2)/hn
Y = 2xnyn/hn
hn とhn+1 で同じ値が連続したら、
X' = (xnxn+1 + Dynyn+1)/hn
Y' = (xnyn+1 + xn+1yn)/hn
X = X'2 + DY'2
Y = 2X'Y'

■ 例
D = 23の場合
nghkxy
-100001
001410
147141
232351
3371194
24 = (52 + 23*12)/2
5 = 2*5*1/2
D = 29の場合
nghkxy
-100001
001510
154251
2351112
32511613
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でもすぐに答えが出てくる。