2020/06/14

え、エラトステネスの篩ですか？

素数とは

素数とは１より大きい自然数で、正の約数が１とそれ自身のみという数字のことで、私は学校での知識としての素数より、映画「博士の愛した数式」で寺尾聰さん演じる博士がもっとも愛した数字としての印象が強いです。この映画では素数の他にも「完全数」「友愛数」なども出てきます。この映画で、数学に興味を持った方は意外と多いのではないでしょうか。

今回はその博士に最も愛された素数についての問題をやってみます。これはこのところずっと勉強しているOCamlのL99の31問目が素数の問題だからですが、実は素数判定にはやり方が色々あるようで、OCamlのL99-31のドキュメントには「より賢い解決法として「エラトステネスの篩(ふるい)」というものがあるから調べてみて」と書かれています。

エラトステネスの篩をWikipediaで調べてみると、手順として

  1. 探索リストに2からxまでの整数を昇順で入れる
  2. 先頭の数を素数リストに移動し、その倍数を探索リストから消す
  3. 手順２を先頭値がxの平方根に達するまで続ける
  4. 探索リストに残った数と素数リストに移動する
という手順で行うそうです。平方根…何となく記憶にあるが…

まずはOCamlで作っていこうと思います。

手順１は以下のような関数を作成しました。


let rec range s e =
  if s > e then []
  else s :: range (s + 1) e
;;

sはスタート、eはエンドの略です。これでsを2とし、eを任意のxとすれば2からxまでの昇順のリストが作成できます。


# range 2 10;;
- : int list = [2; 3; 4; 5; 6; 7; 8; 9; 10]

次に手順２の先頭の数を素数リストに移動し、その倍数を消す、ですが、以下のような感じになりました。


let filter n l = List.filter (fun x -> x = n || x mod n <> 0) l

この関数を実行すると


# filter 2 (range 2 20);;
- : int list = [2; 3; 5; 7; 9; 11; 13; 15; 17; 19]

このようになります。結果をご覧のとおり2以外の2の倍数は除外されています。

これらを何とか組み合わせ、完成したコードが以下のようなものです。平方根云々はちょっと何言ってるか分からないので、何となくスルーしました。ここが高速化のツボのようなので、今後の課題として保留、要勉強。


exception Not_prime_number

let is_prime n =
  let rec range s e = if s > e then [] else s :: range (s + 1) e
  in
  let filter n l = List.filter (fun x -> x = n || x mod n <> 0) l
  in
  let rec make_primes acc = function
    | [] -> List.rev acc
    | h :: t -> make_primes (h :: acc) (filter h t)
  in
  let primes = make_primes [] (range 2 n)
  in
  if n < 2 then raise Not_prime_number
  else List.mem n primes
;;

結果は以下のように一応機能している感じです。


# is_prime 2;;
- : bool = true
# is_prime 3;;
- : bool = true
# is_prime 6;;
- : bool = false
# is_prime 97;;
- : bool = true

まとめ

これは果たしてエラトステネスの篩というやつになっているのだろうか？たぶんなっていないのだろう。流れとしては、まず与えられた数までの素数リストを作成し、そのリスト内に与えられた数があるかどうかを探す。うん、どうやらエラトステネスの篩ではないな。もうちょっと勉強しなければならないようです。平方根のあたりの処理がうまくできたらCommon Lispのほうでもやってみようと思います。ちなみにOCamlのL99-31の解答をCommon Lispでやってみると、以下のようなものになりました。


(defun primep (n)
  (let* ((nn (abs n)))
    (labels ((is-not-divisor (d)
               (or (> (* d d) nn)
                   (and (/= (mod nn d) 0)
                        (is-not-divisor (+ d 1))))))
      (and (/= nn 1)
           (is-not-divisor 2)))))

追記(2020/06/14)

私のプログラムを測定してみたところ、遅い！！

これはどう？合ってますか？というコードが以下


exception Not_prime_number

let rec range s e =
  if s > e then []
  else s :: range (s + 1) e

let filter n l = List.filter (fun x -> x mod n <> 0) l

let is_prime n =
  let list = range 2 n in
  let rec make_primes acc = function
    | [] -> acc
    | h :: t -> 
        if (sqrt (float_of_int h)) < (sqrt (float_of_int n))
        then make_primes (h :: acc) (filter h t)
        else h :: acc
  in
  let primes = make_primes [] list
  in
  if n < 2 then raise Not_prime_number
  else List.mem n primes
;;

平方根を添えて…

追記(2020/06/17)

間違いを修正しました。間違いの箇所は


    | h :: t -> 
        if (sqrt (float_of_int h)) < (sqrt (float_of_int n))
        then make_primes (h :: acc) (filter h t)
        else h :: acc

という部分です。

修正したものが以下


exception Not_prime_number

let rec range s e =
  if s > e then []
  else s :: range (s + 1) e

let filter n l = List.filter (fun x -> x mod n <> 0) l

let is_prime n =
  let list = range 2 n in
  let rec make_primes acc = function
    | [] -> acc
    | h :: t -> 
        if (float_of_int h) < (sqrt (float_of_int n))
        then make_primes (h :: acc) (filter h t)
        else List.append (h :: acc) t
  in
  let primes = make_primes [] list
  in
  if (n < 2) then raise Not_prime_number
  else List.mem n primes
;;

今回のコードでは# is_prime 111111という大きめの桁数でも遅くないので、うまくいっている感じがします。

これをCommon Lispで書いてみると、以下のようになりました。


(defun filter (n l)
  (loop for i in l when (/= (mod i n) 0)
        collect i))

(defun primep2 (n)
  (let ((l (loop for i from 2 upto n collect i)))
    (labels ((make-primes (acc ls)
               (cond
                 ((endp ls) acc)
                 (t
                  (if (< (car ls) (sqrt n))
                      (make-primes (cons (car ls) acc) (filter (car ls) (cdr ls)))
                      (append (cons (car ls) acc) (cdr ls)))))))
      (let ((primes (make-primes '() l)))
        (unless (< n 2)
            (if (eql (find n primes) nil)
                nil
                t))))))

こちらも* (primep2 111111)を大きな桁数でもControl stack ...... disabledとならないので、うまくいっているように思います。

追記(2020/07/02)

上記のprimep2はうまく機能していなかったようです。L99-39で使った時に発覚しました。とりあえずGitHubのファイルはコメントアウトしました。また改善していこう思います。

author: date:
Label: ,