８．Wieferich素数 (Fermat quotient)

定義

任意の素数pと、pと互いに素な整数aに対し、

a^p-1≡1 (mod p)

が成り立つ。これをフェルマーの定理（または、フェルマーの小定理）と呼ぶ。
この式より、a^p-1-1はpで割り切れることがわかり、
q=(a^p-1-1)/p をフェルマー商 (Fermat quotient) と呼ぶ。

では、この q がさらに p で割り切れる、
つまり、a^p-1≡1 (mod p²) が成り立つことがあるだろうか？

Wieferich は、1909年に

x^p+ y^p = z^p (x, y, z, p∈N)

が自然数解を持つならば（これをフェルマーの大定理の第１の場合と呼ぶ）

2^p-1-1≡0 (mod p²)

が成り立つことを証明した。
これにちなみ、a^p-1-1が p²で割り切れる場合、pをWieferich素数と呼ぶ。

Mirimanoff は、フェルマーの大定理の第１の場合が解を持つならば、

2^p-1-1≡0 (mod p²)

と同時に、

3^p-1-1≡0 (mod p²)

も成り立つことを示した。
さらに、A. Granville と M. Monagan は、a=5, 7, 11, ……の場合も同時に成り立つことを示した。

探査プログラム

フェルマーの大定理は既に証明されてしまったが、 Wieferich素数の探査は、依然として問題として残っている。
ここでは、a=2〜999（ただしベキ乗数を除く）について、
a^p-1-1≡0 (mod p²) を満たすような素数 p を探してみることにする。

プログラムは以下のとおり。

 10   print=print+"fermat.txt"
 20   ' fermat.ub
 30   dim T%(999)
 40   for I%=2 to 999:T%(I%)=0:next I%
 50   for I%=2 to 31:T%(I%^2)=1:next I%
 60   for I%=2 to 9:T%(I%^3)=1:next I%
 70   for I%=2 to 5:T%(I%^4)=1:next I%
 80   T%(3^5)=1:T%(2^5)=1:T%(2^7)=1:T%(2^9)=1
 90   J%=0
100   for I%=2 to 999
110     if T%(I%)=0 then inc J%:T%(J%)=I%
120   next I%:' total number is 959
130   ' case p=2
140   for I%=1 to J%
150     if T%(I%)@4=1 then print T%(I%);":";P
160   next I%
170   ' case p>2
180   P=3
190   Pm=(P-1)\2:P2=P^2:P2m=P2-1
200     for I%=1 to J%:A=modpow(T%(I%),Pm,P2)
210       if or{A=1,A=P2m} then print T%(I%);":";P
220     next I%
230   if P<N then P=nxtprm(P):goto 190

40〜80行は、ベキ乗の場合を対象外とするための処理。
100〜120行で、探査対象を配列に格納する。1000以下の場合、総数は959個となる。
130行以降が条件 a^p-1-1≡0 (mod p²) をチェックする部分だが、
少し小技を使っている。

3以上の素数は奇数なので p-1 は偶数となる。
(p-1)乗で1になるということは、(p-1)/2乗で 1 または -1 (= p^2-1) になっているはずである。
(p-1)乗と(p-1)/2乗では、binary method で掛け算が１回減り、
その分、比較の回数が１回増えるが、トータルでは多少早くなるはずである。
（実際にやってみると、数％速度が向上した）

それで、180行〜230行はpが3以上の素数の場合、140行〜160行はp=2の場合となっている。
p=2の場合は、１乗なので、aを4で割って1余るかどうか、直接調べればよい。
p>3の場合は、200行がmodpow計算、210行が判定となっている。
modpowで (p-1)/2, p^2 と、判定 p^2-1 という値が必要になるが、この計算結果を生かしたいので、
全体を p に関するループにし、内側のループは探査対象 a のループとする。
次の素数の生成は ubasic の組込関数 nxtprm(p) を使う。

ただし nxtprm(p) は p<=2^32-5 までしか使えないので、もっと広いレンジ
(例えば 10^10 まで) をターゲットとしたとき、別な手段で素数を生成する必要がある。
オーソドックスな手段としては、最初の n個の素数の積 #Pn をとり、
1 <= a < #Pn の範囲で gcd(#Pn, a)>1 となる数をふるい落して、
#Pn と互いに素となる a を配列に格納しておき、

#Pn + a1, #Pn + a2, ...

を次の素数の候補とする方法がある。
13までの積 #Pn=2*3*5*7*11*13=30030 の場合、互いに素となる数は5760個。
17までの積 510510 の場合、互いに素となる数は92160個。
ubasicではこの個数の配列を取れないので、13までの積のふるいを使う。
プログラムは、120〜130行の間に、

122   L%=2*3*5*7*11*13:dim Diff%(L%):C%=0:' L=30030
124   for I%=1 to L%-1
126     if gcd(I%,L%)=1 then inc C%:Diff%(C%)=I%
128   next I%:' c%=5760

を挿入し、メインは 230行以降に、

240   Base=((2^32-5)\L%)*L%:N=10^11
250   for K%=1 to C%:P=Base+Diff%(K%):if prmdiv(P)<>P then 300
260     Pm=(P-1)\2:P2=P^2:P2m=P2-1
270     for I%=1 to J%:A=modpow(T%(I%),Pm,P2)
280       if or{A=1,A=P2m} then print T%(I%);":";P
290     next I%
300   next K%
310   Base=Base+L%:if Base<N then 250

を追加する。
p の初期値 Base は L%=#Pn の整数倍から始め、Nをリミットとする。
p は素数の候補であり、必ず素数になる訳ではないので、
270〜290行が空回りにならないように、250行にpの素数判定を入れている。
(素数判定を prmdiv(p) ではなくフェルマーテストにすると速くなりそうだが、
  効果の程は実験してみないと不明。実際に2^32越えを探査範囲とするとき、チェックすることにする。）

実行結果

Wilfrid Keller と Jörg Richstein は、

100 < a < 1000 の a について、10⁴ < p < 10¹¹ の範囲まで
a < 100 の a について、2³² < p の範囲

について探査を行い、a=5 に対する解 p=188748146801 を見つけた。

[参考文献]

Wilfrid Keller, Jörg Richstein,
Solutions of the congruence a^p-1≡1 (mod p^r),
Mathematics of Computation, 74(2005), 927-936.

このページでは、現在、Wilfrid Keller, Jörg Richstein の探査範囲を突破し、
1.72 * 10¹²まで計算が進んでいる。 結果はこちら。
8.6×10⁹ 以降は、Richard Fischer (February 23, 2009) の探査結果に基づく。

Wanted list (p ≤ 2.31 * 10¹²の素数については調査済み)

解が見つかっていないのは以下の30個

 47,  72, 186, 187, 200, 203, 222, 231, 304, 311,
335, 347, 355, 435, 454, 542, 546, 610, 639, 662,
760, 772, 798, 808, 812, 858, 860, 871, 983, 986

これらの解は？

以下の数は、a < p となる解が見つかっていない。

21, 28, 29, 50, 51, 61, 73, 82, 89, 99,

116, 126, 129, 132, 145, 154, 181, 188,

202, 214, 218, 226, 229, 236, 237, 249, 259, 261,
268, 270, 276, 293,

301, 306, 309, 327, 330, 334, 338, 341, 350, 351,
356, 357, 360, 373, 377, 379, 386, 388, 394, 399,

402, 412, 414, 417, 420, 437, 458, 466, 469, 478,
482, 489, 496,

506, 541, 543, 567, 570, 582, 593,

601, 604, 605, 609, 612, 615, 623, 628, 629, 630,
640, 654, 665, 668, 673, 674, 677, 690, 692, 694,
696, 697, 698,

701, 711, 714, 718, 721, 724, 725, 727, 730, 733,
743, 751, 753, 759, 769, 775, 777,

800, 807, 811, 829, 834, 836, 837, 838, 842, 843,
845, 846, 851, 859, 867, 869, 878, 883, 898,

905, 912, 914, 919, 920, 931, 942, 944, 946, 954,
972, 978, 980,

これらについて a < p を満たす解は存在するか？