完全数を求める？（Haskell実装）

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

はじめに

本稿では、完全数について取り上げます。タイトルに「？」をつけた理由は、2024年5月現在、知られている完全数の個数は51個¹ですが、探索的に愚直な計算では、現実的な時間に計算できたのが小さい完全数から4個しか求められなかったからです。本稿の後半では、メルセンヌ素数との関係から26個程度までは求められる方法を紹介します。

完全数

完全数の定義

その数を除く正の約数の和がその数に等しくなる自然数

言い換えると「その数の正の約数の総和がその数の２倍になる自然数」とも言えます。約数関数$\sigma _x (n)$を使うと「$\sigma _{1} (a)=2a$を満たす自然数$a$」が完全数です。
Wikipedia『完全数』にも詳しく解説があります。そのなかにも記載がありますが 「偶数の完全数は無数に存在するか？」「奇数の完全数は存在するか？」という問題は未解決問題 です。

それでは、愚直に定義から完全数を求める計算を考えましょう。
与えられた数nが完全数かどうかを判定するisPerfect関数を定義します。nの約数のリストを求める必要がありますが、少し効率的にリスト化するために、$2$〜$\lfloor \sqrt{n} \rfloor$までの約数のリストをlows = filter ((0 ==) . mod n) [2 .. floor (sqrt (fromIntegral n))]でリストアップし、そのリストの数でnを割ったリストが、$1$とその数nを除く約数のリストになります。そのリストにある数の和に$1$を足した値がその数nであれば完全数です。

ghci> isPerfect n = 1 < n && n == (sum (lows ++ [y | x <- lows, let y = n `div` x])) + 1 where lows = filter ((0 ==) . mod n) [2 .. floor (sqrt (fromIntegral n))]
ghci> perfects ns = filter isPerfect ns

調べたい整数リストnsを引数にとってfilter関数で完全数を取り出します。

ghci> perfects [1..] 
[6,28,496,8128    -- この先はかなりの時間が経っても表示されません

8128まではすぐに出力されます。1772年にオイラー（Leonhard Euler）が8番目のメルセンヌ素数を発見しており、8番目の完全数を見つけていることになります。オイラーさんの計算力は尋常でないですね。
さて、我々は既に51番目まで知っているので、試しに既知¹の９番目までの完全数のリストperfsを作り、それが完全数かのチェックをしてみましょう。

ghci> perfs = [6, 28, 496, 8128, 33550336, 8589869056, 137438691328, 2305843008139952128, 2658455991569831744654692615953842176]
ghci> map isPerfect perfs
[True,True,True,True,True,True,True,True,
-- 8個目は待つものの現実的な時間で結果が表示されますが
-- 最後（9個目）のチェックはかなりの時間が経っても表示されません

残念ながら、オイラーさんが発見した8番目のところまでしか確認できませんでした。このまま4番目だけの計算だけでは不甲斐なさすぎます。そこで効率の良い手法を探していたところ、C言語でのプログラムを見つけました。本稿ではHaskell実装を目指していましたが、実装するプログラミング言語を変えて試してみましょう。

% time ./perfects 1 100000000
6
28
496
8128
33550336
./perfects 1 100000000  116.64s user 4.21s system 98% cpu 2:02.37 total

数分間の計算で、5番目の完全数までは計算できました。

奇数の完全数

奇数の完全数が存在するかは未解決問題で、もし奇数の完全数$N$があった場合は $N \gt 10^{1500}$ であることが Ochem and Rao, 2012 により示されています。これまで研究されている内容は、下記サイトに詳しく解説があります。

偶数の完全数

定理１

$2^{N}-1$が素数であるような正の整数$N$に対して、$n=2^{N-1}(2^{N}-1)$は完全数である。

定理２

偶数の完全数は、$2^{N}-1$が素数であるような正の整数$N$を用いて、$n=2^{N-1}(2^{N}-1)$という形で表される。

$2^{N}-1$の形の素数をメルセンヌ素数と言います。つまり、この定理から「メルセンヌ素数」と「偶数の完全数」が一対一に対応しているということになります。この定理の証明は下記のサイトにも紹介があり、高校数学の範囲でもできます。

「メルセンヌ素数は無数に存在するか？」という問題は未解決問題です。偶数の完全数にはそれと一対一対応があるため「偶数の完全数が無数に存在するのか？」という問題は未解決問題です。

メルセンヌ素数

ここからは、メルセンヌ素数を求めることに注力しましょう。まず、次のサイトを参考にフェルマーの小定理の対偶を使います。

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

正の整数$n$と互いに素な整数$a$に対して、$$a^{n-1} \not\equiv 1 \pmod{n}$$を満たすならば、$n$は素数ではない（合成数である）

これを用いた$n$が素数かどうかの（確率的な）判定法を フェルマーテスト と呼びます。²³この判定法を使ってメルセンヌ素数を作っていきましょう。

まず、いつものように素数リストを定義します。

ghci> p=2:3:5:5#p;m#(a:b:x)=[n|n<-[a^2..b^2-2],(mod (n+1) 6)*(mod (n-1) 6)==0,gcd n m<2]++(m*b)#(b:x)

次にフェルマーの小定理の対偶で素数判定を実装してみましょう。

ghci> exp_mod a n = iter a (n-1) n where iter a n base | n == 1 = a | even n = (iter a (n `div` 2) base) ^ 2 `mod` base | otherwise = a * (iter a ((n-1) `div` 2) base) ^ 2 `mod` base
ghci> is_prime n = and $ map (\a -> test n a) [2..(n-1)] where test n a | gcd n a == 1 = (1 == exp_mod a n) | otherwise = False

この関数is_primeでは、$a$が$2$から$n-1$まで対偶を確認しているため素数を判定できていますが、確認する$a$を少なく絞って確率的に素数判定する関数probably_primeにして試してみます。[2,3,5,7,11,13,17,19]の8個までの素数で確認するように変更してみます。

gchi> probably_prime n = and $ map (\a -> test n a) $ takeWhile (<n) [2,3,5,7,11,13,17,19] where test n a | gcd n a == 1 = (1 == exp_mod a n) | otherwise = False
ghci> filter (\p -> probably_prime $ 2^p -1) p
[2,3,5,7,13,17,19,31,61,89,107,127,521,607,1279,2203,2281,3217,4253,4423,9689,9941,11213,19937,21701,23209
-- 26番目までのメルセンヌ素数の生成のための素数が得られました

しばらく待つと、26番目までのメルセンヌ素数の生成のための素数が得られました。ここから先は時間がかかりましたので適当なタイミングでCtrl-Cして処理を止めましょう。
さらに、フェルマーテストよりも効率的とされている「リュカ・レーマーテスト」を使う方法⁴も紹介します。詳しくは、Wikipedia: Lucas–Lehmer primality testやProofWiki: Lucas-Lehmer Testや脚注のサイト⁵を参考にしてください。

ghci> lucasLehmer p | p==2 = True | otherwise = s (2^p-1) (p-1) == 0 where s mp n | n==1 = 4 `mod` mp | otherwise = ((s mp $ n-1)^2-2) `mod` mp
ghci> filter lucasLehmer p
[2,3,5,7,13,17,19,31,61,89,107,127,521,607,1279,2203,2281,3217,4253,4423,9689,9941,11213,19937,21701,23209

こちらのリュカ・レーマーテストでも、26番目までの計算にかなりの時間を要しますので、ご注意ください。

おわりに

上記の方法でメルセンヌ素数を生成する素数を求めていけばよいですが、51個は既知ですからそのリストを使って完全数のリストを定義して終わりにしたいと思います。

ghci> gen_mersenne=[2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423, 9689, 9941, 11213, 19937, 21701, 23209, 44497, 86243, 110503, 132049, 216091, 756839, 859433, 1257787, 1398269, 2976221, 3021377, 6972593, 13466917, 20996011, 24036583, 25964951, 30402457, 32582657, 37156667, 42643801, 43112609, 57885161, 74207281, 77232917, 82589933 ]
ghci> perfects = [2^(n-1) * (2^n - 1) | n <- gen_mersenne ]

26番目程度であれば即座に出力されます。

ghci> take 26 perfects
[6,28,496,8128,33550336,8589869056,137438691328,2305843008139952128,2658455991569831744654692615953842176,191561942608236107294793378084303638130997321548169216,13164036458569648337239753460458722910223472318386943117783728128,14474011154664524427946373126085988481573677491474835889066354349131199152128,...

また、51個目の完全数も次のように出力することができます。画面出力は比較的すぐに始まりますが、数字を出力するだけでも膨大な時間が必要です。

ghci>  perfects !! 50
110847779864148441098176280120360785502254663219598207290510 ...

は、画面出力を諦めて途中で処理を停止しました。コンパイルしてファイル出力するなどして、試すとよいかもしれませんが、くれぐれもディスク容量等のリソースには注意して自己責任でお試しください。

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

本稿の環境

本稿のために使用した環境は以下となります。 macOS: Sonoma 14.4 (chip: Apple M1) GHCup: 0.1.22.0 GHC: 9.6.4

ご一読いただきまして誠に有り難うございます。
（●）（●） Happy Hacking!
/""　＿＿""＼

1. Wikipedia: List_of_Mersenne_primes_and_perfect_numbers、2024年1月現在48番目までは確認されており、49番目以降の順番は検証で確定されておりません。 ↩ ↩²

2. https://math-fun.net/20210225/11311/ ↩

3. https://derbuihan.github.io/posts/haskell_primes/ ↩

4. https://itchyny.hatenablog.com/entry/20120616/1339835082 ↩

5. リュカ–レーマー・テストの証明
   リュカ–レーマー–リーゼル・テスト ↩