連続するSmooth numberをStørmerの定理で求める

Smooth numberは面白い概念なのですがWikipediaが日本語に翻訳されていませんね。Størmerの定理に至ってはほとんど記事が見当たらなかったので詳細に説明してみました。

Smooth numberとは、

Smooth number (Wikipedia)によると、

数論ではN-smooth numberとはすべての素因数がN以下の整数の事を言う。例えば7-smooth numberは素因数が最大7の数のことである。すなわち$49=7^2$、$15750=2 \times 3^2 \times 5^3 \times 7$は7-smoothであり、$ 702 = 2 \times 3^3 \times 13$ は7-smoothではない。

50以下の5-smooth numberは何個ある？

50を超えない$2^x \times 3^y \times 5^z$をすべてリストすると、24個あるのが分かります。もちろんこの数はmaxNが大きくなればそれに伴って増えて行きます。

maxN = 50
acc = []
def smooth(idx, v, primes):
  if idx == len(primes):
    acc.append(v)
    return
  p = primes[idx]
  while v <= maxN:
    smooth(idx+1, v, primes)
    v *= p

smooth(0, 1, [2,3,5])
print(len(acc), sorted(acc))
#24, [1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 50]

連続する5-smooth number (5-smooth pair)

このなかで例えば(2,3)のように(n,n+1)が両方とも5-smooth numberであるペアを5-smooth pairと呼びます。これは$n<50$で9個あります。

sn = [(n,n+1) for n in vs if n+1 in acc]
print(len(sn), sn)
# 9  [(1, 2), (2, 3), (3, 4), (4, 5), (5, 6), (8, 9), (9, 10), (15, 16), (24, 25)]

5-smooth pairは無限にある？

それではmaxNを大きく行ったときに、さらに5-smooth pairが見つかるか試してみます。しかしmaxN=1000にしてみても(80,81)が最後で計10個です。

maxN = 1000
:
print(len(sn), sn)
# 10  [(1, 2), (2, 3), (3, 4), (4, 5), (5, 6), (8, 9), (9, 10), (15, 16), (24, 25), (80, 81)]

OEISのA002071に素数の個数kの時のsmooth pairの数があり、k=3（5-smooth）の時は10となり、それ以上ないようです。

上限はありそうですが、それがいくつか分からないとすべて見つけたかどうかは分かりません。

Størmerの定理

そこでStørmerの定理の出番となります。Størmer's theorem (Wikipedia)によると、

smooth pairの上限を与え、すべてのN-smooth pairをペル方程式をつかって見つける方法を提供する

とあり、まさにこの目的を達成できそうにです。Wikipediaに詳しい手順が書いてありますが、要約すると以下のようになります。

Størmerの定理を用いてすべてのN-smooth pairを見つける手順

1. 与えられた最大素数nに対してそれ以下の素数をリストする。
2. それぞれの素数を2個以上使わないで掛け合わせた数 (squarefree number) qをすべてリストして、2を除く。
3. そのqに対してペル方程式$x^2-2qy^2=1$を解く（解は複数ある）
4. その解の$x$に対して$(x-1)/2$,$(x+1)/2$がペアの候補となるのでこれらが両方ともn-smooth numberかチェックする。
5. 3,4をすべてのqに対して行う

以下に$n=5$の時の手順を詳細に説明します。

1. 素数のリスト

最大が５なので$ [ 2,3,5 ] $です。この個数を$k$とします（この場合は３）
pythonではsympyのprimerangeを使えばいいですね。

2. 素数を2個以上使わないで掛け合わせた数のリスト（２を除く）

$ [ 1,3,5, 6, 10, 15, 30] $の7個です。それそれの素数の指数は0か1なの個数は$2^k-1$になります。この要素をqとします。

3. qに対してペル方程式を解く

\begin{align}
&q = 1の時\\
&x^2-2y^2=1を解く \\ 
&最初の3つの解は\\
&(x,y)=(3,2),(17,12),(99, 70),...\\
\end{align}

4. smooth pairの候補をチェック

\begin{align}
&(n,n+1)=((x-1)/2,(x+1)/2) \\
&=(1,2),(8,9),(49,50),... \\
&49は5-smoothではないので\\
&答えは(1,2),(8,9)となる
\end{align}

5．3,4をすべてのqに対して行う

3,4をすべてのqに対して行うことによって、以下のようにすべてのsmooth number ペアを求めることが出来ました。

q	ペル方程式	(n,n+1)
1	$x^2-2y^2=1$	(1,2),(8,9) ,(49,50),...
3	$x^2-6y^2=1$	(2,3),(24,25) ,(242,243),...
5	$x^2-10y^2=1$	(9,10), ...
6	$x^2-12y^2=1$	(3,4), ...
10	$x^2-20y^2=1$	(4,5),(80,81), ...
15	$x^2-30y^2=1$	(5,6), ...
30	$x^2-60y^2=1$	(15,16), ...

ペル方程式の解をいくつまで調べれば良いか？

ただこれでも、一つのqに対して解が2つ見つかる場合もあるし、1つの事もあるので、ペル方程式の解を何番目までチェックすればよいか分からないと無限の時間がかかってしまいます。そこは幸いにも以下のように上限が示されていて、$N=5$の時は$3$になるので3つめまで調べれば良いことになります。

ペル方程式の解は max(3, (N + 1)/2) 番目まで調べれば良い

ペル方程式を解くコード

pythonのsympyには一般ペル方程式を解くdiop_DNという関数が用意されていますがあまり速くないので以下の記事の中で紹介した連分数からペル方程式の解を発生させるcfp2pell/cfsqrtを使います。

【Project Euler】Problem 94: 正三角形に近いヘロンの三角形

ペル方程式$x^2-2y^2=1$を解くコードはこんな感じです。

q = 1
ansPelEq = cfp2pell(cfsqrt(2*q))
for i in range(1,3+1):
  (x, y) = next(ansPelEq)
  n = (x-1)//2
  print(i,":", (x,y),"->", (n,n+1))
# 1 : (3, 2) -> (1, 2)
# 2 : (17, 12) -> (8, 9)
# 3 : (99, 70) -> (49, 50)

Størmerの定理を用いたコード

これですべての手順はわかったのでこれをコードにします。

def isSmooth(n, prims):
  for p in prims:
    while n % p ==0: n //= p
  return n == 1

N = 5
prims = list(sp.primerange(2,N+1)) #----- prime list p <= N
#----- square free numbers
pset = [1]
for p in prims:
  pset += [sub * p for sub in pset]
pset.remove(2) # remove 2
print(prims, pset)
#------ M : maximum iteration
k = len(prims)
M = max(3, (prims[-1] + 1)//2)
print(f"N = {N}, k = {k}: max_i = {M}")
cnt, maxn = 0, 0
ans = []
for q in pset:
  #print(f"------ q = {q}-------")
  ansPelEq = cfp2pell(cfsqrt(2*q)) # solve pell equation
  for i in range(M+1):
    (x, y) = next(ansPelEq)
    n = (x-1)//2
    if isSmooth(n, prims) and isSmooth(n+1, prims):
      cnt += 1
      if n > maxn: maxn = n
      ans += [(n, n+1)]
    elif i == 0: break
print(f"{N}-smooth number pairs : {sorted(ans)}")
# [2, 3, 5] [1, 3, 6, 5, 10, 15, 30]
N = 5, k = 3: max_i = 3
5-smooth number pairs : [(1, 2), (2, 3), (3, 4), (4, 5), (5, 6), (8, 9), (9, 10), (15, 16), (24, 25), (80, 81)]

N-smooth pairの数と最大値

このコードを使って連N-smooth pairの数（OEIS A002071）と最大値（OEIS A002072）をN=53まで求めてみました。

N	ペアの数	nの最大値
2	1	1
3	4	8
5	10	80
7	23	4374
11	40	9800
13	68	123200
17	108	336140
19	167	11859210
23	241	11859210
29	345	177182720
31	482	1611308699
37	653	3463199999
41	869	63927525375
43	1153	421138799639
47	1502	1109496723125
53	1930	1453579866024

（開発環境：Google Colab）

この考え方はProject Euler、Problem 581: 47-smooth Triangular Numbersを解くのに役に立ちます。