DouKaku? 4865 kozima: (meertens n) で n 桁以下...(Meertens数)

Comment detail

(meertens n) で n 桁以下のものを探します。n=10 のとき、CLISP で1.6秒、SBCL で 0.5秒ほど(Celeron 2.66GHz)。

8～13 で試してみたところ桁数を1増やすごとにほぼ3倍になっているようなので、20桁の場合 SBCL なら単純計算で8時間ぐらいでしょうか。

(defconstant *primes*
  '(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97))

(defun primes (n) (loop repeat n for p in *primes* collect p))
(defun basis (n) (loop repeat n for x = (expt 10 (1- n)) then (/ x 10) collect x))

(defun f (primes basis bound prod sum zero-ok)
  (if primes
      (loop for i from (if zero-ok 0 1) to 9
        as dp = (expt (car primes) i) and ds = (* i (car basis))
        while (<= dp bound) nconc
        (f (cdr primes) (cdr basis) (floor bound dp)
           (* prod dp) (+ sum ds)
           t))
    (if (= prod sum) (list sum) nil)))

(defun meertens (n)
  (loop for i from 1 to n nconc (f (primes i) (basis i) (expt 10 i) 1 0 nil)))

少し効率化して時間が 2/3 程度になりました。寝てる間に20桁まで探させてみたところ約三時間で終わりましたが、結局一つしか見つからなくてちょっと残念です。

(defun f (primes basis bound prod sum zero-ok)
  (declare (optimize (speed 3) (safety 0))
           (type integer bound prod sum))
  (loop repeat (if zero-ok 10 9)
    with p fixnum = (car primes) and prest = (cdr primes)
    and b integer = (car basis) and brest = (cdr basis)
    for newbd integer = (if zero-ok bound (floor bound p)) then (floor newbd p)
    for newprod integer = (if zero-ok prod (* prod p)) then (* newprod p)
    for newsum integer = (if zero-ok sum (+ sum b)) then (+ newsum b)
    while (and (plusp newbd) (<= newprod newsum)) nconc
    (if prest
        (f prest brest newbd newprod newsum t)
      (if (= newprod newsum) (list newsum) nil))))

Haskell の勉強にと移植してみました。GHCi で33秒とあんまり速くはありませんが、わりとすっきりしたコードになったので投稿してみます。ついでに中身について少し説明でも。

上の桁から順に決めていくのですが、その際に桁数を固定することでゲーデル数化する前と後の数(sum, prod)の両方を部分的に計算でき、枝刈りがやりやすくなります。

実際にやっているのは次の二つですが、他にも末尾の 0 の数で刈れたりするかもしれません。

積>和になったら打ち切り(積のほうが速く増えるので)
bound にまだ決まっていない部分の積の上限を保持し、それを超えたら打ち切り

後者のは実質「積が 10^i を超えたら」なのですが、上のようにした方が速いようです。(扱う数が fixnum に収まりやすい？)

primes = flip take [2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
                     31, 37, 41, 43, 47, 53, 59, 61, 67, 71]
basis n = map (\x -> 10^x) [n-1,n-2..0]

f :: [Integer] -> [Integer] -> Integer -> Integer -> Integer -> Integer -> [Integer]
f [] [] _ prod sum _
    | prod == sum = [sum]
    | otherwise   = []
f (p:ps) (b:bs) bound prod sum start
    | start == 0 && prod > sum = []
    | otherwise =
        [x | i <- [start..9], p^i <= bound,
             x <- f ps bs (div bound $ p^i) (prod * p^i) (sum + b * i) 0]

meertens n = [x | i <- [1..n], x <- f (primes i) (basis i) (10^i) 1 0 1]

あーそうか。桁数を固定すればよかったのか。
そうすればかなりの量が枝狩りできますね。

桁数可変で探索してたから、枝がまったく狩れなかった。

ちょっと書き直してみよう。

SMLにも移植。爆速になったけど、MLtonでもSBCLより遅い。。。
10桁で0.79秒、20桁で9時間56分。

  ./110 10  0.79s user 0.00s system 99% cpu 0.795 total
  ./110 20  35746.36s user 4.25s system 99% cpu 9:56:05.18 total

val p = [
   2,  3,  5,  7, 11, 13, 17, 19, 23, 29,
  31, 37, 41, 43, 47, 53, 59, 61, 67, 71
  ] : IntInf.int list

fun basis n = (rev o map (fn x => IntInf.pow (10, x)) o List.tabulate) (n, fn x => x)

fun meertens n = let
  val lst = List.tabulate (n, fn x => x + 1)

  open IntInf

  fun loop _ [] _ prod sum _ =
    if prod = sum then print (toString sum ^ "\n") else ()
    | loop (p::ps) (b::bs) bound prod sum start = let
      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
        else ()
      end

      open Int
      val lst = List.tabulate (10 - start, fn x => x + start)
    in
      app f lst
    end
in
  app (fn i => loop p (basis i) (IntInf.pow (10, i)) 1 0 1) lst
end

val _ = (meertens o valOf o Int.fromString o hd o CommandLine.arguments) ()

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

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

       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時間までになりました。

--- 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 の繰り返しを再帰呼び出しにするとまだ速くなります。

--- 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

DouKaku?

Comment detail

Index

Feed

Other

Link

Pathtraq