Project Euler 解答

Project Euler Problem 100

Project Euler Problem 100

```import Data.Ratio
import Data.List
import Data.Maybe

-- チャンドラグプタの方法
gupta1 :: Num t => t -> (t, t) -> (t, t)
gupta1 n (a,b) = (a*a+n*b*b,2*a*b) -- gupta identity self application
gupta2 :: Num t => t -> (t, t) -> (t, t) -> (t, t)
gupta2 n (a,b) (c,d) = (a*c+n*b*d,a*d+b*c)
lis :: [(Integer, Integer)]
lis = iterate (gupta2 2 (3,1)) (3,1)

-- chakravala methodの実装
chak :: Integral t => t -> [(t, t)]
chak n = iterate (gupta2 n (numerator a',numerator b')) (numerator a',numerator b')
where
(a',b',c') = imp seed
lim = fromJust \$ find (\x -> (x^2)>n) [1..]
seed = findfirst n
gupta (a,b,k) = (a*a+(n%1)*b*b,2*a*b,k*k) -- gupta identity self application
imp x@(a,b,k) | k == 1%1 = (a,b,k)
| k == (-1)%1 = gupta x
| otherwise = imp \$ findabkm n x
-- 自明な初期値を求める。b=1固定でa,kを設定
findfirst n = (a%1,1%1,k%1)
where
a = snd \$ minimum [(abs(x^2 - n),x)|x<-[1..(lim+1)]]
k = a^2 - n
-- chakravala method繰り返し部分
findabkm n (a,b,k) = (aa,bb,kk)
where
list = filter (\x -> (numerator (a+b*(x%1)))`rem`(numerator k)==0) [1..(lim+1)]
m = snd \$ minimum [(abs(x^2 - n),x%1)|x<-list]
aa = (a*m+(n%1)*b)/(abs k)
bb = (a+b*m)/(abs k)
kk = (m*m-(n%1))/k
--(3,1)が初期解のペル値を調べる
pells :: [((Integer, Integer), Integer)]
pells = takeWhile ((/=(3,1)).fst) \$ map (\(l,r) -> (head (chak l),r)) \$ filter (\(x,y) -> x /= (round \$ sqrt \$ fromIntegral x)^2) (zip [2..] [2..])
-- n=8となる
-- 後はrが求まったときにbを求める関係式を求めれば終わり
lis2 :: [(Integer, Integer)]
lis2 = iterate (gupta2 8 (3,1)) (3,1)
-- (2086968960817,313506783024) = (blue,red)
blues :: [Integer]
blues = tail \$ map (\(t,r) -> (t - 2*r + 1)`quot`2) lis2
anslist :: [(Integer, Integer)]
anslist = zip blues (map snd lis2)
ans1imp :: (Integer, Integer)
ans1imp = anslist !! ((fromJust \$ findIndex (\(l,r) -> l+r>10^12) anslist))
ans1 :: Integer
ans1 = let (a,_) =ans1imp in a
-- 756872327473
main :: IO ()
main = print ans1
```
• はじめに
• プロジェクトオイラー問題
• リンク等