【Matlab/Octave】素数を列挙するアルゴリズム

令和6年度秋期応用情報技術者試験午後問題問3「素数を列挙するアルゴリズム」を Matlab/Octave で記載してみた。

関数 prime1

2 以上の自然数 N に対して，N 以下の素数を列挙する関数 prime1 のプログラムを示す。

function prime1(N)
  primes = [];
  d = 2;
  calc_num = 0;
  while (d <= N)
    isPrime = true;
    t = 2;
    while (t < d)
      if (mod(d, t) == 0)
        isPrime = false;
      endif
      t = t + 1;
      calc_num += 1; 
    endwhile
    if (isPrime == true)
      primes = [primes d];
    endif
    d = d + 1;
  endwhile
  primes
  calc_num
endfunction

N = 100;
prime1(N)

関数 prime1 のコードを実行した結果を示す。

primes =

 Columns 1 through 23:

    2    3    5    7   11   13   17   19   23   29   31   37   41   43   47   53   59   61   67   71   73   79   83

 Columns 24 and 25:

   89   97

calc_num = 4851

N = 100 の場合，実行回数は 4,851 回である。

関数 prime2

素数の定義によって，2 以上の自然数 s について，s 自身を除く s の正の倍数 u は，1 と u 以外にも s も約数に含むので素数ではない。この性質を利用して関数 prime1 を改良し，次の手順で素数を列挙する関数 prime2 を考える。

(1) 2 以上 N 以下の自然数について，全て “素数” であるとマークする。

(2) 2 以上 N 以下の自然数 d について，次の a.，b. を行う。

d が “素数でない” とマークされている場合，何もしない。
d が “素数である” とマークされている場合，次の処理を行う。
1. d が素数であることを確定させる。
2. d 以上の自然数 x について，d を x 倍した数を “素数ではない” とマークする。

関数 prime2 のプログラムを示す。

function prime2(N)
  primes = [];
  isPrime = [false];
  c = 2;
  d = 2;
  calc_num = 0;
  while (c <= N)
    isPrime = [isPrime true];
    c = c + 1;
  endwhile
  while (d <= N)
    if (isPrime(1,d) == true)
      primes = [primes d];
      t = d * d;
      while (t <= N)
        isPrime(1,t) = false;
        t = t + d;
        calc_num += 1;
      endwhile
    endif
    d = d + 1;
  endwhile
  primes
  calc_num
endfunction

N = 100;
prime2(N)

関数 prime2 のコードを実行した結果を示す。

primes =

 Columns 1 through 23:

    2    3    5    7   11   13   17   19   23   29   31   37   41   43   47   53   59   61   67   71   73   79   83

 Columns 24 and 25:

   89   97

calc_num = 104

N = 100 の場合，実行回数は 104 回である。

関数 prime3

4 以上の偶数は全て 2 の倍数であるので素数ではない。したがって，2 以外の素数を列挙するためには奇数だけを考慮すればよい。この性質を利用して，関数 prime2 に次の変更を加えた関数 prime3 を考える。

関数の戻り値として素数の一覧が格納される primes にあらかじめ 2 を格納しておく。
いずれのループも奇数についてだけ実行されるようにする。
3 以上の自然数 2k+1 が素数か否かを isPrime[k] で表すようにする。

関数 prime3 のプログラムを示す。

function prime3(N)
  primes = [2];
  isPrime = [];
  c = 3;
  d = 3;
  calc_num = 0;
  while (c <= N)
    isPrime = [isPrime true];
    c = c + 2;
  endwhile
  while (d <= N)
    if (isPrime(1,int8((d - 2) / 2)) == true)
      primes = [primes d];
      t = d * d;
      while (t < N)
        isPrime(1,int8((t - 2) / 2)) = false;
        t = t + 2 * d;
        calc_num += 1;
      endwhile
    endif
    d = d + 2;
  endwhile
  primes
  calc_num
endfunction

N = 100;
prime3(N)

関数 prime3 のコードを実行した結果を示す。

primes =

 Columns 1 through 23:

    2    3    5    7   11   13   17   19   23   29   31   37   41   43   47   53   59   61   67   71   73   79   83

 Columns 24 and 25:

   89   97

calc_num = 28

N = 100 の場合，実行回数は 28 回である。