こんにちは｜こんばんは。カエルのアイコンで活動しております @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!
/"" ＿＿""＼