Project Euler Problem 30

問題30: 2以上の自然数について、各桁を5乗した数の和が元の数と一致する数をすべて求め、その和を答えよ。[=>問題文]
　上限が与えられない問題です。左辺の値いくつまで調べればよいのか。注目するのは「桁数」です。
　k 桁の数を考えます。各桁の最大値は 9 なので右辺 a のとりうる最大値 a-max は次のようになります。

[右辺 a 最大値]
a-max(1): 9^5
a-max(2): 9^5 + 9^5
a-max(3): 9^5 + 9^5 + 9^5
:
a-max(k-1): 9^5 + 9^5 + 9^5 ... 9^5
a-max(k)  : 9^5 + 9^5 + 9^5 ... 9^5 + 9^5

　つまり、a の最大値は桁が増える毎に 9^5 ずつ増えていきます。次は左辺 b の最小値 b-min について考えます。

[左辺 b 最小値]
b-min(1): 1 (0)
b-min(2): 10
b-min(3): 100
:
b-min(k-1): 1000...00  (k-1桁)
b-min(k)  : 1000...000 (k桁)

　k=2 の時は2桁の最小値なので10、k=3の bの最小値は 100 です。b の最小値は桁が増えるごとに 10 倍になります。1桁の最小値は 0 ですが便宜的にここでは 1 とします。
　もう一度 b-min と a-max を並べて書いてみます。

k b-min  a-max
1    1    59049 (9^5)
2   10   118098 (9^5 + 9^5)
3  100   177147 (9^5 + 9^5 + 9^5)
4 1000   236196 (9^5 + 9^5 + 9^5 + 9^5)
:

　このように、はじめのうちは b-min(左辺最小）< a-max (右辺最大)ですが、桁が増える毎の値の増加量は b-min が指数的増加量（というか指数関数そのものですが）であるのに対して、a-max は単調増加です。だから、k を増やしていけばいずれどこかで b-min > a-max となる瞬間があるはずです。そして大小が逆転した後はその差は開いていく一方です。この問題は b = a となる数を探すわけですから、b-min > a-max となった後は調べる必要がありません。
　b-min > a-max となった瞬間の a-max が、調べる値の上限となります。

問題30の解答その1

;;Emacs Lisp
(require 'cl)

(defun problem030a (n)
  (let
      ;; 0^n 〜 9^n を格納した配列。
      ;; n乗計算を省略するためにメモ化しておく。
      ((expt-n-vect
        (do ((v (make-vector 10 0))
             (i 1 (1+ i)))
            ((= i 10) v)
          (aset v i (expt i n))))
       ;; 計算対象の上限。
       ;; k 桁での左辺最小値を bk。 右辺最大値を ak とし、
       ;; bk が ak を超えたとき、その時点の ak の値を
       ;; 計算すべき上限値とする。 
       (limit (do ((k 1 (1+ k))        ; 実は不要だけど、一応
                   (a-max (expt 9 n))  ; a(k)の最大値
                   (b-min 1)           ; b(k)の最小値
                   (a-inc (expt 9 n))) ; a増加量
                  ((<= a-max b-min) a-max)
                (incf a-max a-inc)
                (incf b-min (* b-min 10)))))
    (let*
        ;; 桁分割リスト作成関数。 123 => (1 2 3)
        ((split-digit
          (lambda (x &optional lis)
            (if (zerop x)
                lis
              (funcall split-digit (/ x 10)
                       (cons (mod x 10) lis)))))
         ;; 左辺 b から 右辺 a を計算する関数
         (calc-a
          (lambda (b)
            (reduce
             #'(lambda (acc y)
                 (+ acc (aref expt-n-vect y)))
             (funcall split-digit b)
             :initial-value 0))))
      ;; メインルーチン
      ;; 左辺 b (limit 〜 2) に対する 右辺 a を計算し
      ;; b = a となる b をリストに保存し、合計する。
      (do ((lis) (b limit (1- b)))
          ((= b 1) (apply #'+ lis))
        (when (= b (funcall calc-a b))
          (push b lis))))))

(problem030a 5)

　やってることは、最初の説明のとおりです。最後のメインルーチンでは、左辺 b の値を上限 limit から 2 までデクリメントし、すべての b に対する a を計算して b = a となるものを得ています。
　このやり方だとn=5の場合、約47万回の右辺計算をしています。しかもその中には b = 1234,1342,2341... のように右辺 a の値が同じになるものが多数含まれています。この無駄をなんとかしてみましょう。

問題30の解答その2

(require 'cl)

(defun problem030b (n)
  ;; --- 前半は、problem030a と同じ ----------
  (let ((expt-n-vect                    ; 0^n 〜 9^n を格納した配列
         (do ((v (make-vector 10 0))
              (i 1 (1+ i)))
             ((= i 10) v)
           (aset v i (expt i n))))
       ;; 計算対象の上限。
       (limit (do ((k 1 (1+ k))        ; 実は不要だけど、一応
                   (a-max (expt 9 n))  ; a(k)の最大値
                   (b-min 1)           ; b(k)の最小値
                   (a-inc (expt 9 n))) ; a増加量
                  ((<= a-max b-min) a-max)
                (incf a-max a-inc)
                (incf b-min (* b-min 10)))))
    (let
        ;; 桁分割リスト作成関数。 123 => (1 2 3)
        ((split-digit
          (lambda (x &optional lis)
            (if (zerop x)
                lis
              (funcall split-digit (/ x 10)
                       (cons (mod x 10) lis)))))
         ;; --- ここまで、problem030a と同じ ----------
         ;; リスト要素のn乗和計算関数。 (1 2 3) => 1^n+2^n+3^n
         (sum-of-expt-n
          (lambda (lis) (reduce #'(lambda (acc y)
                                    (+ acc (aref expt-n-vect y)))
                                lis
                                :initial-value 0))))
      (let*
          ;; limitの桁数
          ((limit-digits (1+ (truncate (log limit 10))))
           ;; 2〜limit までの数を桁分割し sortしたリストを key として
           ;; 保存した hash tableを作成。hash table の key なので重複は
           ;; 除去される。
           (ht (do ((ht (make-hash-table :test #'equal))
                    (x limit (1- x))
                    (border (expt 10 (1- limit-digits)))
                    (zero-list))
                   ((= x 1) ht)
                 ;; すべてのリストが同じ長さになるよう zero-list を
                 ;; append して調整する。
                 (let ((lis (append zero-list
                                    (funcall split-digit x))))
                   (puthash (sort lis #'<) t ht))
                 ;; x の桁が減ったら、zero-list に 0 を追加
                 (when (= x border)
                   (push 0 zero-list)
                   (setq border (/ border 10)))))
           ;; 題意に該当する数を保存する為のリスト
           (correspond-numbers))
        ;; 作った hash table (ht) のすべてのキーに対して題意のテスト
        ;; を行い合致するものを correspond-numbers に保存する
        (maphash
         #'(lambda (a-list v)
             ;; a-list は右辺の基数となる数値のリスト。
             ;; a-list の要素の n乗和を 左辺 bとし、
             ;; bを桁にばらして sortしたリストを b-list とする.
             (let* ((b (funcall sum-of-expt-n a-list))
                    (b-list (sort (funcall split-digit b) #'<))
                    (b-list-len (length b-list)))
               ;; a-list の桁数はすべて limit-digits なので
               ;; b-list の桁数も 0 埋めしてあわせる
               (if (< b-list-len limit-digits)
                   (setq b-list
                         (append
                          (make-list (- limit-digits b-list-len) 0)
                          b-list)))
               ;; a-list と b-list が一致したら、それが求める数。
               (if (equal a-list b-list)
                   (push b correspond-numbers))))
         ht)
        ;; correspond-numbers には 1 も含まれるので、
        ;; sort して car を除いてから合算。
        (apply #'+ (cdr (sort correspond-numbers #'<))))))))

(problem030b 5)

　1234,1342,2341...のような重複要素を除去してからメインの計算ループに入ろうという発想です。上限 limit までの数をすべて (1 2 3 4) (1 3 4 2) (2 3 4 1) のようなリストに変換し、sort した上で hash table に key として保存しています。ソートすることにより、上記 3 つのような数はすべて同じパターンとなり hash table によって重複要素は除去されます。
　n=5の場合、すべてのリストを ht に登録し重複除去された時点でパターン数は 4542個まで絞り込まれました。右辺 a の計算回数は、problem030a の 100分の1にまで減っています。その代わり前準備がかなり複雑になっています。

プロファイル

Function     Avg (sec)
===========  =========
problem030a  4.156500
problem030b  2.421500

　hash table 登録や リストのソートなど、オーバーヘッドになりそうな処理をしているにもかかわらず、後者の方式のほうが速いようです。ループ回数 100倍の差が効いているのでしょうか？