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

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

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
(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桁まで探させてみたところ約三時間で終わりましたが、結局一つしか見つからなくてちょっと残念です。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
(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))))

Index

Feed

Other

Link

Pathtraq

loading...