😸

Gauche で素数で遊ぶ

3 min read

Gauche のmath.prime モジュールで遊びます。

はじめに

以下のモジュールを use します。

(use math.prime)
(use scheme.list) ; take-while

math.prime モジュールのドキュメントには、各関数の入力範囲について言及がありますが、本記事では小さな正整数のみ扱うこととし、大きな数については考慮しませんのでご注意ください。

N 以下の素数を求める

gosh> (take-while (^p (<= p 100)) (primes))
(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97)

小さな値から順に N 弧の素数を求める

gosh> (take (primes) 25)
(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97)

N 番目の素数を求める

gosh> (list-ref (primes) 0)
2
gosh> (list-ref (primes) 1)
3
gosh> (list-ref (primes) 3)
7

N を素因数分解する

gosh> (mc-factorize 10)
(2 5)
gosh> (mc-factorize 2600)
(2 2 2 5 5 13)

N が素数かどうか調べる

gosh> (small-prime? 13)
#t
gosh> (small-prime? 1024)
#f

N 以上の最小の素数を求める

gosh> (find (^p (> p 30)) (primes))
31

N 以下の最大の素数を求める

gosh> (last (take-while (^p (< p 30)) (primes)))
29

N 以上の素数を、小さな順から M 弧求める

gosh> (take (drop-while (^p (< p 30)) (primes)) 3)
(31 37 41)

約数の個数を求める

n を素因数分解して次のように表せるとき、

n = p_1^{e_1} \cdot p_2^{e_2} \cdots p_m^{e_m}

n の約数の個数は、次式で求めることが出来ます。

(e_1 + 1) \cdot (e_2 + 1) \cdots (e_m + 1)
(define (number-divisors n)
  ;; (素数 . 個数) のペアを返す
  (define (fac n)
    (let ((ht (make-hash-table)))
      (dolist (prime (mc-factorize n))
        (hash-table-put! ht prime (+ 1 (hash-table-get ht prime 0))))
      (hash-table->alist ht)))

  (fold * 1 (map (^a (+ (cdr a) 1)) (fac n))))

実行結果です。

gosh> (number-divisors 100)
9
gosh> (number-divisors 1024)
11

約数の総和を求める

n を素因数分解して次のように表せるとき、

n = p_1^{e_1} \cdot p_2^{e_2} \cdots p_m^{e_m}

n の約数の総和は次式で表せます。

(1 + p_1 + p_1^2 + \cdots + p_1^{e_1}) \cdot (1 + p_2 + p_2^2 + \cdots + p_2^{e_2}) \cdot \cdots \cdot (1 + p_m + p_m^2 + \cdots + p_m^{e_m})

等比数列の和の公式より、上の式は

\frac{p_1^{e_1 + 1} - 1}{p_1 - 1} \cdot \frac{p_2^{e_2 + 1} - 1}{p_2 - 1} \cdot \cdots \cdot \frac{p_m^{e_m + 1} - 1}{p_m - 1}

となります。

(define (sum-divisors n)
  ;; (素数 . 個数) のペアを返す
  (define (fac n)
    (let ((ht (make-hash-table)))
      (dolist (prime (mc-factorize n))
        (hash-table-put! ht prime (+ 1 (hash-table-get ht prime 0))))
      (hash-table->alist ht)))

  (fold * 1 (map (^a (let ((p (car a))
                           (e (cdr a)))
                       (/ (- (expt p (+ e 1)) 1)
                          (- p 1))))
                 (fac n))))

実行結果です。

gosh> (sum-divisors 12)
28
gosh> (sum-divisors 150)
372

約数を列挙する

;; n の約数を求める
(define (divisors n)
  (let ((ht (make-hash-table)))
    (do ((i 1 (+ i 1)))
        ((> (* i i) n))
      (when (zero? (modulo n i))
        (hash-table-put! ht i #t)
        (hash-table-put! ht (/ n i) #t)))
    (sort (hash-table-keys ht))))

実行結果です。

gosh> (divisors 12)
(1 2 3 4 6 12)
gosh> (divisors 13)
(1 13)
gosh> (divisors 1024)
(1 2 4 8 16 32 64 128 256 512 1024)

Discussion

ログインするとコメントできます