Problem 12

素数の時間です（謎

Problem 12 - PukiWiki

三角数の数列は自然数の和で表わされ、7番目の三角数は 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28 である。三角数の最初の10項は
1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

となる。

最初の7項について、その約数を列挙すると、以下のとおり。

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

これから、7番目の三角数である28は5つ以上の約数をもつことが分る。

では、501 個以上の約数をもつ最初の三角数はいくらか。

501個以上の約数って！！死ねる。

とりあえずナイーブな奴でお試し。

(define (problem12 n)
  (define (triangle-sum n)
    (apply + (iota n 1)))
  (define (divisor n)
    (filter (lambda (m)
              (= (modulo n m) 0))
            (iota n 1)))
  (define (iter x)
    (if (>= (length (divisor (triangle-sum x))) n)
        x
        (iter (+ x 1))))
  (iter 2))

(problem12 5) ; 7
(problem12 100) ; 384
(problem12 501) ; 終わらず

(実際に求めるのは三角数)

そうだよね。終わるはずないよね〜。中に素数(またはそれに近いもの)が含まれてると終わらないよね〜orz

ということで。ずっと前に作った素数キャッシュ付き素因数分解 + 約数を利用することに。

(define *primers-cache* '(2))

(define (factor n)
  (define (primers-filter n p)
    (remove
      (lambda (x)
        (zero? (modulo x n)))
      p))
  (define (make-primers p)
    (fold-right
      (lambda (a b)
        (primers-filter a b))
      (iota (car p) (+ (car p) 1))
      p))
  (define (factor-iter m p1 p2)
    (define (next p)
      (factor-iter m (cons (car p) p1) (cdr p)))
    (cond ((= m 1) '())
          ((zero? (modulo m (car p1)))
           (cons (car p1)
                 (factor-iter (/ m (car p1)) p1 p2)))
          ((null? p2) 
           (next (let ((c (make-primers p1)))
                       (set! *primers-cache* (append (reverse p1) c))
                       c)))
          (else
            (next p2))))
  (factor-iter n '(2) (cdr *primers-cache*)))

(require (lib "list.ss"))

(define (divisor n)
  (define (unique list)
    (fold-right
      (lambda (a b)
        (if (member a b)
            b
            (cons a b)))
      '() list))
  (define (subsets s)
    (if (null? s)
        (list '())
        (let ((rest (subsets (cdr s))))
             (append rest (map
                            (lambda (x)
                              (cons (car s) x))
                            rest)))))
  (sort
    (unique (map (lambda (m)
                   (fold * 1 m))
                 (subsets (factor n))))
    <))

(define (problem12 n)
  (define (triangle-sum n)
    (apply + (iota n 1)))
  (define (iter x)
    (if (>= (length (divisor (triangle-sum x))) n)
        (triangle-sum x)
        (iter (+ x 1))))
  (iter 2))

(problem12 6) ; 28
(problem12 501) ; 76576500

素数生成を頑張りすぎて、肝心の約数が微妙だったり（汗

1分位かかるので、少ない数で練習してキャッシュを貯めておくとイライラせずに済む。

factorしてみると。

% factor 76576500
76576500: 2 2 3 3 5 5 5 7 11 13 17

すくな！！

11個に素因数分解されて、約数は500個を越える。

つまり、素因数分解であぁ、コラ無理かも・・・と思われたら、大部分は無視して良さげな感じもする。この問題用に最適化すれば、かなり早くなると思う。

いろいろな回答を見て回ると

すげー落ち込む。

頑張らねば。

追記

バグがあったので修正。

高速化のためのメモ

約数の個数について。

• 約数 - Wikipedia