challenge Meertens数

お題#100「正整数のゲーデル数化?」で定義した goedel を適用すると自分自身になるような数,すなわち goedel (n) == n となるような正整数 n を見つける関数を定義してください.

このような数のことをMeertens数と言うそうです.

32bitsの符号なし整数(あるいは10進10桁整数)までの範囲で探すのにどのくらい計算時間がかかったかをCPUのスペックとともに教えてください.また,その実装で64bit符号なし整数(あるいは10進20桁整数)までの範囲で探すのにどのくらい計算時間がかかりそうか見積ってください(もちろん実際に計算して計算時間を示していただくのでもかまいません).

補足:

自分でやってみたHaskellでは30行ほどのプログラムで、10進10桁の範囲までの探索は ghci で 4.56秒(CPU Intel Core 2 Duo 2.6GHz)でした。

Posted feedbacks - diff

とりあえず ToString 使わないように改造してみた。30秒速くなった。変化がねえ(;´Д`)
 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
--- meertens.old.cs     Wed Dec 19 01:16:26 2007
+++ meertens.cs Wed Dec 19 01:17:20 2007
@@ -18,9 +18,10 @@
     }
     static double Goedel(long num, int[] primes) {
         double goedel = 1;
-        string numstr = num.ToString();
-        for(int i = 0; i < numstr.Length; i++)
-            goedel *= Math.Pow(primes[i], numstr[i] - '0');
+        for(int i = Length(num) - 1; 0 <= i; i--) {
+            goedel *= Math.Pow(primes[i], num % 10);
+            num /= 10;
+        }
         return goedel;
     }
     static int[] Primes(int count) {
@@ -34,5 +35,14 @@
         SKIP:;
         }
         return primes.ToArray();
+    }
+    static int Length(long num) {
+        long num2 = 10;
+        int cnt = 1;
+        while(num2 <= num) {
+            cnt++;
+            num2 *= 10;
+        }
+        return cnt;
     }
 }

面白くはないけど、これでちょっと速くなりました。

SBCL との差は……あまり自信ありませんが、loop の中でクロージャを作って呼び出すオーバーヘッドとか?

1
2
3
4
5
6
7
8
9
       end

       open Int
-      val lst = List.tabulate (10 - start, fn x => x + start)
+      val lst = if start = 0 then [0,1,2,3,4,5,6,7,8,9]
+                else [1,2,3,4,5,6,7,8,9]
     in
       app f lst
     end

独自のpowでさらに枝刈り。20桁での時間は約4時間までになりました。

 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
32
33
34
35
--- 110.sml     2007-12-21 19:45:26.606130399 +0000
+++ 110-3.sml   2007-12-21 20:58:19.679130289 +0000
@@ -13,16 +13,30 @@
   fun loop _ [] _ prod sum _ =
     if prod = sum then print (toString sum ^ "\n") else ()
     | loop (p::ps) (b::bs) bound prod sum start = let
+      fun pow (x, y) = let
+        fun loop r 0 = 1
+          | loop r 1 = r
+          | loop r n = let
+              val z = r * x
+            in
+              if z <= bound then loop z (n - 1)
+              else bound + 1
+            end
+      in
+        loop x y
+      end
+
       fun f x = let
         val pp = pow (p, x)
       in
         if pp <= bound then
-          loop ps bs (bound div pp) (prod * pp) (sum + b * fromInt x) 0
+          loop ps bs (bound div pp) (prod * pp) (sum + b * x) 0
         else ()
       end

       open Int
-      val lst = List.tabulate (10 - start, fn x => x + start)
+      val lst = if start = 0 then [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
+                else [1, 2, 3, 4, 5, 6, 7, 8, 9]
     in
       app f lst
     end

そっか、そこ刈れてなかったんですね。それなら f の繰り返しを再帰呼び出しにするとまだ速くなります。

 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
--- 4926.sml    2007-12-22 15:37:01.000000000 +0900
+++ meertens.sml        2007-12-22 15:36:18.000000000 +0900
@@ -26,11 +26,12 @@
         loop x y
       end

-      fun f x = let
+      fun f [] = ()
+        | f (x::xs) = let
         val pp = pow (p, x)
       in
         if pp <= bound then
-          loop ps bs (bound div pp) (prod * pp) (sum + b * x) 0
+          (loop ps bs (bound div pp) (prod * pp) (sum + b * x) 0; f xs)
         else ()
       end

@@ -38,7 +39,7 @@
       val lst = if start = 0 then [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
                 else [1, 2, 3, 4, 5, 6, 7, 8, 9]
     in
-      app f lst
+      f lst
     end
 in
   app (fn i => loop p (basis i) (IntInf.pow (10, i)) 1 0 1) lst

Index

Feed

Other

Link

Pathtraq

loading...