読者です 読者をやめる 読者になる 読者になる

Problem 78

少し前から一部で Project Euler というのが流行しているようだ。
http://projecteuler.net/index.php?section=problems
日本語訳を試みているサイトもある。
http://odz.sakura.ne.jp/projecteuler/index.php?Project%20Euler
言うなればプログラミングの問題集と言ってよいだろう。 しかし、ただのプログラミング演習ではなく、 Euler の名を冠するだけあって数学的な性質を利用しないと解けないものもある。 Problem 78 などは特に難しいという記事 (id:umeaji:20080919:p1) を見掛けて私も Gauche でやってみることにした。
分配関数 p について p(n) が百万で割り切れる n を求めよという問題であるが、とりあえず百万で割るのは置いて、 p を私なりに愚直に実装するとこうなった。

(define (p n)
  (rlet1 t 0
    (let loop ((m n)(r n))
      (if (zero? r)
          (inc! t)
          (do ((c (min r m)(dec! c)))
              ((zero? c))
            (loop c (- r c)))))))

激しく遅い。 p(80) を算出するだけで1分を越えてしまう。
全く話にならないので、途中経過をハッシュテーブルに格納するようにして高速化してみたのが以下だ。

(define p
  (let1 ht (make-hash-table 'equal?)
    (lambda(n)
      (let loop ((m n)(r n))
        (if (zero? r)
            1
            (begin
              (set! m (min r m))
              (if-let1 e (hash-table-get ht (cons m r) #f)
                e
                (rlet1 t 0
                  (do ((c (min r m) (- c 1)))
                      ((zero? c))
                    (set! t (+ t (loop c (- r c)))))
                  (hash-table-put! ht (cons m r) t)
                  ))))))))

これなら p(80) 程度は一瞬で終わる。 かなりの高速化だ。 p(1) から p(1000) までを算出しても私の環境では18分ほどで完了した。
とは言っても、これでもまだまだ話にならない遅さである。 計算時間の増加傾向を見ると、 p(2000) とかにしたらやはりとんでもない時間がかかりそうだ。 やはり小手先の工夫ではどうにもならない。
と、いうわけで他にどんな方法で解いている人がいるのかとウェブ検索してみたところ、 Python での例を発見した。
http://hi.baidu.com/patrolsun/blog/item/3592a51c6dfefe8d86d6b629.html
Python を使ったことはないが、それほど難しくはないので Gauche にそのままベタ移植してみたのが以下だ。

(define p
  (let1 ht (make-hash-table 'eqv?)
    (lambda(n)
      (cond ((zero? n) 1)
            ((negative? n) 0)
            ((hash-table-get ht n #f) => values)
            (else
             (let loop ((k 1) (count 0))
               (let ((n1 (- n (/ (- (* 3 k k) k) 2)))
                     (n2 (- n (/ (+ (* 3 k k) k) 2))))
                 (let ((count
                        (+ count
                           (* (if (even? k) -1 1)
                              (+ (p n1) (p n2))))))
                   (if (<= n2 0)
                       (begin
                         (hash-table-put! ht n count)
                         count)
                       (loop (+ k 1) count))))))))))

なんという爆速。
読んでもどうしてこれで問題が解けるのかわからないのだが、数学的な性質を用いることでここまで高速な計算が可能というのは驚きである。
これで安心して p(n) が百万で割切れるものを探すことが出来る。

(do ((x 1 (+ 1 x)))
    ((zero? (modulo (p x) 1000000)) x))

1分ちょっとで終了した。 答えは n=55374 であった。
最初の愚直な実装ではとてもじゃないが到達できない値である。
数学重要。
Document ID: 1eb01dff7df57563b9a1d0143aa32a84