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

2010年12月

2010年12月30日 (木)

Haskell で「エラトステネスの篩」 その2

 昨日に引き続き、「エラトステネスの篩」です。

 昨日の記事のコードで十分すぎるくらいの速さが出たので、今度は "mapM_" などを使ってコードをコンパクトにまとめる努力をしてみました。(昨日書いた記事のコードが、試行錯誤中に書いたバグ入りのコードだったことに今日気づきました。現在は正しいコードに直してあります)
 元のコードに比べて、コードがかなり短くなった割には、速度はあまり落ちなかったので、自分としては結構満足しています。

import Data.Array.ST import Data.Array.Unboxed import Data.Array.Base (unsafeRead, unsafeWrite) import Control.Monad sieve :: Int -> UArray Int Bool sieve n = runSTUArray $ do arr <- newArray (0, n) True mapM_ (setFalse arr) $ [0, 1] ++ [4, 6 .. n] mapM_ (loop arr) [3, 5 .. floor $ sqrt $ fromIntegral n] return arr where loop arr k = do v <- unsafeRead arr k when v $ mapM_ (setFalse arr) [k * k, k * (k + 2) .. n] setFalse arr k = unsafeWrite arr k False primeList :: Int -> [Int] primeList n = [p | (p, bool) <- assocs $ sieve n, bool]

2010年12月29日 (水)

Haskell で「エラトステネスの篩」

 Haskell では通常の配列は書き換えのコストが高いので、 Ruby で Project Euler を解くときに作った"prime_list"のような、配列をバリバリ書き換えていくような「エラトステネスの篩」は使い物になりませんでした。
 たぶん STArray あたりを使えば、書き換えのコストが低くなるので実現できそうな気はしていましたが、モナドの知識が不足していて手が出せませんでした。

 そんな時、STArray を使った「エラトステネスの篩」を実際にコーディングしている記事を見つけました。
 STArray 以外にも unsafeRead や unsafeWrite なども使われていて、とてつもなく速いコードでした。
 ただ、よく見ると「sieve 関数」のアルゴリズムに若干改良の余地があったので、ちょっと手を加えてみました。

import Control.Monad import Data.Array.ST import Data.Array.Unboxed import Data.Array.Base(unsafeRead, unsafeWrite, unsafeAt) sieve :: Int -> UArray Int Bool sieve n = runSTUArray $ do arr <- newArray (0, n) True unsafeWrite arr 0 False unsafeWrite arr 1 False let end = floor $ sqrt (fromIntegral n) remove_multiples arr 1 4 loop arr end 3 return arr where loop arr end k | k > end = return () | otherwise = do v <- unsafeRead arr k when v $ remove_multiples arr k (k * k) loop arr end (k + 2) remove_multiples arr p k | k > n = return () | otherwise = do unsafeWrite arr k False remove_multiples arr p (k + p + p)

 ほんのちょっとの手直しなので、どこが変わったか分からないかもしれませんね。
 でも、このちょっとの改良でけっこう速くなったんですよ。家のパソコンで一億までの素数の和を求めたとき、改良前が 3.68s だったのが、改良後は 2.53s になりました。

 それから、この「sieve 関数」を使って素数のリストを返す関数も考えてみました。

primeList :: Int -> [Int] primeList n = [p | (p, bool) <- assocs $ sieve n, bool]

 素数の上限を指定しなければならないという制限はありますが、自作の「primes 関数」よりかなり高速です。

※追記 : 12/29 に書いた記事のコードが間違ってたため、 12/30 にコードの部分だ けを書き換えました。

2010年12月 9日 (木)

Project Euler : Problem 57 ~ 漸化式

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

 

 考え方に関してはこちらをご覧ください。

 今回は、漸化式の分子と分母をタプルにした無限リスト "root2" を定義してみました。

import ForEuler (dexToList) root2 :: [(Integer, Integer)] root2 = iterate func (3, 2) where func (a, b) = (a + 2 * b, a + b) problem057 :: Int problem057 = length $ filter check $ take 1000 root2 where check (a, b) = size a > size b size x = length $ dexToList x main :: IO () main = print problem057

« 2010年11月 | トップページ | 2011年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            
フォト

最近のトラックバック

無料ブログはココログ