Project Euler: Problem 12 #1

The sequence of triangle numbers is generated by adding the natural numbers. So the 7^(th) triangle number would be 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28. The first ten terms would be:

1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

Let us list the factors of the first seven triangle numbers:

1: 1
3: 1,3
6: 1,2,3,6
10: 1,2,5,10
15: 1,3,5,15
21: 1,3,7,21
28: 1,2,4,7,14,28

We can see that 28 is the first triangle number to have over five divisors.

hat is the value of the first triangle number to have over five hundred divisors?

http://projecteuler.net/index.php?section=problem&id=12

自然数 n に対して、1 から n までの合計を求めたものを三角数と呼ぶ。三角数の数列の中から 500 個の因数に分解できる最初の三角数を求めよとのこと。

とりあえず、Gauche で書いてみたけど実行が終わらない。現状は三角数を順番に求めて因数に分解しているけど、そんな方法で答えを求めようにも時間ばっかりかかってどうしようもないのだろう。ちょっくらアルゴリズムを練り直す必要があるが…どうすりゃいいかなぁ。

#!/usr/bin/env gosh

(use srfi-1)
(use util.combinations)

;; n 以下の素数のリストを返す
(define (primes n) ; sieve of eratosthenes (list version)
  (case n
    [(1) '()]
    [else
      (let loop [(i 0)
                 (m (floor->exact (sqrt n)))
                 (p (iota (- n 1) 2))]
        (if (< i m)
          (let* [(e (list-ref p i))
                 (p (filter (lambda (x)
                              (or (= x e)
                                  (not (= (remainder x e) 0))))
                            p))]
            (loop (+ i 1) m p))
          p))]))

;; n 番目の三角数を求める
(define (triangle-number n)
  (let loop [(n n)
             (r 0)]
    (case n
      [(1) (+ r 1)]
      [else
        (loop (- n 1) (+ r n))])))

;; n を素因数分解した結果のリストを返す
(define (prime-factors n)
  (case n
    [(1) '()]
    [(2) '(2)]
    [else
      (let loop [(num n)
                 (plis (primes n))
                 (flis '())] ; factors
        (if (< 1 num)
          (if (pair? plis)
            (if (zero? (remainder num (car plis)))
              (loop (quotient num (car plis))
                    plis
                    (cons (car plis) flis))
              (loop num
                    (cdr plis)
                    flis))
            (error "no prime left for dividing"))
          flis))])) 

;; n の全ての因数のリストを返す
(define (factors n)
  (case n
    [(1) '(1)]
    [else
      (let* [(f (cons '(1)
                      (cons (list n)
                            (power-set (prime-factors n)))))
             (g (map (cut fold * 1 <>)
                     f))
             (h (delete-duplicates (sort g)))]
        h)]))

(define (main args)
  (let loop [(n 1)]
    (let* [(tr (triangle-number n))
           (rf (factors tr))]
      (if (< 500 (length rf))
        [begin
          (newline)
          (newline)
          (print #`"the answer is ,tr ,rf.")]
        [begin
          (format #t "~16@a: ~a\n" tr (string-join (map number->string rf) ", "))
          (loop (+ n 1))]))))

つづく。