Chinaunix

标题: Sieve of Eratosthenes 的 Haskell/Scheme 实现 [打印本页]
作者: MMMIX    时间: 2008-09-21 00:31
标题: Sieve of Eratosthenes 的 Haskell/Scheme 实现
Sieve of Eratosthenes 算法描述见 http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes

Haskell 实现如下:

  2. sieve :: [Integer] -> [Integer]
  3. sieve (x:xs) = x : sieve xs'
  4.              where xs' = filter (\y -> y `mod` x /= 0) xs

  6. primes = 2 : sieve [3,5..]
复制代码


加了平方优化后的版本

  2. sieve :: [Integer] -> [Integer]
  3. sieve (x:xs) = x : sieve (xs' ++ xs''')
  4.              where
  5.                 (xs', xs'') = break (>= x*x) xs
  6.                 xs''' = filter (\y -> y `mod` x /= 0) xs''

  8. primes = 2 : sieve [3,5..]
复制代码


实际上，第二个版本有严重的性能问题(break 会导致大量的内存操作)以至于还没有第一个快。

真正优化的版本(by win_hate):

  2. isPrime :: Integer -> [Integer] -> Bool
  3. isPrime p (x:xs)
  4.     | x*x > p           = True
  5.     | p `mod` x == 0    = False
  6.     | otherwise         = isPrime p xs

  8. primes = 2 : [ p | p <- [3,5..], isPrime p primes ]
复制代码


BTW，这个还可以再稍加优化。

[ 本帖最后由 MMMIX 于 2009-5-20 14:15 编辑 ]
作者: win_hate    时间: 2008-09-21 01:20
lazy scheme


  2. (define (integers n)
  3.   (cons n (integers (+ n 1))))

  5. (define (sieve seq)
  6.   (let ((p (car seq)))
  7.   (cons p (sieve
  8.            (filter (lambda (x) (not (= 0 (modulo x p)))) (cdr seq))))))
  9.            
  10. (define primes (sieve (integers 2)))
复制代码

> (list-ref primes 50)
233
作者: MMMIX    时间: 2008-09-21 01:59
标题: 回复 #2 win_hate 的帖子
我觉着 integers 的定义可以优化下
PLT Scheme 实现(需要在 DrScheme 中将语言设为 Module)

  2. #lang lazy

  4. (define (integers n)
  5.   (cons n (integers (+ n 2))))

  7. (define (sieve seq)
  8.   (let ((p (car seq)))
  9.   (cons p (sieve
  10.            (filter (lambda (x) (not (= 0 (modulo x p)))) (cdr seq))))))
  11.            
  12. (define primes (cons 2
  13.                      (sieve (integers 3))))
复制代码


  2. > (list-ref primes 50)
  3. 233
复制代码

[ 本帖最后由 MMMIX 于 2008-9-21 10:51 编辑 ]
作者: win_hate    时间: 2008-09-21 10:15
一个带平方优化的版本


  2. (define odd
  3.   (letrec ((f (lambda (n)
  4.              (cons n (f (+ n 2))))))
  5.     (f 3)))

  7. (define primes
  8.   (cons 2 (filter primes? odd)))

  10. (define (primes? n)
  11.   (define (f ps)
  12.     (let ((p (car ps)))
  13.       (cond ((> (* p p) n) #t)
  14.             ((= (modulo n p) 0) #f)
  15.             (else (f (cdr ps))))))
  16.   (f primes))
复制代码
作者: win_hate    时间: 2008-09-21 10:17
原帖由 MMMIX 于 2008-9-21 01:59 发表
我觉着 integers 的定义可以优化下
PLT Scheme 实现

#lang lazy

(define (integers n)
  (cons n (integers (+ n 2))))

(define (sieve seq)
  (let ((p (car seq)))
  (cons p (sieve
         ...


说得对，用奇数列相当于在选材上先 filter 了一次。

那个 #lang lazy 在 mzscheme 中好用吗？我这里试了不行。

[ 本帖最后由 win_hate 于 2008-9-21 10:34 编辑 ]
作者: MMMIX    时间: 2008-09-21 10:50
原帖由 win_hate 于 2008-9-21 10:17 发表
那个 #lang lazy 在 mzscheme 中好用吗？我这里试了不行。

mzscheme 不清楚(应该也行)，我向来都是直接用 DrScheme 的
在 DrScheme 中应该把语言设为 Module，然后通过 #lang 选择具体要用的语言，这也是推荐的作法。
作者: win_hate    时间: 2008-09-21 11:21
谢谢！


用模块之后可以了，但对整个体系还是很模糊。要找个时间把 PLT 的文档好好看一看才行。


  2. (module nth (lib "lazy.ss" "lazy")
  3.   
  4.   (define odd
  5.     (letrec ((f (lambda (n)
  6.                   (cons n (f (+ n 2))))))
  7.       (f 3)))
  8.   
  9.   (define primes
  10.     (cons 2 (filter primes? odd)))
  11.   
  12.   (define (primes? n)
  13.     (define (f ps)
  14.       (let ((p (car ps)))
  15.         (cond ((> (* p p) n) #t)
  16.               ((= (modulo n p) 0) #f)
  17.               (else (f (cdr ps))))))
  18.     (f primes))
  19.   
  20.   (define (nth-prime n)
  21.     (! (list-ref primes n)))
  22.   
  23.   (provide nth-prime))
复制代码
作者: MMMIX    时间: 2008-09-21 16:21
原帖由 win_hate 于 2008-9-21 11:21 发表
谢谢！


用模块之后可以了，但对整个体系还是很模糊。要找个时间把 PLT 的文档好好看一看才行。


其实用 #lang 也是可以的，不过直接用 (list-ref primes 50) 是没有输出的，但是可以用 (display (list-ref primes 50)) 显示结果。
作者: win_hate    时间: 2008-09-21 17:44
我不是很明白你的意思。

如果直接选 lazy scheme，是可以运行的。我回的第一贴用的就是这个环境



如果把语言选为 module，则得到一条错误。



可能我模块那里没写对？

[ 本帖最后由 win_hate 于 2008-9-21 17:49 编辑 ]
作者: MMMIX    时间: 2008-09-21 18:28
原帖由 win_hate 于 2008-9-21 17:44 发表


如果把语言选为 module，则得到一条错误。


估计是版本的问题，我看你用的还是 3.72, 我测试时用的是 4.1

collects/lazy/lang/reader.ss 这个文件在我这是有的，内容很简单：

  2. (module reader syntax/module-reader
  3.   lazy)
复制代码

[ 本帖最后由 MMMIX 于 2008-9-21 18:31 编辑 ]
作者: win_hate    时间: 2008-09-21 19:10
我下个官方的 svn 版本捣鼓一下。你用的是那个版本吗？

就我现在这个版本，只要在代码的开头加入


  2. (require (lib "lazy.ss" "lazy"))
复制代码


则 mzscheme 就变成 lazy 的了。

[ 本帖最后由 win_hate 于 2008-9-21 19:15 编辑 ]
作者: MMMIX    时间: 2008-09-21 19:23
原帖由 win_hate 于 2008-9-21 19:10 发表
我下个官方的 svn 版本捣鼓一下。你用的是那个版本吗？


我用的是正式发布的 4.1 版, http://download.plt-scheme.org/drscheme/
作者: win_hate    时间: 2008-09-21 19:26
没有我可用的版本(XUbuntu)，x86_64 那个是红帽的，估计不好用。
作者: MMMIX    时间: 2008-09-21 19:33
原帖由 win_hate 于 2008-9-21 19:26 发表
没有我可用的版本(XUbuntu)，x86_64 那个是红帽的，估计不好用。

那估计你要用的话只能自己从源码编译了。
作者: win_hate    时间: 2008-09-21 22:17
svn check out 的版本(11826)，4.1.0.3

现在用 lang 可以用了，谢谢！

[ 本帖最后由 win_hate 于 2008-9-21 22:25 编辑 ]
作者: MMMIX    时间: 2008-09-21 22:51
原帖由 win_hate 于 2008-9-21 22:17 发表
svn check out 的版本(11826)，4.1.0.3

现在用 lang 可以用了，谢谢！

不客气
作者: MMMIX    时间: 2008-09-24 15:30
原帖由 win_hate 于 2008-9-21 10:15 发表
一个带平方优化的版本

这个你有没有做过性能测试？反正我的改进版我测了一下，反而没有改进之前的快 当然，算法本身没有问题，是实现不当造成的。
作者: win_hate    时间: 2008-09-24 16:10
n 比较小的时候用这种方法可能没什么好处。但用大一点的 n，比如 1000， 在我这里测，差别很大。你看看我的截图：

s1.png (44.89 KB, 下载次数: 86)

s1.png

s2.png (42.41 KB, 下载次数: 80)

s2.png
作者: win_hate    时间: 2008-09-24 16:20
学习 Haskell，你的第一个代码我能看明白，第二个就看不懂了。帮我解释一下吧


  2. sieve :: [Integer] -> [Integer]
  3. sieve (x:xs) = x : sieve (xs' ++ xs''')            
  4.              where
  5.                 (xs', xs'') = break (>= x*x) xs      
  6.                 xs''' = filter (\y -> y `mod` x /= 0) xs''

  8. primes = 2 : sieve [3,5..]
复制代码


作者: MMMIX    时间: 2008-09-24 16:43
原帖由 win_hate 于 2008-9-24 16:20 发表
学习 Haskell，你的第一个代码我能看明白，第二个就看不懂了。帮我解释一下吧


xs'++xs'''，中的 ++ 是什么意思？

++ 是 list append 操作，例如 "abc" ++ "def" => "abcdef"

在 Haskell 中，"abc" 是 ['a', 'b', 'c'] 的缩写

(xs', xs'') = break (>= x*x) xs，这里是给 xs' 和 xs'' 绑定值吗？如果是，分别绑定了什么？

是给 xs' 和 xs'' 绑定值(通过 pattern matching)。break p xs 会将 xs 从不满足 p 的那一点分成两部分，分别绑定到 xs' 和 xs''，例如
> break (> 3) [1,2,4,1,5]
([1,2], [4,1,5])
作者: win_hate    时间: 2008-09-24 20:47
有点明白了

按我对第二段代码的理解，设当前素数为 p，你把后面的数分成前后两段。只用 p 来筛后段，而前段让 p 之前的素数来筛。

但是你分段的时候，是通过计算平方来判断的，这与直接拿 p 来模相比，效率不会有明显改进。而且，用 p 来模可以去掉一些合数，计算平方并不能去掉合数，仅仅是把淘汰的工作留给 p 之前的素数了。基于这些分析，第二段代码性能可能不如第一段。
作者: MMMIX    时间: 2008-09-24 21:23
原帖由 win_hate 于 2008-9-24 20:47 发表

按我对第二段代码的理解，设当前素数为 p，你把后面的数分成前后两段。只用 p 来筛后段，而前段让 p 之前的素数来筛。

没错，前半段是由 p 之前的素数来筛的。
但是你分段的时候，是通过计算平方来判断的，这与直接拿 p 来模相比，效率不会有明显改进。而且，用 p 来模可以去掉一些合数，计算平方并不能去掉合数，仅仅是把淘汰的工作留给 p 之前的素数了。基于这些分析，第二段代码性能可能不如第一段。

分段之后，一些取模运算可以变为比较，运算量应该会小些。

我观察到的变慢的情况主要是大量内存操作引起的，实际上，profile 数据显示，80% 甚至更多的时间都是在做垃圾回收(GC).

真正优化后的版本应该是:

  2. sieve :: [Integer] -> [Integer]
  3. sieve (x:xs) = x : sieve xs'
  4.              where
  5.                 xs' = filter isPrime xs
  6.                 isPrime y
  7.                     | y < x*x           = True
  8.                     | y `mod` x /= 0    = True
  9.                     | otherwise         = False

  11. primes = 2 : sieve [3,5..]
复制代码


这是执行结果(编译时使用 -O2 优化选项)
lee@debian:~/tmp/sieve$ time ./sieve0 10000    －－ 未加平方优化的版本
104743

real        0m26.741s
user        0m26.698s
sys        0m0.024s
lee@debian:~/tmp/sieve$ time ./sieve3 10000    －－ 上述加了平方优化的版本
104743

real        0m19.842s
user        0m19.733s
sys        0m0.020s

[ 本帖最后由 MMMIX 于 2008-9-24 21:48 编辑 ]
作者: win_hate    时间: 2008-09-24 22:44
我说的平方优化和你的还不大一样，学着用 Haskell 写了一个(逻辑上跟我前面的 scheme 是一致的)。


  2. isprime :: Integer->[Integer]->Bool
  3. isprime n (x:xs) | x*x > n = True
  4.                  | n `mod` x == 0 = False
  5.                  | otherwise = isprime n xs
  6.                               

  8. primes = 2: filter (\x->isprime x primes) [3..]
复制代码

[ 本帖最后由 win_hate 于 2008-9-24 22:46 编辑 ]
作者: flw    时间: 2008-09-25 09:18
在我机器上 win_hate 这个版本比前面 MMMIX 那个在计算 10000 时快两个数量级。

  1. flw@debian:~/study$ ghc --make -O2 3mix_2.hs
  2. flw@debian:~/study$ time ./3mix_2 10000
  3. 104743

  5. real    0m4.027s
  6. user    0m3.728s
  7. sys     0m0.052s
  8. flw@debian:~/study$ ghc --make -O2 win_hate.hs
  9. flw@debian:~/study$ time ./win_hate 10000
  10. 104743

  12. real    0m0.103s
  13. user    0m0.092s
  14. sys     0m0.008s
  15. flw@debian:~/study$
复制代码


[ 本帖最后由 flw 于 2008-9-25 09:20 编辑 ]
作者: MMMIX    时间: 2008-09-25 09:31
原帖由 win_hate 于 2008-9-24 22:44 发表
我说的平方优化和你的还不大一样，学着用 Haskell 写了一个(逻辑上跟我前面的 scheme 是一致的)。

明白了，这个确实可以提前过滤掉许多数字。
作者: MMMIX    时间: 2008-09-25 09:39
原帖由 flw 于 2008-9-25 09:18 发表
在我机器上 win_hate 这个版本比前面 MMMIX 那个在计算 10000 时快两个数量级。

嗯，我那个实现在数字增大的时候性能会急剧下降，占用的内存太多了。
作者: MMMIX    时间: 2008-09-25 19:37
原帖由 win_hate 于 2008-9-24 22:44 发表
我说的平方优化和你的还不大一样，学着用 Haskell 写了一个(逻辑上跟我前面的 scheme 是一致的)。

帮你润饰下

  2. isPrime :: Integer -> [Integer] -> Bool
  3. isPrime p (x:xs)
  4.     | x*x > p           = True
  5.     | p `mod` x == 0    = False
  6.     | otherwise         = isPrime p xs

  8. primes = 2 : [ p | p <- [3,5..], isPrime p primes ]
复制代码
作者: win_hate    时间: 2008-09-25 21:36
恩， list comprehension 的确是很漂亮的表达方式。
作者: MMMIX    时间: 2009-05-21 10:29
HackageDB 上前不久上传了一个叫做 primes 的 package, 下面是其描述：
primes: Efficient, purely functional generation of prime numbers
This Haskell library provides an efficient lazy wheel sieve for prime generation inspired by Lazy wheel sieves and spirals of primes by Colin Runciman and The Genuine Sieve of Eratosthenes by Melissa O'Neil.


http://hackage.haskell.org/cgi-bin/hackage-scripts/package/primes



欢迎光临 Chinaunix (http://bbs.chinaunix.net/) Powered by Discuz! X3.2