ウォルステンホルム素数（Haskellで出力する）

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

はじめに

本稿は「ウォルステンホルム素数（Wolstenholme prime）」を出力する処理を、Haskellで実装するShortエントリーです。すでにヴィーフェリッヒ素数を出力するエントリをご紹介しましたが、本稿でも似たような処理をご紹介します。

紹介する処理は、非常に時間を要する処理です。お試しになる際には、お時間があるときをお勧めします。

ウォルステンホルム素数

ウォルステンホルム素数は、@tsujimotterさんのYouTubeでも紹介されています。

2007年時点では$10^9$までの範囲で調べられており、既知のウォルステンホルム素数は、$16843$と$2124679$のみです。

いつくかの同値の定義がありますが、二項係数による定義を次に紹介します。

ウォルステンホルム素数の定義

${}_{2p-1} C_{p-1} \equiv 1 \pmod{p^{4}}$ を満たす素数$p$

ウォルステンホルム素数は無限個存在すると予想されており、また$x$以下のウォルステンホルム素数の個数は約$\ln{\ln{x}}$個と予想されておりますが、いずれも未解決問題です。

それでは、ウォルステンホルム素数をHaskellで出力する処理を作っていきましょう。
まずは、定義にコンビネーション（組合せ）が出てきますので、${}_n C_m$の計算を記述しましょう。${}_n C_m = \frac{{}_n P_m}{{}_m P_m}$であることを踏まえて実装します。

-- 順列（permutation）の計算
ghci> perm n m | n<m = 1 | n<=0 = 1 | m<=0 = 1 | otherwise = foldr (*) 1 [(n-m+1)..n]
-- 組合せ（combination）の計算
ghci> comb n m | n<m = 1 | n<=0 = 1 | m<=0 = 1 | otherwise = perm n m `div` (perm m m)

次に、エラトステネスの篩を踏まえて素数列を実装します。

-- 素数列 p
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)

ウォルステンホルム素数の定義を素直に実装すると、次のようにウォルステンホルム素数列 w ができます。

-- ウォルステンホルム素数
ghci> w=[n|n<-p, mod (comb (2*n-1) (n-1)) (n^4) == 1]

定義から実装を進めているだけですので、ここまではすぐに処理が終わります。さて、実際の計算を行ってみましょう。ウォルステンホルム素数は現在２個まで知られていますので、無邪気にウォルステンホルム素数列から最初の２つを取り出す処理として実行してみましょう。

ghci> take 2 w
[16843

だいぶ時間はかかりましたが、最初の$16843$は出力されました。その次の値を待ちます。しかし、かなりの時間（数時間）放置しましたが、残念ながら$2124679$まで到達しませんでしたので、Ctrl-Cで処理を中断しました。Wikipediaによると、この２つめの素数は、ほぼ20年間発見されなかったそうです。やはり手元のPCではなかなか酷な計算のようです。

おわりに

ご一読いただきまして有り難うございます。
試しに以下も試しましたが、一晩おいても結果が出力されませんでした。チャレンジして出力されるようでしたら、教えて頂きたいです。

ghci> mod (comb (2*2124679-1) (2124679-1)) (2124679^4)

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

本稿の環境

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

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