素数列生成と素数判定のHaskellプログラム

こんにちは｜こんばんは。カエルのアイコンで活動しております @kyamaz です。

はじめに

次のエントリーにまとめてありますが、は少し前まで素数についてのエントリーを幾つか投稿してきました。

そこで素数生成と素数判定のHaskellのワンラインプログラムを使いました。それぞれ参考にしたサイトがありますので、はじめにプログラムと参照したサイトのリンクを挙げておきます。

素数生成のプログラム

素数生成のワンライナー

p=2:3:5#p where n#x@(m:q:y)=[n|gcd m n<2]++(n+2)#last(x:[m*q:y|q^2-3<n])

素数判定のプログラム

強引にワンラインにしていますので一行が長いです。

素数判定のワンライナー

is_prime' n | n<=1 = False | n==2 || elem n pp = True | otherwise = and $ map (\a -> t n a) pp where t m a | gcd m a < 2 = (1 == exmod a m) | otherwise = False; pp = [3,5,7,11,13,17,19]; exmod a m = itr a (m-1) m where itr a x b | x == 1 = a | even x = (itr a (x `div` 2) b) ^ 2 `mod` b | otherwise = a * (itr a ((x-1) `div` 2) b) ^ 2 `mod` b

参考サイトにもプログラムの処理は詳しく書かれていません。本稿では、このプログラムの処理について説明することと致します。

素数生成

「はじめに」に掲げた処理は、無理やりワンラインにしてありますので、これを改行つきに整形すると次のようになります。

素数生成プログラムを改行つきに整形

p=2:3:5#p             -- (A)
 where                -- (B)
   n#x@(m:q:y)        -- (C)
     = [n|gcd m n<2]  -- (D)
     ++(n+2)#last(x:[m*q:y|q^2-3<n])  -- (E)

各行に(A)〜(E)までのIDを振りました。それぞれの行について説明をしていきましょう。

(A)
素数列 p を定義します。
[2,3]++(5#p)と同値処理で、p が再帰的配列。二項演算子 (#) は (B) 以降で定義されます。

(B)
where 以降で、二項演算子 (#) を定義しています。

(C)
二項演算子 (#) を定義します。
この演算子の型は (#) :: Integral a => a -> [a] -> [a] となります。
(A) より各引数の初期値は次のとおり。
n <- 5 x <- p(素数列) m <- 2 q <- 3 y <- 5#p

(D)
[n|gcd m n<2]の処理は、対象の n が２つ前の素数 m と互いに素なら配列に含めます。
初期値の動作としては、5 は 2 と互いに素なので配列に含めます。

(E)
二項演算子 (#) を再帰的に使い、(n+2)#last(x:[m*q:y|q^2-3<n]) を配列に追加します。
初期値の動作としては、(#) 7 (last(p:[2*3:(5#p)|3^2-3<5]))となり、lastの中は[2*3:(5#p)|3^2-3<5]が条件部を満たさないため p です。
(#) 7 pは次の[n|gcd m n<2]の処理ができて 7 は 3 と互いに素なので 7 が配列に含まれます。
さらに再帰的処理を続けて見てみましょう。次の(n+2)#last(x:[m*q:y|q^2-3<n])が評価は、(#) 9 (last(p:[2*3:(5#p)|3^2-3<7]))となり、lastの中は [2*3,5,(7#p)] が決まり、(#) 9 [2*3,5,(7#p)]の[n|gcd m n<2]が処理できて 9 は 2*3 と互いに素ではないので、配列に含まれません。
n < 23 までは、2*3 と互いに素かで評価しますが、n <- 25のときに、2*3*5 と互いに素かで評価することになります。つまり、最後のlastのところではnより小さい素数の積を先頭の項にして、それ以降の項がnより大きい素数列を作っている式になります。

素数判定

次に、素数判定の処理について取り上げましょう。その前に、重要な定理であるフェルマーの小定理を紹介します。

フェルマーの小定理

すべての素数 $p$ において、$p$ とは互いに素である任意の $a \in \mathbb{Z}$ に対して、$a^{p-1} \equiv 1 \pmod{p} $ が成り立つ

この定理に対して、次のように対偶を考える。この対偶は正しい。

フェルマーの小定理の対偶

$n$ と互いに素な整数 $a$ が、$a^{n-1} \not\equiv 1 \pmod{n} $ の条件を満たせば $n$ は合成数である

この対偶は、『$n$ は合成数である』ことの十分条件であり、必要条件ではありません。つまり『$n$ と互いに素な整数 $a$ （ただし $2 \le a \lt n$）全てにおいて１つでも　$a^{n-1} \equiv 1 \pmod{n} $ を満たすとすると、その $n$ は合成数ではない』とは言えません。しかし、合成数ではない確率は高いということが知られており、確率的な素数判定の方法に用いられております。これを「フェルマーテスト」と呼びます。

このテストでは確率的に素数であろうということで『確率的素数』と呼ばれます。また、確率的素数と判定されるものを『擬素数』といい、擬素数は検証する整数 $a$ を変えると　$a^{n-1} \not\equiv 1 \pmod{n} $ が成り立つものですが、なかでも任意の $a$ に対して $a^{n-1} \equiv 1 \pmod{n} $ が成立するものもあり、『カーマイケル数』または『絶対擬素数』といいます。$561$ が最小のカーマイケル数であることは有名です。

フェルマーテストは確率的素数判定法ですが、決定的素数判定法もあります。詳しくはWikipediaを参考にしてください。（いずれプログラムを紹介するエントリも寄稿しよと思います。）

さて、それでは「はじめに」に掲げた処理を、次のように改行つきに整形して、各行を説明していきましょう。各行に(01)〜(15)までのIDを振りました。

素数判定プログラムを改行つきに整形

is_prime' n                                -- (01)
  | n<=1 = False                           -- (02)
  | n==2 || elem n pp = True               -- (03)
  | otherwise = and $ map (\a -> t n a) pp -- (04)
  where                                    -- (05)
    t m a                                  -- (06)
      | gcd m a < 2 = (1 == exmod a m)     -- (07)
      | otherwise = False;                 -- (08)
    pp = [3,5,7,11,13,17,19];              -- (09)
    exmod a m = itr a (m-1) m              -- (10)
      where                                -- (11)
        itr a x b                          -- (12)
          | x == 1 = a                     -- (13)
          | even x = (itr a (x `div` 2) b) ^ 2 `mod` b             -- (14)
          | otherwise = a * (itr a ((x-1) `div` 2) b) ^ 2 `mod` b  -- (15)

(01)
is_prime'関数の定義です。引数に検査する整数(Integer) n をとります。
確率的な判定関数ですので、関数名にプライム(')をつけております。

(02)
is_prime'関数のガード式です。
n<=1のときには False を返します。

(03)
is_prime'関数のガード式。
n==2 または、(09)で定義される素数列 pp に含まれていたら True を返します。

(04)
is_prime'関数のガード式。
and $ map (\a -> t n a) pp は、素数列 pp に対して テスト関数 t を行い、その結果が全て真（True）かを and をとって検証します。

(05)
where 以降で、関数 t、素数列 pp、関数 exmod を定義しています。

(06)
関数 t の定義です。引数に検査する整数(Integer) m、ある整数(Integer) a をとります。

(07)
関数 t のガード式で、条件である左辺 gcd m a < 2 は a と m が互いに素かを調べています。互いに素ならば、$a^{n-1} \equiv 1 \pmod{n} $ の真偽を返します。

(08)
関数 t のガード式で、a と m が互いに素でなければ False を返します。

(09)
素数列 pp = [3,5,7,11,13,17,19] は検証する整数 $a$ に相当します。全て既知の素数ですので、n とは互いに素になります。このリスト長を変えることにより、処理時間が左右しますし、擬素数を素数と判定してしまうかが左右されます。（後述）

(10)
関数 exmod の定義です。引数に検査する整数(Integer) m、ある整数(Integer) a をとります。右辺は、関数 itr です。

(11)〜(15)
where 以降で、関数 itr を定義しています。
この処理は、$a^{n-1} \pmod{n} $ を計算しています。計算量の工夫として、例えば、$a^{10} = (a^5)^2 = (a * a^4)^2 = (a * (a^2)^2)^2 = (a * (a*a)^2)^2$ と高々二乗を使うだけの計算にしています。また展開した途中で mod も計算されているので、計算量が軽減されています。

カーマイケル数

前出の参照サイトでも検証されていますが、カーマイケル数として以下が挙げられています。

[1729,2821,6601,8911,15841,29341,41041,46657,52633,63973,75361,101101]

参照サイトでは(03)のガード式がなくpp=[2,3,5]のプログラムになります。ここで紹介したプログラムの(09)を変更して擬素数が合成数と判定されるかを見てみましょう。

pp=[3,5]のとき

ghci> filter is_prime' [1729,2821,6601,8911,15841,29341,41041,46657,52633,63973,75361,101101]
[1729,2821,6601,8911,15841,29341,41041,46657,52633,63973,75361,101101]
-- 参照サイトと同様の結果になり、リストにある擬素数を素数と判定してしまう

pp=[3,5,7]のとき

ghci> filter is_prime' [1729,2821,6601,8911,15841,29341,41041,46657,52633,63973,75361,101101]
[29341,46657,75361]
-- 7を検査に加えると、擬素数を素数と判定する数は3つに減少

pp=[3,5,7,11]のとき

ghci> filter is_prime' [1729,2821,6601,8911,15841,29341,41041,46657,52633,63973,75361,101101]
[29341,46657]
-- 7,11を検査に加えると、擬素数は2つ

pp=[3,5,7,11,13]のとき

ghci> filter is_prime' [1729,2821,6601,8911,15841,29341,41041,46657,52633,63973,75361,101101]
[]
-- 7,11,13を検査に加えると、10 万以下の擬素数と判定される数は0個となりました

pp の素数列の長さを大きくすると擬素数と判定される数は減るようです。ここで紹介した例では、pp を20以下の既知の素数列としましたが、擬素数が混在することは変わりません。合成数 $n$ について、$n$ 未満の任意の素数 $p$ について、$a^{n-1} \equiv 1 \pmod{n} $ の真偽を調べるように変更するとよく、前出の素数列 p を使ってプログラムの(09)を pp = takeWhile (\q-> q*q<n) p に変更するとカーマイケル数の判定も上手くいきますが、ppの素数列が長くなるため計算時間とはトレードオフがあります。

ppを変更したワンライナー

 is_prime' n | n<=1 = False | n==2 || elem n pp = True | otherwise = and $ map (\a -> t n a) pp where t m a | gcd m a < 2 = (1 == exmod a m) | otherwise = False; pp = takeWhile (\q-> q*q<n) p; exmod a m = itr a (m-1) m where itr a x b | x == 1 = a | even x = (itr a (x `div` 2) b) ^ 2 `mod` b | otherwise = a * (itr a ((x-1) `div` 2) b) ^ 2 `mod` b

検証は以下。

ghci> filter is_prime' [561,1105,1729,2465,2821,6601,8911,10585,15841,29341,41041,46657,52633,62745,63973,75361,101101,115921,126217,162401,172081,188461,252601,278545,294409,314821,334153,340561,399001,410041,449065,488881,512461,530881,552721,656601,658801,670033,748657,825265,838201,852841,997633,1024651,1033669,1050985,1082809,1152271,1193221,1461241,1569457,1615681,1773289,1857241,1909001,2100901,2113921,2433601,2455921,2508013,2531845,2628073,2704801,3057601,3146221,3224065,3581761,3664585,3828001,4335241,4463641,4767841,4903921,4909177,5031181,5049001,5148001,5310721,5444489,5481451,5632705,5968873,6049681,6054985,6189121,6313681,6733693,6840001,6868261,7207201,7519441,7995169,8134561,8341201,8355841,8719309,8719921,8830801,8927101,9439201,9494101,9582145,9585541,9613297,9890881]
[]

おわりに

いかがでしたでしょうか。ご一読いただきまして有り難うございます。

最後に、本稿を記載するために検証したHaskell環境を記しておきます。お手元の環境で検証する際に、動作が異なるときには参考になるかもしれません。

本稿の環境

本稿のために使用した環境は以下となります。

• macOS: Sonoma 14.6.1 (chip: Apple M1)
• GHCup: 0.1.30.0
• GHC: 9.6.6

（●）（●） Happy Hacking!
/""　＿＿""＼