http://projecteuler.net/index.php?section=problems&id=44
まず、できるだけ数学を使わない方法で書いた。Plをlが小さいほうから調べ、Pkとの和が五角数であるものを探す。ただし、kはPk+1-Pk ≤ Plを満たす限り。その中でPl + 2Pkが五角数であるものを探す。きれいに書けた。
int_sqrt :: Int -> Int
int_sqrt n = floor (sqrt (fromIntegral n))
is_pentagonal n = m * m == k && mod (1 + m) 6 == 0 where
k = 1 + 24 * n
m = int_sqrt kpentagonal = [ div (n * (3 * n - 1) ) 2 | n <- [1..] ]
queue p (q:qs) | (head qs) - q > p = []
| otherwise = [(p,q)] ++ (queue p qs)
queues (p:ps) (q:qs) = (queue p (q:qs)) ++ (queues ps qs)
is_valid (p,q) = is_pentagonal (p + q) && is_pentagonal (p + q * 2)
d = head (filter is_valid (queues pentagonal (tail pentagonal)))
main = print(d)
しかし返ってこなかった。
仕方がないので少し数学的にアプローチを。
Pj - Pk = Pl
j(3j - 1) - k(3k - 1) = l(3l - 1)
(j - k)(3(j + k) - 1) = l(3l - 1)
lと3l-1は互いに素なので、個別に素因数分解して結合する。そうすると簡単にその約数が出るので、この条件を満たすj,kが求められる。
primes = 2:[3,5..]
div_pow n p | mod n p /= 0 = (n, 0)
| otherwise = (\(n, e) -> (n, e + 1)) (div_pow (div n p) p)
factorize 1 ps = []
factorize n (p:ps)
| p * p > n = [(n,1)]
| otherwise = if e > 0 then [(p, e)] ++ f else f where
f = factorize m ps
(m, e) = div_pow n p
fact n = factorize n primes
divs [] = [(1,1)]
divs ((p,e):fs) = [ (p^e1 * m, p^(e-e1) * n) |
e1 <- [0..e], (m,n) <- (divs fs) ]pentagonal n = div (n * (3 * n - 1)) 2
convert (d1,d2) = (pentagonal k, pentagonal (k - d1))
where k = div (d1 + (div (d2 + 1) 3)) 2
is_valid (d1,d2) = mod (d1 + d2) 2 == 1
&& mod d2 3 == 2 && d2 + 1 - d1 * 3 > 0
a = map (\n -> filter is_valid
(divs (fact n ++ fact (3 * n - 1)))) [1..]
b = map convert (concat a)int_sqrt :: Int -> Int
int_sqrt n = floor (sqrt (fromIntegral n))
is_pentagonal n = m * m == k && mod (1 + m) 6 == 0 where
k = 1 + 24 * n
m = int_sqrt k
c = filter (\(x,y) -> is_pentagonal (x + y)) b
main = print(head (map (\(x,y) -> x - y) c))
