教你用Ruby算命！

本文又名《看看我的破机器能算多少个梅森数出来》

代码如下，

mersennes=[]

def is_prime?(n)
  # 这里是用了费马小定理，很慢很慢！
  (2..n-1).each{|y| y**(n-1)%n==1?true:(return false)}
end

(1..13).each do |n|
  m=(2**n-1)
  mersennes<<m if is_prime?(m)
end

#mersennes=(1..13).inject([]){|arr,x|is_prime?(2**x-1)?arr:(arr<<2**x-1)}
#上面这个写法，2**x要计算两次，写法好看，但性能很低

p mersennes

我算到(1..14)这里就不敢再增加range的范围了，(1..15)我的台式机跑了半天也没有跑出结果，看来我的电脑太慢了。

最后我得到的梅森数数组是：

[1, 3, 7, 31, 127, 8191]

我只算到了M4，1不算梅森数。

对照梅森数的发现史，我发现历史上从M5开始的发现者才被历史记录，M4及以前的发现者都被遗忘了，或者估计其它的成就没有多少，碌碌无为。以前算到M4的那哥们是个无名氏，苦命的孩儿啊。我估计大部分人也都差不多会碌碌无为终了一生。目前我只能算到这一位，等有空用Erlang改进下看看，看看自己能否进入被历史记住的人物名单。

我说完了，你们自己也算算吧，看你们排行老几。

PS：我那个is_prime?这个判断素数的算法是很浪费性能的，不过它不会错过一个漏网之鱼，比Miller－Rabin算法好在这里，但是太慢！

梅森素数列表：http://zh.wikipedia.org/w/index.php?title=梅森素数&variant=zh-cn





《纯粹娱乐，别太计较，不过很好玩^-^》

评论

61 楼 oCameLo 2009-04-21  

-_-b

浪费电。。

60 楼 CharlesCui 2009-04-21  

Erlang一个晚上的结果,从M(2)到M(9941)

M(9941) => 34608828249085121524296039576741331672262866890023854779048928344500622080983411446436437554415370753366448674763505018641470709332373970608376690404229265789647993709760358469552319045484910050304149809818540283507159683562232941968059762281334544739720849260904855192770626054911793590389060795981163838721432994278763633095377438194844866471124967685798888172212033000821469684464956146997194126921284336206463313859537577200462442029064681326087558257488470489384243989270236884978643063093004422939603370010546595386302009073043944482202559097406700597330570799507832963130938739885080198416258635194522913042562936679859587495721031173747796418895060701941717506001937152430032363631934265798516236047451209089864707430780362298307038193445486493756647991804258775574973833903315735082891029392359352758617185019942554834671861074548772439880729606244911940066680112823824095816458261761861746604034802056466823143718255492784779380991749580255263323326536457743894150848953969902818530057870876229329803338285735419228259022169602665532210834789602051686546011466737981306056247480055071718250333737502267307344178512950738594330684340802698228963986562732597175372087295649072830289749771358330867951508710859216743218522918811670637448496498549094430541277444079407989539857469452772132166580885754360477408842913327292948696897496141614919739845432835894324473601387609643750514699215032683744527071718684091832170948369396280061184593746143589068811190253101873595319156107319196071150598488070027088705842749605203063194191166922106176157609367241948160625989032127984748081075324382632093913796444665700601391278360323002267434295194325607280661260119378719405151497555187549252134264394645963853964913309697776533329401822158003182889278072368602128982710306618115118964131893657845400296860012420391376964670183983594954112484565597312460737798777092071706710824503707457220155015899591766244957768006802482976673920392995410164224776445671222149803657927708412925555542817045572430846389988129960519227313987291200902060882060733762075892299473666405897427035811786879875694315078654420055603469625309399653955932310466430039146465805452965014040019423897552675534768248624631951431493188170905972588780111850281190559073677771187432814088678674286302108275149258477101296451833651979717375170900505673645964696355331369819296000267389583289299126738345726980325998955997501176664201042888546085699446442834195232948787488410595750197438786353119204210855804692460582533832967771946911459901921324984968810021189968284941331573164056304725480868921823442538199590383852412786840833479611419970101792978355653650755329138298654246225346827207503606740745956958127383748717825918527473164970582095181312905519242710280573023145554793628499010509296055849712377978984921839997037415897674154830708629145484724536724572622450131479992681684310464449439022250504859250834761894788889552527898400988196200014868575640233136509145628127191354858275083907891469979019426224883789463551

59 楼 CharlesCui 2009-04-20  

同意。

目前的算法大概都一样的吧，都是用卢卡斯莱莫的那个方法，而且都只用一台计算机做计算。

我们该想想怎么把这个东东写成多个计算机分布计算才NB。

58 楼 night_stalker 2009-04-20  

这样的比较没什么意思……

对大规模的问题，语言的差别撑死也就常数级，算法复杂度才是性能的主宰。
假设 Intel C 平均操作比 Ruby 快 100 倍，那么同样的算法就差 100 倍左右。

看看算法复杂度的比较例子：

时间复杂度 O(n^5) 的算法，是 O(n^4) 耗时的 n 倍。
如果 n 是 100 位整数，那么区别就是：
1亿亿亿亿亿亿亿亿亿亿亿亿亿亿亿亿亿亿亿亿亿亿亿亿亿 倍

57 楼 CharlesCui 2009-04-20  

Erlang版正好快到15分钟也计算到了M(3217),

看来对于这段算法，Erlang和Groovy这种运行在虚拟机上的编译型语言速度都差不多啊。

我那个机器是个虚拟机，2核2G。

56 楼 CharlesCui 2009-04-20  

lcllcl987 写道

15分钟算到M4000:
f(2) = 3
f(3) = 7
f(5) = 31
f(7) = 127
f(13) = 8191
f(17) = 131071
f(19) = 524287
f(31) = 2147483647
f(61) = 2305843009213693951
f(89) = 618970019642690137449562111
f(107) = 162259276829213363391578010288127
f(127) = 170141183460469231731687303715884105727
f(521) = 6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151
f(607) = 531137992816767098689588206552468627329593117727031923199444138200403559860852242739162502265229285668889329486246501015346579337652707239409519978766587351943831270835393219031728127
f(1279) = 10407932194664399081925240327364085538615262247266704805319112350403608059673360298012239441732324184842421613954281007791383566248323464908139906605677320762924129509389220345773183349661583550472959420547689811211693677147548478866962501384438260291732348885311160828538416585028255604666224831890918801847068222203140521026698435488732958028878050869736186900714720710555703168729087
f(2203) = 1475979915214180235084898622737381736312066145333169775147771216478570297878078949377407337049389289382748507531496480477281264838760259191814463365330269540496961201113430156902396093989090226259326935025281409614983499388222831448598601834318536230923772641390209490231836446899608210795482963763094236630945410832793769905399982457186322944729636418890623372171723742105636440368218459649632948538696905872650486914434637457507280441823676813517852099348660847172579408422316678097670224011990280170474894487426924742108823536808485072502240519452587542875349976558572670229633962575212637477897785501552646522609988869914013540483809865681250419497686697771007
f(2281) = 446087557183758429571151706402101809886208632412859901111991219963404685792820473369112545269003989026153245931124316702395758705693679364790903497461147071065254193353938124978226307947312410798874869040070279328428810311754844108094878252494866760969586998128982645877596028979171536962503068429617331702184750324583009171832104916050157628886606372145501702225925125224076829605427173573964812995250569412480720738476855293681666712844831190877620606786663862190240118570736831901886479225810414714078935386562497968178729127629594924411960961386713946279899275006954917139758796061223803393537381034666494402951052059047968693255388647930440925104186817009640171764133172418132836351
f(3217) = 259117086013202627776246767922441530941818887553125427303974923161874019266586362086201209516800483406550695241733194177441689509238807017410377709597512042313066624082916353517952311186154862265604547691127595848775610568757931191017711408826252153849035830401185072116424747461823031471398340229288074545677907941037288235820705892351068433882986888616658650280927692080339605869308790500409503709875902119018371991620994002568935113136548829739112656797303241986517250116412703509705427773477972349821676443446668383119322540099648994051790241624056519054483690809616061625743042361721863339415852426431208737266591962061753535748892894599629195183082621860853400937932839420261866586142503251450773096274235376822938649407127700846077124211823080804139298087057504713825264571448379371125032081826126566649084251699453951887789613650248405739378594599444335231188280123660406262468609212150349937584782292237144339628858485938215738821232393687046160677362909315071

别的语言算到M(607)也是瞬间的，我跑下Erlang看看到15分钟算到哪里，再贴给你。

55 楼 lcllcl987 2009-04-20  

15分钟算到M4000:
f(2) = 3
f(3) = 7
f(5) = 31
f(7) = 127
f(13) = 8191
f(17) = 131071
f(19) = 524287
f(31) = 2147483647
f(61) = 2305843009213693951
f(89) = 618970019642690137449562111
f(107) = 162259276829213363391578010288127
f(127) = 170141183460469231731687303715884105727
f(521) = 6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151
f(607) = 531137992816767098689588206552468627329593117727031923199444138200403559860852242739162502265229285668889329486246501015346579337652707239409519978766587351943831270835393219031728127
f(1279) = 10407932194664399081925240327364085538615262247266704805319112350403608059673360298012239441732324184842421613954281007791383566248323464908139906605677320762924129509389220345773183349661583550472959420547689811211693677147548478866962501384438260291732348885311160828538416585028255604666224831890918801847068222203140521026698435488732958028878050869736186900714720710555703168729087
f(2203) = 1475979915214180235084898622737381736312066145333169775147771216478570297878078949377407337049389289382748507531496480477281264838760259191814463365330269540496961201113430156902396093989090226259326935025281409614983499388222831448598601834318536230923772641390209490231836446899608210795482963763094236630945410832793769905399982457186322944729636418890623372171723742105636440368218459649632948538696905872650486914434637457507280441823676813517852099348660847172579408422316678097670224011990280170474894487426924742108823536808485072502240519452587542875349976558572670229633962575212637477897785501552646522609988869914013540483809865681250419497686697771007
f(2281) = 446087557183758429571151706402101809886208632412859901111991219963404685792820473369112545269003989026153245931124316702395758705693679364790903497461147071065254193353938124978226307947312410798874869040070279328428810311754844108094878252494866760969586998128982645877596028979171536962503068429617331702184750324583009171832104916050157628886606372145501702225925125224076829605427173573964812995250569412480720738476855293681666712844831190877620606786663862190240118570736831901886479225810414714078935386562497968178729127629594924411960961386713946279899275006954917139758796061223803393537381034666494402951052059047968693255388647930440925104186817009640171764133172418132836351
f(3217) = 259117086013202627776246767922441530941818887553125427303974923161874019266586362086201209516800483406550695241733194177441689509238807017410377709597512042313066624082916353517952311186154862265604547691127595848775610568757931191017711408826252153849035830401185072116424747461823031471398340229288074545677907941037288235820705892351068433882986888616658650280927692080339605869308790500409503709875902119018371991620994002568935113136548829739112656797303241986517250116412703509705427773477972349821676443446668383119322540099648994051790241624056519054483690809616061625743042361721863339415852426431208737266591962061753535748892894599629195183082621860853400937932839420261866586142503251450773096274235376822938649407127700846077124211823080804139298087057504713825264571448379371125032081826126566649084251699453951887789613650248405739378594599444335231188280123660406262468609212150349937584782292237144339628858485938215738821232393687046160677362909315071

54 楼 lcllcl987 2009-04-20  

算了， 上一个groovy版本的，瞬间算算到M1000:

def f = 1g 
(2g..1000g).each { n -> 
f = 2 * f + 1 
if (f.isProbablePrime(100)) 
println "f($n) = $f" 
} 

输出结果如下：
f(2) = 3
f(3) = 7
f(5) = 31
f(7) = 127
f(13) = 8191
f(17) = 131071
f(19) = 524287
f(31) = 2147483647
f(61) = 2305843009213693951
f(89) = 618970019642690137449562111
f(107) = 162259276829213363391578010288127
f(127) = 170141183460469231731687303715884105727
f(521) = 6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151
f(607) = 531137992816767098689588206552468627329593117727031923199444138200403559860852242739162502265229285668889329486246501015346579337652707239409519978766587351943831270835393219031728127

53 楼 ray_linn 2009-04-11  

night_stalker 写道

唔，由于延迟求值特性，虽然尾递归了，但是还有一个巨大的 data thunk ..

需要用 $! 强制求值。

不会爆栈的双线程版本如附件。

+RTS -N2 似乎也没有改善，还是接近50% 60%

-threaded 能让haskell用native thread？

52 楼 night_stalker 2009-04-10  

ray_linn 写道

我觉得Erlang和C# Plinq比并没有什么优势.... haskell我看不太懂，能确认算法是完全一样的么？

晚上找一个3.5的ilasm改一个带tail的C#版本来。

算法都一样，各人线程数不同。

linq 语法其实挺 functional 的，F# 就更像 haskell 了。

并发数多，并且涉及大量共享数据才能体现出 Erlang 优于经典线程模型的特点。但这里不是那种服务器处理 n 用户连接的情况…… 为了充分利用多核 cpu，两个线程保持 busy 就够了，线程太多反而可能会因为内存分页问题搞慢速度。

51 楼 CharlesCui 2009-04-10  

我也觉得Erlang和C#甚至Java比优势并不明显，那不应该叫优势，应该叫特色而已。

不能笼统的说面向函数编程比面向过程编程好，只能说个别地方有特色而已。

就如同Erlang的轻量级进程、高并发就是特色，但其他的和Java、C#比，没觉得有什么牛的。

50 楼 ray_linn 2009-04-10  

我觉得Erlang和C# Plinq比并没有什么优势.... haskell我看不太懂，能确认算法是完全一样的么？

晚上找一个3.5的ilasm改一个带tail的C#版本来。

49 楼 CharlesCui 2009-04-10  

楼上的同学们，看下我修改好的Erlang版的计算Mersenne(LucasLehmer)第二版
跟第一版比较，重写了math:pow/2方法，这个方法可以支持很大的幂运算。


http://charlescui.iteye.com/admin/blogs/365159

-module(mersenne).  
-export([run/1,mersenne/1,lucasLehmer_sequence/2,pow/2]).  
  
-spec pow(integer(), non_neg_integer()) -> integer()  
      ; (float(), non_neg_integer()) -> float().  
  
pow(X, N) when is_integer(N), N >= 0 -> pow(X, N, 1).  
  
pow(_, 0, P) -> P;  
pow(X, N, A) -> pow(X, N-1, A*X).  
  
lucasLehmer_sequence (0, M) -> 4 rem M;  
lucasLehmer_sequence (N, M) when N>0 ->  
        erlang:round((pow(lucasLehmer_sequence(N-1,M),2)-2)) rem M.  
  
mersenne(N) when N>0 ->  
   erlang:round(pow(2,N))-1.  
  
lucasLehmer_test(N) ->  
   lucasLehmer_sequence(N - 2, mersenne(N)) =:=0.  
  
run(2) ->  
    case lucasLehmer_test(2) of  
        true ->  
                io:format("M(~w) => ~w~n", [2,mersenne(2)]);  
        false ->  
                false  
    end;  
run(N) ->  
    case lucasLehmer_test(N) of  
        true ->  
                io:format("M(~w) => ~w~n", [N,mersenne(N)]);  
        false ->  
                false  
    end,  
    run(N-1).  

计算梅森数可以算好多，只是效率还有待提高，因为算到M(2281)还是花了一点时间的，不知道Haskell和C#算到M(2281)花了多少时间，ps：这个Erlang版本还是单进程，只用到CPU的一个核。 有空再修改下，让Erlang并行运算的优势发挥下。

48 楼 night_stalker 2009-04-10  

BigInteger 写着太麻烦，给伪码(和迭代很接近了)：

//tail call sequence
int tcs(int n, int m, int result) {
  if(n == 0){
    return result;
  }
  return tcs((n - 1), m, ((result ^ 2 - 2 ) % m));
}

//L-L sequence
int seq(int n, int m) {
  return tcs(n, m, 4 % m);
}

吃内存也不少，算M(20996011)的时候 280 多M，由于爱惜本本，没算出来就给我停了……

47 楼 ray_linn 2009-04-10  

night_stalker 写道

唔，由于延迟求值特性，虽然尾递归了，但是还有一个巨大的 data thunk ..

需要用 $! 强制求值。

不会爆栈的双线程版本如附件。

我看不太懂haskell，你的尾递归写成java是如何的？

这个版本消耗多少内存？？

46 楼 night_stalker 2009-04-10  

唔，由于延迟求值特性，虽然尾递归了，但是还有一个巨大的 data thunk ..

需要用 $! 强制求值。

不会爆栈的双线程版本如附件。

45 楼 ray_linn 2009-04-10  

night_stalker 写道

ray_linn 写道

没有 +RTS -K来制定栈空间的话，hasksell也崩溃了.

很奇怪它的CPU利用率很低，才50%，没办法全力进行计算...如果算法完全一样的话，hasksell应该会比较慢。。。


hasksell好像还没算完第一个数。。。。大数的时候，都是龟速。

编译使用 -O2 和 -threaded 参数(-O3 似乎更牛？)

ghc -O2 -threaded --make messen.hs

双核的使用2个本地线程运行可以吃到 100%：

messen +RTS -N2

我用的是  -O3 ,不过运算速度真的快不起来，这还是个中等长度的数。。

44 楼 night_stalker 2009-04-10  

ray_linn 写道

没有 +RTS -K来制定栈空间的话，hasksell也崩溃了.

很奇怪它的CPU利用率很低，才50%，没办法全力进行计算...如果算法完全一样的话，hasksell应该会比较慢。。。


hasksell好像还没算完第一个数。。。。大数的时候，都是龟速。

编译使用 -O2 和 -threaded 参数(-O3 似乎更牛？)

ghc -O2 -threaded --make messen.hs

双核的使用2个本地线程运行可以吃到 100%：

messen +RTS -N2

43 楼 CharlesCui 2009-04-10  

night_stalker 写道

Ruby 没有尾递归优化的，何况递归的位置在数列那里……

tail call 版本：

ll_seq' 0 m res = res
ll_seq' n m res = ll_seq' (n - 1) m ((res ^ 2 - 2) `mod` m)

ll_seq  n m = ll_seq' n m (4 `mod` m)

你上面给ruby那个版本做的优化叫什么？
这个方法很好，直接减少了堆栈的长度，否则原始版本还会一个函数一个函数的往里面扔。

42 楼 ray_linn 2009-04-10  

没有 +RTS -K来制定栈空间的话，hasksell也崩溃了.

很奇怪它的CPU利用率很低，才50%，没办法全力进行计算...如果算法完全一样的话，hasksell应该会比较慢。。。


hasksell好像还没算完第一个数。。。。大数的时候，都是龟速。