教你用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





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

我算到了第15个（偷懒借用了别人的代码，rabin-miller.rb如附件）

require 'rabin-miller.rb'
(1..10000).each do |n|
  if n.prime? # 如果n为合数，2**n-1必为合数
    m = 2**n -1
    puts m if m.prime?
  end
end

补充：改进了下算法，算到第17个了……

再补充：#18
25911708601320262777624676792244153094181888755312542730397492316187401926658636
20862012095168004834065506952417331941774416895092388070174103777095975120423130
66624082916353517952311186154862265604547691127595848775610568757931191017711408
82625215384903583040118507211642474746182303147139834022928807454567790794103728
82358207058923510684338829868886166586502809276920803396058693087905004095037098
75902119018371991620994002568935113136548829739112656797303241986517250116412703
50970542777347797234982167644344666838311932254009964899405179024162405651905448
36908096160616257430423617218633394158524264312087372665919620617535357488928945
99629195183082621860853400937932839420261866586142503251450773096274235376822938
64940712770084607712421182308080413929808705750471382526457144837937112503208182
61265666490842516994539518877896136502484057393785945994443352311882801236604062
62468609212150349937584782292237144339628858485938215738821232393687046160677362
909315071

第一天大概能算20几个……

写了个 Haskell 版本，用 Lucas-Lehmer 方法检验，速度比 Rabin-Miller 更快：
http://night-stalker.iteye.com/blog/363154

昨晚只看到你的第一条回复，然后我就想着今天写个haskell版本的超过你那个ruby的，结果一来公司发现你已经写了个haskell的。。。

恭喜你，你已经被历史记住了^-^

等我想想办法啊，那我用erlang写个吧。

相当强悍的卢卡斯-莱莫算法，我用ruby实现了这个算法之后，发现算到M(19)轻轻松松。
不过算到下一个梅森数的时候就完全不行了，不是算法不行，是ruby不行。

不过我的命运已经可以写进历史了O(∩_∩)O哈哈~

下面是我的代码：

require 'benchmark'

def lucasLehmer_sequence(n)
   n>=0?(n==0?4:lucasLehmer_sequence(n-1)**2-2):(raise ArgumentError,"lucasLehmer seq must larger than 0.")
end

def lucasLehmer_test(n)
  lucasLehmer_sequence(n-2)%mersenne(n)==0
end

def mersenne(n)
  n>=1?(2**n-1):(raise ArgumentError,"Mersenne seq must larger than 0.")
end

stat=nil

(2..70).each do |n|
  rt=Benchmark::realtime{stat=lucasLehmer_test(n)}
  p ("M(#{n}) => #{mersenne(n)} | realtime : #{rt}") if stat
end

输出：

"M(3) => 7 | realtime : 0.0"

"M(5) => 31 | realtime : 0.0"

"M(7) => 127 | realtime : 0.0"

"M(13) => 8191 | realtime : 0.0"

"M(17) => 131071 | realtime : 0.0309998989105225"

"M(19) => 524287 | realtime : 0.562000036239624"

F:/MySummary/ICE/cachePerf/lib/Lucas–Lehmer.rb:4: warning: in a**b, b may be too big

F:/MySummary/ICE/cachePerf/lib/Lucas–Lehmer.rb:4: warning: Bignum out of Float range

只要在生成 Lucas-Lehmer 数列的时候，每轮都 mod 一次梅森数，就能算下去。

另外我用haskell写了个多线程版本，只是本机貌似跑不出100% cpu的效果来……

还有，M(19)不是第7个么…… haskell算到 M(11213)了……

你看我代码，不就是每轮都mod一次梅森数么？

def lucasLehmer_test(n)  
  lucasLehmer_sequence(n-2)%mersenne(n)==0  
end  

ruby的BigNum长度据说只和内存大小有关，看样子不是啊！

我的意思是，将数列的递推公式改成：
lucasLehmer_sequence(n) = (lucasLehmer_sequence(n-1)**2-2) % mersenne(p)

然后 test就变成
lucasLehmer_sequence(p) == 0

这样便可以控制队列中的数的大小。（本质就是是在 N/p 里面做运算）

lucasLehmer_sequence(n)不是卢卡斯莱莫队列，而是队列中的第N个元素，所以不管在那里求余，复杂度都没有下降。

的确是第7个，因为第8个就报那个Bignum的错了，不知道是不是ruby大数计算不行了。

你M(11213)好NB啊！