« 2011年1月 | トップページ | 2011年3月 »

2011年2月

2011年2月27日 (日)

Project Euler : Problem 59 ~ 暗号解読

 問題はこちらをご覧ください。
 また、自作の "ForEuler module" に関してはこちらをご覧ください。

 

 今回の問題の解法は、以前 Ruby で解いた時と同じです。考え方についてはこちらをご覧下さい。

{- * キーワードを探す : problem059_find_key.hs -} import Data.Bits (xor) import Data.List (transpose, sort, group, sortBy) import Data.Ord (comparing) import Data.Char (chr, ord) main :: IO () main = do ss <- readFile "cipher1.txt" print $ map f $ map top5 $ divide 3 $ read ("[" ++ ss ++ "]") where f ns = deemChar 'e' ns -- cs を n 個おきにグループに分る -- ex : divide 3 [1..10] => [[1,4,7,10],[2,5,8],[3,6,9]] divide :: Int -> [Int] -> [[Int]] divide n cs = (transpose . f) cs where f [] = [] f xs = as : f bs where (as, bs) = splitAt n xs -- ns の要素のうち出現頻度の高い順に 5 個をリストにして返す top5 :: [Int] -> [Int] top5 ns = map fst $ take 5 $ reverse $ sortBy (comparing snd) lst where lst = [(head xs, length xs) | xs <- (group . sort) ns] -- 鍵文字を c と仮定して、ns を復元する(結果は String として返す) deemChar :: Char -> [Int] -> [String] deemChar c ns = map f ns where f n = [chr $ xor n $ ord c]
{- * 暗号文を復元する : problem059_decode.hs * コンパイル後のファイル名を "decode"、キーワードを "***" と仮定すると、 > decode *** で、暗合文を復元する。 -} import Data.Bits (xor) import Data.Char (ord, chr) import System (getArgs) main :: IO () main = do args <- getArgs ss <- readFile "cipher1.txt" print $ decode (read ("[" ++ ss ++ "]")) (head args) {- -- *** をキーワードに置き換える main :: IO () main = do ss <- readFile "cipher1.txt" print $ decode (read ("[" ++ ss ++ "]")) "***" -} decode :: [Int] -> String -> String decode ns pw = map chr $ zipWith (xor) ns ms where ms = map ord $ concat $ repeat pw

2011年2月23日 (水)

基数ソート

 こちらのブログを見て知ったのですが、「基数ソート」というソート法があるんですね。
 面白そうだったので、 Wikipedia を参考に自分でもコードを書いてみました。(基本方針は、Wikipedia に従って、バケットソートと組み合せることにしました。また、ソートする要素は十進数の整数に限定しました)

 バケットソートに関して二通りの実装を思いついたのですが、まずはタイピング量が少いほうから……。

import Data.List (replicate) radixSort1 :: Integral a => [a] -> [a] radixSort1 xs = loop 1 xs where maxData = maximum xs loop d xs | d > maxData = xs | otherwise = loop (d * 10) $ bucketSort d xs $ replicate 10 [] bucketSort _ [] yss = concatMap reverse yss bucketSort d (x : xs) yss = bucketSort d xs $ ass ++ (x : bs) : bss where i = fromIntegral (rem (div x d) 10) (ass, bs : bss) = splitAt i yss
 試しに 0 〜 999 の範囲の数がランダムに 15000 個入っているリストをソートさせてから合計を求めさせたら、「Core2 Quad, 2.83GHz, Windows 7」という環境で
>ghc -O -o radixSort1 radixSort1.hs --make ghc -O -o radixSort1 radixSort1.hs --make Linking radixSort1.exe ... >radixSort1 +RTS -sstderr radixSort1 +RTS -sstderr radixSort1 +RTS -sstderr 7418991 24,968,552 bytes allocated in the heap 5,379,064 bytes copied during GC 695,808 bytes maximum residency (6 sample(s)) 25,708 bytes maximum slop 3 MB total memory in use (0 MB lost due to fragmentation) Generation 0: 40 collections, 0 parallel, 0.02s, 0.02s elapsed Generation 1: 6 collections, 0 parallel, 0.00s, 0.00s elapsed INIT time 0.00s ( 0.00s elapsed) MUT time 0.02s ( 0.02s elapsed) GC time 0.02s ( 0.02s elapsed) EXIT time 0.00s ( 0.00s elapsed) Total time 0.03s ( 0.04s elapsed) %GC time 50.0% (57.9% elapsed) Alloc rate 1,600,548,205 bytes per MUT second Productivity 50.0% of total user, 41.0% of total elapsed
という結果になりました。
 「あまり速くないのでは?」と思っていたのですが、以前書いた記事の「クイックソート」や「マージソート」に迫る速さだったので、ビックリしました。
 ただ、"ass ++ (x : bs) : bss" といったところや、"(ass, bs : bss) = splitAt i yss" といったところで頻繁にリスト操作が行われているので、これをなんとかすればもっと速くなる可能性があります。

 そこで、冗長なバージョンのバケットソートを試してみました。

radixSort2 :: Integral a => [a] -> [a] radixSort2 xs = loop 1 xs where maxData = maximum xs loop d xs | d > maxData = xs | otherwise = loop (d * 10) $ bucketSort d xs [] [] [] [] [] [] [] [] [] [] bucketSort _ [] ys0 ys1 ys2 ys3 ys4 ys5 ys6 ys7 ys8 ys9 = concatMap reverse [ys0, ys1, ys2, ys3, ys4, ys5, ys6, ys7, ys8, ys9] bucketSort d (x : xs) ys0 ys1 ys2 ys3 ys4 ys5 ys6 ys7 ys8 ys9 | y == 0 = bucketSort d xs (x : ys0) ys1 ys2 ys3 ys4 ys5 ys6 ys7 ys8 ys9 | y == 1 = bucketSort d xs ys0 (x : ys1) ys2 ys3 ys4 ys5 ys6 ys7 ys8 ys9 | y == 2 = bucketSort d xs ys0 ys1 (x : ys2) ys3 ys4 ys5 ys6 ys7 ys8 ys9 | y == 3 = bucketSort d xs ys0 ys1 ys2 (x : ys3) ys4 ys5 ys6 ys7 ys8 ys9 | y == 4 = bucketSort d xs ys0 ys1 ys2 ys3 (x : ys4) ys5 ys6 ys7 ys8 ys9 | y == 5 = bucketSort d xs ys0 ys1 ys2 ys3 ys4 (x : ys5) ys6 ys7 ys8 ys9 | y == 6 = bucketSort d xs ys0 ys1 ys2 ys3 ys4 ys5 (x : ys6) ys7 ys8 ys9 | y == 7 = bucketSort d xs ys0 ys1 ys2 ys3 ys4 ys5 ys6 (x : ys7) ys8 ys9 | y == 8 = bucketSort d xs ys0 ys1 ys2 ys3 ys4 ys5 ys6 ys7 (x : ys8) ys9 | y == 9 = bucketSort d xs ys0 ys1 ys2 ys3 ys4 ys5 ys6 ys7 ys8 (x : ys9) where y = rem (div x d) 10
 いや〜ひどいですね。"y == ..." が 10 個も並んでいます。
 こちらも "radixSort1.hs" と同様のテストをしてみました。
>ghc -O -o radixSort2 radixSort2.hs --make ghc -O -o radixSort2 radixSort2.hs --make Linking radixSort2.exe ... >radixSort2 +RTS -sstderr radixSort2 +RTS -sstderr radixSort2 +RTS -sstderr 7418991 3,834,816 bytes allocated in the heap 1,103,176 bytes copied during GC 123,072 bytes maximum residency (1 sample(s)) 19,984 bytes maximum slop 2 MB total memory in use (0 MB lost due to fragmentation) Generation 0: 7 collections, 0 parallel, 0.00s, 0.00s elapsed Generation 1: 1 collections, 0 parallel, 0.00s, 0.00s elapsed INIT time 0.00s ( 0.00s elapsed) MUT time 0.00s ( 0.01s elapsed) GC time 0.00s ( 0.00s elapsed) EXIT time 0.00s ( 0.00s elapsed) Total time 0.00s ( 0.01s elapsed) %GC time 0.0% (11.1% elapsed) Alloc rate 3834,816,000,000 bytes per MUT second Productivity 100.0% of total user, 0.0% of total elapsed
 なんとこちらは、私の書いた「クイックソート」や「マージソート」より速くなってしまいました。

 ソートする要素の桁数が増えれば、その分「基数ソート」は時間がかかります。そのためソートする要素を選ぶソート法かもしれませんが、ツボに嵌れば強力なソート法だということが分りました。

2011年2月13日 (日)

Haskell で行列の計算をしてみる

 みなさんは『オイラーの贈物』という本を御存じですか?
 「オイラーの公式」を導くために必要な、代数、幾何、解析の知識がこの一冊で学べるという本です。高校から大学で習う数学の知識があれば、素晴しい数学の世界を垣間見ることができます。いろいろな定義や公式が導き方とともに載っているので、数学的な知識の確認をしたい時に重宝します。

 ここ数日、この本の「ベクトルと行列」の章を読み返していました。
 読んでいて、ふと「Haskell で行列の計算をするためのモジュールを作ってみよう」と思ました。
 Haskell の行列演算モジュールなど、絶対誰かが作っているはずですが、ちょっと頭の体操代わりに自分で作ってみたくなったんです。
 で、できたのがこれです。

{- Matrix.hs -- ベクトル及び行列の計算 coding by Tsumuji -} module Matrix where import Data.List (transpose) -- ベクトルの演算 -- * ベクトルは数値のリストの形で表す。 -- ex : 3次元のベクトル [1,2,3] type Vector a = [a] -- Scalar * Vector (*/) :: Num a => a -> Vector a -> Vector a (*/) n xs = map (* n) xs infixl 7 */ -- Vector + Vector (/+/) :: Num a => Vector a -> Vector a -> Vector a (/+/) xs ys = zipWith (+) xs ys infixl 6 /+/ -- Vector - Vector (/-/) :: Num a => Vector a -> Vector a -> Vector a (/-/) xs ys = zipWith (-) xs ys infixl 6 /-/ -- Vector * Vector (inner product) (/*/) :: Num a => Vector a -> Vector a -> a (/*/) xs ys = sum $ zipWith (*) xs ys infixl 7 /*/ -- Vector の絶対値 vectorAbs :: Floating a => Vector a -> a vectorAbs xs = sqrt (xs /*/ xs) -- n 次元のゼロベクトル zeroVector :: Num a => Int -> Vector a zeroVector n = take n $ repeat 0 -- 行列の演算 -- * 行列は行ベクトルのリストの形で表す。 -- ex : 2行3列の行列 [[1,2,3],[4,5,6]] -- * 行と列の整合性は調べていないので注意!! type Matrix a = [Vector a] -- 行列の m 行目 n 列目の要素(通常のリストと同様に 0 から数え始める) matrixRef :: Num a => Int -> Int -> Matrix a -> a matrixRef m n xss = (xss !! m) !! n -- Scalar * Matrix (*|) :: Num a => a -> Matrix a -> Matrix a (*|) n xss = map (n */) xss infixl 7 *| -- Matrix + Matrix (|+|) :: Num a => Matrix a -> Matrix a -> Matrix a (|+|) xss yss = zipWith (/+/) xss yss infixl 6 |+| -- Matrix - Matrix (|-|) :: Num a => Matrix a -> Matrix a -> Matrix a (|-|) xss yss = zipWith (/-/) xss yss infixl 6 |-| -- Matrix * Matrix (|*|) :: Num a => Matrix a -> Matrix a -> Matrix a (|*|) xss yss = [[xs /*/ ys | ys <- yss'] | xs <- xss] where yss' = transpose yss infixl 7 |*| -- m 行 n 列のゼロ行列 zeroMatrix :: Num a => Int -> Int -> Matrix a zeroMatrix m n = take m $ repeat $ zeroVector n -- n 行 n 列の単位行列 unitMatrix :: Num a => Int -> Matrix a unitMatrix n = [take n $ f x | x <- [0 .. n - 1]] where f x = replicate x 0 ++ [1] ++ repeat 0 -- Matrix ^ n (正方行列のみ可) (|^) :: Num a => Matrix a -> Int -> Matrix a (|^) xss 0 = unitMatrix (length xss) (|^) xss n = power (|*|) xss n infixr 8 |^ -- 高速累乗計算 -- ex : 2^100 == power (*) 2 100 power :: (t -> t -> t) -> t -> Int -> t power op a n = loop (n - 1) a a where loop 0 _ r = r loop n x r = loop (div n 2) (op x x) r' where r' = if odd n then op r x else r
 ベクトルや行列の演算では、要素にランダムにアクセスすることは無いので、データの型をベクトルは「数字のリスト」、行列は「ベクトルのリスト」としました。
 ベクトルや行列のサイズや整合性をチェックしていない手抜きバージョンなので、使うときにはちょっと注意が必要ですが……。

 さらに、『オイラーの贈物』には行列の計算でフィボナッチ数列を求める方法も載っていたので、これもコーディングしてみました。

{- fib_matrix.hs coding by Tsumuji -} import Matrix -- フィボナッチ数列の第n項を求める (n > 0) fibonacci :: Int -> Integer fibonacci n = matrixRef 1 0 $ [[1,1],[1,0]] |^ (n - 1) |*| [[1],[1]] main :: IO () main = print $ map fibonacci [1..20]
>ghc -O -o fib fib_matrix.hs --make ghc -O -o fib fib_matrix.hs --make Linking fib.exe ... >fib fib [1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181,6765]

2011年2月 6日 (日)

Project Euler : Problem 58

 問題はこちらをご覧ください。
 また、自作の "ForEuler module" に関してはこちらをご覧ください。

 

 辺の長さと四隅の数の関係を調べてみると、次のようなことが分ります。

辺の長さ : 四隅の数 : 交差 1 : 1 : 0 3 : 9, 7, 5, 3 : -2 5 : 25, 21, 17, 13 : -4 7 : 49, 43, 37, 31 : -6
 つまり四隅の数は、辺の長さを L とすると「初項 L^2, 公差 1-L の等差数列」と見ることができます。
 ということで次のようなコードを書いてみました。
import ForEuler (isPrime) problem058 :: Integer problem058 = head [s | (s, _, _) <- dropWhile check $ tail ts] where ts = scanl f (1, 0, 1) [3, 5 ..] check (_, p, t) = 10 * p >= t f (_, p, t) s = (s, p + length (pNums s), t + 4) pNums x = filter isPrime [x ^ 2 + (1 - x) * y | y <- [0 .. 3]] main :: IO () main = print problem058
 対角線上の数の個数の合計を t 、そのうちの素数の個数の合計を p とした時、dropWhile に与える条件を「p / t >= 0.1」とすると、計算の度に数値の型を Int から Float に変換しなければいけません。
 そこで条件式を変形して「10 * p >= t」とすることで Int だけで計算できるようにしてみました。

« 2011年1月 | トップページ | 2011年3月 »

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            
フォト

最近のトラックバック

無料ブログはココログ