« 2009年11月 | トップページ | 2010年1月 »

2009年12月

2009年12月27日 (日)

Project Euler : Problem 10 〜 200万以下の素数の和

 問題はこちらをご覧ください。

 

 問題の解答の説明に入る前に、以前の記事の補足から……。
 "Haskell で素数"の記事に rst76 さんが付けてくださったコメントにもあるように、primes の型を [Int] にするとかなり速度が上がることが分かりました。
 Project Euler の問題を解くだけであれば、必要となる素数は Int の範囲で充分だとは思いますが、今回の問題のように最終的な解答が Int の範囲を越えることはありえます。型が [Int] のままでは over flow によって足元をすくわれてしまうこともありそうです。
 そこでまず intPrimes :: [Int] で [Int] 型の素数列を求めて、それを [Integer] に変換する事で primes :: [Integer] を求めることにしてみました。

 では、問題の解答です。

intPrimes :: [Int] intPrimes = [2, 3] ++ sieve [] (tail intPrimes) 5 where sieve divs (x : xs) n = ps ++ sieve (divs ++ [x]) xs (x * x + 2) where isPrime (y : ys) m | rem m y /= 0 = isPrime ys m | otherwise = False isPrime _ _ = True ps = filter (isPrime divs) [n, n + 2 .. x * x - 2] primes :: [Integer] primes = map fromIntegral intPrimes problem010 :: Integer problem010 = sum $ takeWhile (<= 2000000) primes main = print problem010
 Core Solo T1300 1.66GHz(2MB) のパソコンで約 0.67 秒で答えが出ました。primes :: [Integer] の時には約 1.5 秒かかっていたので、型を変えただけでかなりの高速化になっています。

 「エラトステネスの篩」の条件限定版も、素数を生成する時の型を [Int] にするとやはり高速化できます。

primeList :: Int -> [Integer] primeList n = map fromIntegral $ 2 : sieve [3, 5 .. n] where isqrt = (truncate . sqrt . fromIntegral) n sieve (p : xs) | p > isqrt = p : xs | otherwise = p : sieve [x | x <- xs, rem x p /= 0] main = print $ sum $ primeList 2000000
 Core Solo T1300 1.66GHz(2MB) のパソコンで約 1.38 秒で答えが出ました。

 Ruby や Scheme を使っている時には「型」について悩むことはありませんでしたが、Haskell では「型」についてもっとしっかり考えていく必要があるようです。

Project Euler : Problem 9 〜 ピタゴラス数

 問題はこちらをご覧ください。

 

 一番簡単な方法は、リスト内包表記を使って組み合わせを見つけることだと思います。
 a + b + c = n かつ a < b < c という条件から 0 < a < (n/3), a < b < (n/2) ということが分かるので、これを利用して、探索範囲をちょっとだけ絞り込みました。

-- Brute force pythagoreanNums :: Integral a => a -> [(a, a, a)] pythagoreanNums n = [(a, b, c) | a <- [1 .. div n 3], b <- [a + 1 .. div n 2], c <- [n - a - b], a ^ 2 + b ^ 2 == c ^ 2] problem009 :: Integral a => a -> a problem009 n = a * b * c where (a, b, c) = head $ pythagoreanNums n main = print $ problem009 1000
 最適化オプションを付けてコンパイルすると、約 0.06 秒で答えが出ました。
 しかしこの解き方では、a と b の二つの変数を二重ループにして条件に合う組み合わせを探すので、答えが出るまでに時間がかかります。

 以前 Ruby でこの問題を解いたときにも使った方法では、さらに速く答えが出せます。(詳しい解説はこちらをご覧ください)

-- ピタゴラス数 -- a + b + c = n, a^2 + b^2 = c^2, a < b < c の組を探す。 pythagoreanNums :: Integral a => a -> [(a, a, a)] pythagoreanNums n | odd n = [] | otherwise = [(a, b, n - a - b) | a <- [1 .. div n 3], check a, b <- [calc a], a < b] where check a = rem (a * n) (2 * (n - a)) == 0 calc a = div (n * (n - 2 * a)) (2 * (n - a)) problem009 :: Integral a => a -> a problem009 n = a * b * c where (a, b, c) = head $ pythagoreanNums n main = print $ problem009 1000
 これだと、実際に動かす変数は a の一つだけなので、約 0.005 秒で答えが出ます。

Project Euler : Problem 8

 問題はこちらをご覧ください。

 

 1000 桁の数を一行に並べるのに抵抗があったのと、各桁の数字をリストに入れて処理するつもりだったので、まず、1000 桁の数を 20 桁ずつ区切ってリストに入れておきました。
 さらにそのリストをもとに、各桁の数字を要素とするリストを作って、処理をしました。

import Data.Char (digitToInt) nums :: [Integer] nums = [73167176531330624919225119674426574742355349194934, 96983520312774506326239578318016984801869478851843, 85861560789112949495459501737958331952853208805511, 12540698747158523863050715693290963295227443043557, 66896648950445244523161731856403098711121722383113, 62229893423380308135336276614282806444486645238749, 30358907296290491560440772390713810515859307960866, 70172427121883998797908792274921901699720888093776, 65727333001053367881220235421809751254540594752243, 52584907711670556013604839586446706324415722155397, 53697817977846174064955149290862569321978468622482, 83972241375657056057490261407972968652414535100474, 82166370484403199890008895243450658541227588666881, 16427171479924442928230863465674813919123162824586, 17866458359124566529476545682848912883142607690042, 24219022671055626321111109370544217506941658960408, 07198403850962455444362981230987879927244284909188, 84580156166097919133875499200524063689912560717606, 05886116467109405077541002256983155200055935729725, 71636269561882670428252483600823257530420752963450] partProducts :: [Int] -> [Int] partProducts xs'@(x : xs) | length xs' < 5 = [] | otherwise = prd : partProducts xs where prd = product $ take 5 xs' problem008 :: Int problem008 = maximum $ partProducts ns where ns = concatMap dexToList nums dexToList n = map digitToInt $ show n main = print problem008

2009年12月25日 (金)

Project Euler : Problem 7 ~ 10001 番目の素数

 問題はこちらをご覧ください。

 

 前回の記事に出てきた素数列を作る関数を使います。

 まずは「エラトステネスの篩」から。

primes :: [Integer] primes = 2 : sieve [3, 5 ..] where sieve (x : xs) = x : sieve [y | y <- xs, rem y x /= 0] main = print $ primes !! 10000
 約 6.1 秒と、ちょっと時間がかかります。

 

 次に、自分が改良した関数。

primes :: [Integer] primes = [2, 3] ++ sieve [] (tail primes) 5 where sieve divs (x : xs) n = ps ++ sieve (divs ++ [x]) xs (x * x + 2) where isPrime (y : ys) m | rem m y /= 0 = isPrime ys m | otherwise = False isPrime _ _ = True ps = filter (isPrime divs) [n, n + 2 .. x * x - 2] main = print $ primes !! 10000
 約 0.05 秒で答えが出ます。「エラトステネスの篩」の約 120 倍のスピードです。

Haskell で素数

 "Project Euler" の問題には「素数」に関連したものが時々出てきます。そこで今回は、Haskell でどうやって素数列を扱っていけば良いかを考えてみます。

 Haskell で有名なものに「エラトステネスの篩」があります。(以下のコードは最初から 2 以外の偶数を省いたものです)

primes :: [Integer] primes = 2 : sieve [3, 5 ..] where sieve (p : xs) = p : sieve [x | x <- xs, rem x p /= 0]
 「簡単なコードで素数の無限数列を扱える」ということで Haskell のコードとしては有名(?)なものですが、実は非常に遅いです。
 このコードは、"ghc" で最適化のオプションを付けてコンパイルしても、「20万以下の素数の和」を求めるのに約 40 秒もかかってしまいます。(これ以降、時間に関する記述は、最適化オプションを付けてコンパイルしたコードのものです。ちなみに私の使っているノートパソコンのスペックは "Pentium M, 1500MHz" です)
 Project Euler の Problem 10 のように「200万以下の素数の和」を求めるのにはかなりの時間がかかることが予想されます。

 そこで、自分なりに考えたのが次のコードです。やっていることは、単に素数の候補が本当に素数かどうかを判定しているだけです。
 ただし、素数の候補はただの奇数ではなく、"6 * n ± 1" の形をとる奇数だけを集めることによって、ちょっとだけスピードアップを狙っています。

primes :: [Integer] primes = [2, 3] ++ filter isPrime ns where ns = [x + y | x <- [6, 12 ..], y <- [-1, 1]] isqrt = truncate . sqrt . fromIntegral isPrime n = all ((/= 0) . (rem n)) $ takeWhile (<= isqrt n) primes
 このコードでは、「200万以下の素数の和」が約 4.8 秒で出ます。

 

 自分ではこれ以上いいアイデアを考えつかなかったので、ネット上をいろいろと探してみたところ、こちらのホームページで素晴らしいコードを見つけました。(実際の考え方やアルゴリズムについては、直接こちらのホームページをご覧ください)

primes5 = 2 : primes' where primes' = 3 : sieve 0 5 sieve i x = filter isPrime [x,x+2..p*p-2] ++ sieve (i+1) (p*p+2) where (ps,p:_) = splitAt i primes' isPrime x = all ((/=0).rem x) ps
 たぶん、私の知る限りではもっとも優れたアルゴリズムではないかと思います。実際にコンパイルして「200万以下の素数の和」を求めてみると、約 2.18 秒で答えがでました。

 人のコードをそのまま使うのはちょっと悔しい気がしたので、「このコードを土台にして何とかもっといいコードを書けないか」とかなりの時間考えました。
 まずやったことは、"isPrime" 関数で使われている "ps" を毎回 "primes" から切り出しているのが、少し時間の無駄遣いの用な気がしたので、アルゴリズム自体は変えずに、"sieve" に "ps" に相当する数列を保持させるようにしてみました。

primes :: [Integer] primes = [2, 3] ++ sieve [] (tail primes) 5 where sieve divs (x : xs) n = ps ++ sieve (divs ++ [x]) xs (x * x + 2) where isPrime m = all ((/= 0) . (rem m)) divs ps = filter isPrime [n, n + 2 .. x * x - 2]
 実際にコンパイルしてみると……実行時間はほとんど変わりませんでした。毎回 "divs ++ [x]" としているのが処理に負担をかけているのでしょうか?

 次に、"isPrime" で関数を合成して "all" を使って処理しているところを敢えて再帰関数にしてみました。

primes :: [Integer] primes = [2, 3] ++ sieve [] (tail primes) 5 where sieve divs (x : xs) n = ps ++ sieve (divs ++ [x]) xs (x * x + 2) where isPrime (y : ys) m | rem m y /= 0 = isPrime ys m | otherwise = False isPrime _ _ = True ps = filter (isPrime divs) [n, n + 2 .. x * x - 2]
 すると……なんと「200万以下の素数の和」が約 1.95 秒で出ました。抽象化をするとそれなりに代償がいるということでしょうか?
 もしかしたら、再帰関数の方に末尾再帰の最適化が行われたせいで若干速くなったのでしょうか?

 

 ちなみに冒頭で「遅くて使えない」と言った「エラトステネスの篩」ですが、ちょっと手を加えて「n 以下の素数列を求める」といった使い方に限定するとかなり使えるものに変わります。

primeList :: Integer -> [Integer] primeList n = 2 : sieve [3, 5 .. n] where isqrt = (truncate . sqrt . fromIntegral) n sieve (p : xs) | p > isqrt = p : xs | otherwise = p : sieve [x | x <- xs, rem x p /= 0]
 "main = print $ sum $ primeList 2000000" としてコンパイルすると、約 2.7 秒で答えが出ます。

2009年12月19日 (土)

Project Euler : Problem 6

 問題はこちらをご覧ください。

 今回はただ素直に計算しただけですね。

problem006 :: Int -> Int problem006 n = sumSquare n - squareSum n where sumSquare n = (sum [1 .. n]) ^ 2 squareSum n = sum [x ^ 2 | x <- [1 .. n]] main = print $ problem006 100

Project Euler : Problem 5 ~ 最小公倍数

 問題はこちらをご覧ください。

 要するに最小公倍数を求める問題です。
 Haskell には標準ライブラリに「lcm」があるので、

problem005 :: Integral a => a -> a problem005 n = foldr1 lcm [1 .. n] main = print $ problem005 20
で、終わってしまいます。

 

 これではあまりにも簡単過ぎるので、もうちょっと遊んでみましょう。

 では、自分で最小公倍数を求める関数を作ってみます。
 整数 a, b の最大公約数を gcd(a, b)、最小公倍数を lcm(a, b) とすると、lcm(a, b) = a * b / gcd(a, b) となります。そこで、最大公約数を求める "gcd'" を「ユークリッドの互除法」を使って実装し、それを使って "lcm'" を実装してみます。

gcd' m 0 = m gcd' m n = gcd' n (mod m n) lcm' m n = div (m * n) (gcd' m n) main = print $ foldr1 lcm' [1 .. 20]

 

 次に、ちょっと違う方向から考えてみましょう。
 2 〜 20 までの数を素因数分解した場合、すべての素因数をひとつずつ(同じ素因数が複数個ある場合はもっとも次数の大きいもの)を集めていくと答えが出るはずです。

import Data.List (group) -- 素因数分解 factorize :: Integral a => a -> [(a, Int)] factorize 1 = [(1, 0)] factorize n = format $ factorize' n divs where divs = 2 : 3 : [x + y | x <- [6, 12 ..], y <- [-1, 1]] factorize' n xs'@(x : xs) | n < x * x = [n] | rem n x == 0 = x : factorize' (div n x) xs' | otherwise = factorize' n xs format = map (\xs -> (head xs, length xs)) . group maxFactors :: Integral a => [(a, Int)] -> [(a, Int)] -> [(a, Int)] maxFactors xs [] = xs maxFactors xs ((f, i) : ys) | null bs = maxFactors ((f, i) : xs) ys | i > (snd $ head bs) = maxFactors (as ++ ((f, i) : tail bs)) ys | otherwise = maxFactors xs ys where (as, bs) = break (\ x -> fst x == f) xs problem005 :: Integral a => a -> a problem005 n = product $ map (\x -> (fst x) ^ (snd x)) xs where xs = foldr1 maxFactors $ map factorize [2 .. n] main = print $ problem005 20
 やたら長くなりましたが、ちゃんと答えは出ます。
 でも、やっぱり一番最初のコードが一番いいですね。

2009年12月14日 (月)

Project Euler : Problem 4 ~ 回文数

 問題はこちらをご覧ください。

 一番簡単な方法は、すべての 3 桁の数の積を求めていく方法だと思います。

-- 回文数か? isPalindromic :: Show a => a -> Bool isPalindromic n = str == reverse str where str = show n problem004 :: Integer problem004 = maximum [x | x <- prod, isPalindromic x] where prod = [y * z | y <- [100 .. 999], z <- [y .. 999]] main = print problem004

 しかし、この方法では計算量が多いため、時間がかかります。
 コンパイルして実行した場合、私のパソコンでは次のようになりました。

real 0m0.224s user 0m0.188s sys 0m0.008s

 以前、Ruby 版のコードを載せた記事にも書いたのですが、実際に求めるのは最大値だけなので、「3 桁の数の積を求めるための二重ループを作るときに、大きい数から積を求めていって、その時点で分かっている回文数の最大値より積の数が小さくなったらそれ以上計算をしない」という方法をとると、計算量が減らせるのでかなりのスピードアップが望めます。
 そこで、自分が以前書いた Ruby 版や Scheme 版を基にして Haskell 版のコードを書こうとしたのですが、積を求めるための二重ループをどうやって関数の形にするかで結構苦労しました。

-- 回文数か? isPalindromic :: Show a => a -> Bool isPalindromic n = str == reverse str where str = show n nss :: [[Int]] nss = map prd [999, 998 .. 100] where prd x = [x * y | y <- [x, x - 1 .. 100]] problem004 :: Int problem004 = foldl findMax 0 nss where -- その時点の最大値より小さい値は調べない findMax m ns | null ns' = m | otherwise = max m (head ns') where ns' = filter isPalindromic $ takeWhile (> m) ns main = print problem004

 このコードをコンパイルして実行すると、

real 0m0.013s user 0m0.012s sys 0m0.000s
という結果になりました。最初のコードの約 14 倍のスピードです。

 さらに、ちょっと違う視点からの解法を考えてみます。
 それは「まず、大きい順に回文数を作っていき、3 桁の数の積に分解できるものを探す」という方法です。
 999 x 999 = 998001 なので、997799 以下の回文数を探せばいいことになります。

-- 奇数桁の回文数を作る oddPalNum :: Int -> Int oddPalNum n = read (ns ++ tail (reverse ns)) where ns = show n -- 偶数桁の回文数を作る evenPalNum :: Int -> Int evenPalNum n = read (ns ++ reverse ns) where ns = show n -- 3 桁 * 3 桁に分解できるか? findDivs :: Int -> Bool findDivs n = not $ null ns where isqrt = truncate . sqrt . fromIntegral ns = [x | x <- [isqrt n, isqrt n - 1 .. 100], rem n x == 0, div n x >= 100, div n x <= 999] problem004 :: Int problem004 = head [x | x <- pals1 ++ pals2, findDivs x] where pals1 = map evenPalNum [997, 996 .. 100] pals2 = map oddPalNum [997, 996 .. 100] main = print problem004

 Project Euler は同じ問題でもアプローチの仕方を変えるといろいろな解法が考えられるので、一つの問題でいろいろなコードを書くことができます。
 コードを書く練習としてはなかなか良いのでは?と思います。

追記 (2010/05/16):以前書いたコードが冗長に思えたので、一部書き直しました。

2009年12月 3日 (木)

Project euler : Problem 3 ~ 素因数分解

 問題はこちらをご覧ください。

 いわゆる素因数分解の問題ですね。
 というわけで Haskell でも素因数分解をする関数を作ってみました。

 一番単純な試し割り法ですが、そこそこ速いので Project Euler を解くにはこの位で十分かなと思っています。
 後々のことを考えて、結果は素因数と次数をペアにしたものをリストにまとめる形にしました。

import Data.List (group) -- 素因数分解 factorize :: Integral a => a -> [(a, Int)] factorize 1 = [(1, 0)] factorize n = format $ factorize' n divs where divs = 2 : 3 : [x + y | x <- [6, 12 ..], y <- [-1, 1]] factorize' n xs'@(x : xs) | n < x * x = [n] | rem n x == 0 = x : factorize' (div n x) xs' | otherwise = factorize' n xs format = map (\xs -> (head xs, length xs)) . group main = print $ factorize 600851475143

 Ruby や Scheme で作ったものと比べると、コードの長さが約半分位になっています。

Project euler : Problem 2 〜 フィボナッチ数

 問題はこちらをご覧ください。

 皆様ご存知の Fibonacci 数の問題ですが、この問題で提示されているフィボナッチ数は "1, 1, 2, 3, 5 ..." といった一般的なものではなく、"1, 2, 3, 5 ..." といったふうに、一つずれた形になっています。(まあ、答えには影響はありませんが……)
 まずは、再帰的な定義によるものから……同じ計算を繰り返しやっているので遅いです。

fibonacci :: Integer -> Integer fibonacci 1 = 1 fibonacci 2 = 2 fibonacci n = fibonacci (n - 2) + fibonacci (n - 1) problem002 :: Integer -> Integer problem002 n = sum [x | x <- fibList, even x] where fibList = takeWhile (< n) $ map fibonacci [1 ..] main = print $ problem002 4000000

 次は反復的にフィボナッチ数列を求めるもの……これはそこそこ速いのでは?

fibonacci :: [Integer] fibonacci = fib 1 2 where fib a b = a : fib b (a + b) problem002 :: Integer -> Integer problem002 n = sum [x | x <- takeWhile (< n) fibonacci, even x] main = print $ problem002 4000000

 Haskell では比較的簡単な定義で無限の数列を扱えるのが魅力的ですね。

2009年12月 2日 (水)

Project Euler : Problem 1 ~ 今度は Haskell で

 先日まで、Project Euler を Scheme でやり直そうと思っていたんですが、Haskell の練習ということで、しばらくは Haskell でやり直していきます。ではさっそく第 1 問から……

 問題はこちらをご覧ください。

 問題自体は簡単なので、いろんな方法を考えてみました。

-- 補助関数によるループ main = print $ problem001 999 problem001 :: Integer -> Integer problem001 n = iter n 0 where iter m sum | m == 0 = sum | rem m 3 == 0 = iter (m - 1) (sum + m) | rem m 5 == 0 = iter (m - 1) (sum + m) | otherwise = iter (m - 1) sum

-- List によるループの抽象 main = print $ problem001 999 problem001 :: Integer -> Integer problem001 n = sum $ filter check [1 .. n] where check x = mod x 3 == 0 || mod x 5 == 0

-- リスト内包表記 main = print $ problem001 999 problem001 :: Integer -> Integer problem001 n = sum [x | x <- [1..n], (mod x 3 == 0 || mod x 5 == 0)]

-- 等差数列を利用 main = print $ problem001 999 problem001 n = sum1 + sum2 - sum3 where sum1 = sum [0, 3 .. n] sum2 = sum [0, 5 .. n] sum3 = sum [0, 15 .. n]

 最初は「なんだか、難しそうだなぁ……」と思っていましたが、実際にコードを書いてみると、(「純粋」とは言えないまでも)関数言語の仲間である Scheme と似た感覚でコードがかけたので、少しだけ苦手意識がなくなりました。
 Haskell と Scheme では lazy と eager という違いがあるので、まったく同じようにはいかないとは思いますが、Scheme で書いたコードを参考にすれば、少しは楽に Haskell のコードが書けるような気がします。

« 2009年11月 | トップページ | 2010年1月 »

2016年7月
          1 2
3 4 5 6 7 8 9
10 11 12 13 14 15 16
17 18 19 20 21 22 23
24 25 26 27 28 29 30
31            
フォト

最近のトラックバック

無料ブログはココログ