こんにちは｜こんばんは。カエルのアイコンで活動しております @kyamaz 
です。
はじめに
本稿は、「ヴィーフェリッヒ素数1（Wieferich prime）」を出力する処理を、Haskellを使ってワンライナーで実装するネタ的2なShortエントリーです。
ヴィーフェリッヒ素数とは
@tsujimotterさんのはてなブログエントリ『数学者の名前がついた素数』にヴィーフェリッヒ素数の解説があります。また、次のYouTube動画でも紹介されています。
これらの情報を参考にヴィーフェリッヒ素数の定義を以下に記します。
- ヴィーフェリッヒ素数の定義
- $2^{p-1} \equiv 1 \pmod{p^{2}}$ を満たす素数$p$
ここで、関係する次の『フェルマーの小定理』が一般的に知られています。
- フェルマーの小定理
- すべての素数$p$において、任意の$a \in \mathbb{Z}$（$a$と$p$は互いに素）とすると、
$a^{p-1} \equiv 1 \pmod{p} $ が成り立つ
このフェルマーの小定理の$a=2$において、$mod$を$p^2$に置き換えた条件が成立する素数をヴィーフェリッヒ素数と呼びます。実は、ヴィーフェリッヒ素数は1093と3511の２つしか見つかっていません。また、これ以上あるかも分かっておらず、いわゆる未解決問題です。
素数列を生成する
テーマとしては素数を扱いますので、まずはHaskellで素数列の生成を実装しましょう。素数を生成するプログラムは色々なサイトで紹介されており、ChatGPTでも教えてもらえますが、このはてなブログで紹介されている手法を参考にしてワンライナーで生成します。
お試しするには、GHCでプログラムファイルをコンパイルするのではなく、GHCiでプロンプトから対話で進めます。
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)
この処理は、単にエラトステネスの篩を使って処理しているだけですが、この記述は少し分かりづらいかもしれません。生成された素数が適切かをいくつかピックアップして確認しておきましょう。
生成した素数列を確認する
確認のため10,000番目、100,000番目、1,000,000番目の素数を表示してみます。検索できる素数リストやChatGPTで調べられる数と同じであることを確認しております。
ghci> p!!9999
104729
ghci> p!!99999
1299709
ghci> p!!999999
15485863
ヴィーフェリッヒ素数列を生成する
pが素数列ですのでpを使って、その中からヴィーフェリッヒ素数の条件をもつものを次のように内包表記でリストwを記述します。
ghci> w=[n|n<-p, mod (2^(n-1)) (n^2) == 1]
一般的にコンピューターで扱う整数は、ハードウェアのCPU性能に依存しており、64ビットCPUのコンピューターでは符号なしで$2^{64}$までという制約があります。つまり上記のmodの引数である2^(n-1)はコンピューターで扱うにはかなり大きな整数になりますが、HaskellのInterger型は多倍長整数を扱えるため、オーバーフロー等の制約を気にしないで、数学的な定義どおり式で表せば動作可能なプログラムとなります。
10桁の最小素数での計算（例）
計算がエラーなく処理されるかを試してみます。次のような値での処理は、さほど待たずに現実的な処理時間で出力が得られます。但しこの値が正しいかを確認するための参照値がなく若干不安ではありますが...
ghci> a=1000000007
ghci> mod (2^(a-1)) (a^2)
489841860428893000
ヴィーフェリッヒ素数を出力する
ヴィーフェリッヒ素数列wは、次のようにすればインタラクティブに表示できます。
ghci> w
[1093,3511
素数列の生成は効率的なプログラムではありませんので、2022年12月現在に調べられている範囲（$2^{64}$（約$18×10^{18}$））まで素数列を進めて、発見されていない３番めのヴィーフェリッヒ素数を発見するような範囲まで処理を進めるには、途方も無い時間を要することになるでしょう。ですから、3511が表示されたらCtrl-Cで処理を止めましょう。
現時点でヴィーフェリッヒ素数は２つしか確認されていませんので、次のように指定すれば、最初の２つを表示するようにでき、このようにした方がよいでしょう。
ghci> take 2 w
[1093,3511]
いずれにしても[1093,3511]が得られています。念のため、この結果[1093,3511]のそれぞれに対して、$2^{p-1} \equiv 1 \pmod{p^{2}}$を確認しましょう。
ghci> 2^1092
53058536290991634737396540170283858539198978395771271474551549598742698935712572215646062825029533563666321339663466341699370021653261445199826360714064955944866263669556212233526813063143284104557957658305559283253118889734882642786544086063293282737405311454375777285526890640894984855797452638426498665888535738081213096656896
ghci> mod (2^1092) (1093^2)
1
ghci> 2^3510
4123678330405087685498374221366984954089317057900104262462784875123790425105688284686404828183633224942349909964891086086676925759584441340507088115762723035704344160664907942351632995309234873230720724730672756917176305067668035542444896101177345052548810184785328926774725935717835798651444434276670015032998880347558303319286200881282911325381616533722737115450145930916961537131913757215277371800137275864121652932255403611649883620607541869985422611380925136858353400304606133084255357430301930613665490583193975084451066213800348446463658847351044167257072998067401824974454641569423021718445571724956556610925411248121271018867103767845126877315661086297321467468234450388231735471475672432407097826847947552577946732477743671085802532103177658853957960953937628759295649081969101068141306491454515657017075167637032474324430969643933568914915807923425951587598028215460799776130796600200956964482108313216473786634870222110938349411949520063457665149056813193039601852026593109271569256781038420076239187040954693149511676801395787300912507360641024
ghci> mod (2^3510) (3511^2)
1
この通り、[1093,3511]がヴィーフェリッヒ素数の条件を満たしていることが分かります。
おわりに
ご一読いただきまして有り難うございます。最近は毎週の投稿にチャレンジしていますが、少しネタ切れ感も出てきておりまして、少しショートな記事になりましたがご容赦ください。本稿は数学的な話題には触れていませんが、先の@tsujimotterさんやせきゅーんさんのブログ3にて関連した話題が詳しく解説されています。興味が沸いた方は、それらの情報を参考にすることをオススメします。
最後に、本稿を記載するために検証したHaskell環境を記しておきます。お手元の環境で検証する際に動作が異なる場合などは、参考になるかもしれません。
本稿の環境
本稿のために使用した環境は以下となります。
macOS: Sonoma 14.4 (chip: Apple M1)
GHCup: 0.1.22.0
GHC: 9.6.4
（●）（●） Happy Hacking!
/"" ＿＿""＼