参考：什么是梅森数

http://zh.wikipedia.org/w/index.php?title=%E6%A2%85%E6%A3%AE%E7%B4%A0%E6%95%B0&variant=zh-cn

最近看见有人用ruby算梅森数。

鄙人用groovy也来玩玩这个命题。

一：首先，给一个原始的实现：
boolean isPrime(n) { if(n==1) return false; if(n==2) return true; for(i in 2..Math.round(n**(0.5))){if(n%i==0) return false} return true; } mersennes = []; (2..31).each{n=2**(it)-1;if(isPrime(n)) mersennes << n} println mersennes;
 

我的机器耗时2秒算到n=31的梅森数2147483647
,达到欧拉的水平. O^O
之后就很难了， 等了5分钟也没算出n=22的梅森数。

二：继续优化：

算法分析：
f(n)=2**n-1
f(n+1)=2**(n+1)-1=(2**n)*2-1=(2**n)*2-2+1=2(2**n-1)+1=2f(n)+1

==>
f(n+1)=2f(n)+1

eg:
f(1)=2**1-1=1
f(2)=2*f(1)+1=3
f(3)=2f(2)+1=7
f(4)=2f(3)+1=15
...

代码如下，注意：加入了一个cache， 存放已经运算出来的f(n):
boolean isPrime(n) { if(n==1) return false; if(n==2) return true; for(i in 2..Math.round(n**(0.5))){if(n%i==0) return false} return true; } cache = new HashMap(); def function(n) { if(cache.getAt(n)!=null) return cache.getAt(n) value = (n==1?1:2*function(n-1)+1) cache.putAt(n, value); return value; } (2..31).each{v=function(it);if(isPrime(v)) println "f"+it+"="+v}

还是算到f31， 用时200毫秒，速度提高了一个数量级。
f61还是出不来。
大数据在groovy中默认为double的， 看来要换成BigInteger试一下。

三：终极优化

采用BigInteger的方法isProbablePrime来判断是否为素数，速度瞬间上升：

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

15分钟算到M4000，达到1957年的世界纪录，令人鼓舞：
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

发表于 @ 2009年04月20日　10:31:00 | 评论( loading... ) | 编辑| 举报| 收藏

旧一篇:山寨版工作流：groovy控制的责任链 | 新一篇:一个cool Map