素数のお勉強 4 - 約数の和

数学ガールを再読中デス。2章のミルカさんの宿題より。

正の整数nが与えられていたとき、nの約数の和を求めよ。

紙に書いて写経中。

・・・

アレ?なんで掛け算ナンダ。

(a + b)(c + d) = a(c + d) + b(c + d) = ac + ad + bc + bd

と。これが分配法則。中学くらいでやった。

4の約数[1,2,4]の和、3の約数[1,3]の和で掛け算すると、

(1 + 2 + 4)(1 + 3) = 1(1 + 2 + 4) + 3(1 + 2 + 4) = 1 + 2 + 4 + 3 + 6 + 12

スゲー。12の約数が出た!!不思議過ぎる。足したら28。

分配法則スゲー。(積の法則と言うらしい)

・・・

で、変形できるよね。

1 + 2 + 4 = 2^0 + 2^1 + 2^2

等比数列の和!!

素数pを0乗からn乗した和は。(pは1ではない。素数なので定義しなくても良さそうだけど)

p^0 + p^1 + p^2 + ... + p^n = (p^(n+1) - 1)/(p - 1)

っと。

・・・・

いや、数学ガールを読み込めてませんでした。反省。


1024と36の場合。

1024 = 2^10の約数の和なら等比数列の和から 2^11 - 1 = 2047
36 = 2^2 * 3^2なら、(2^3 - 1)/1 * (3^3 - 1)/2 = 96

素晴らしい。


素因数分解して、p^nの形にする。

(define (factor n)
  (letrec ((iter (lambda (i m)
                   (cond ((= m 1) '())
                         ((zero? (modulo m i)) (cons i (iter i (/ m i))))
                         (else
                            (iter (+ i
                                     (if (= i 2) 1 2))
                                  m))))))
          (iter 2 n)))

(factor 1024) ; (2 2 2 2 2 2 2 2 2 2)
(factor 36)   ; (2 2 3 3)

(define (count-list l)
  (if (null? l)
      '()
      (receive (a b)
               (partition (lambda (x) (equal? (car l) x)) l)
               (cons (cons (car l) (length a)) (count-list b)))))

(count-list (factor 1024)) ; ((2 . 10))
(count-list (factor 36))   ; ((2 . 2) (3 . 2))

一気にやっちゃってもいいんだけど2段がまえで。


素因数分解したリストに対して等比数列の和(p^n - 1)/(p - 1)の積で約数の和。

(define (divisor-sum n)
  (fold (lambda (l acc)
          (* acc
             (/ (- (expt (car l) (+ (cdr l) 1)) 1)
                (- (car l) 1))))
        1
        (count-list (factor n))))

(divisor-sum 1024) ; 2047
(divisor-sum 36) ; 91

以前書いた奴より断然はえぇ。(素因数分解した段階でp^nの形になってないと早くはないか)


等比数列の和から色々発展していくので、きっちり覚えておきたい。もう忘れないぞぉ〜。

おぉ、

著者からスター貰った。ちょっと嬉しい。

さて続き。

完全数を出してみよう」

factorをちょっと速度改善して、(> i (sqrt m))をぶち込んでおく。


完全数メルセンヌ素数から導き出される。

メルセンヌ数は、

メルセンヌ数 = 2^p - 1

pが素数の時メルセンヌ素数が出てくる(全てではない)。素数から素数が出てくるなんて素敵だ。


完全数メルセンヌ素数M(p)を元に、

完全数 = 2^(p-1) * メルセンヌ素数M(p) = p^(p-1) * (2^p - 1)

完全数の約数の和が元の数の2倍になったら完全数だ。


8番目のメルセンヌ素数M(31)で挑戦する。

M(31) = 2^31 - 1 = 2147483647

である。メルセンヌ素数M(31)から完全数を求め、素因数分解して約数の和を求め、完全数かどうか確かめるのだ!!

(define perfect-31 (* (expt 2 (- 31 1)) (- (expt 2 31) 1)))

perfect-31                       ; 2305843008139952128
(count-list (factor perfect-31)) ; ((2 . 30) (2147483647 . 1))
(/ (divisor-sum perfect-31) 2)   ; 2305843008139952128

うほぉ。一致。

8番目の完全数だと言うのに桁数はやたらでかい。完全数が 2のp-1乗 * メルセンヌ素数M(p) からしか生まれないと言うのは不思議な感じだ。

あ、そうそう。

これが言いたかった。

完全数、2305843008139952128」