素数41は強い？

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

はじめに

突然ですが、わたし が好きな素数は $41$ です。
今月は $41$ という数字を見る機会が多い月でした。乃木坂46 の次のシングルが 41枚目。@tsujimotter さんが、日曜数学会ミニin北海道で「$\sqrt{41}$の概算」という魅力的なタイトルでご発表。参加していない者としてはブログ記事を楽しみにしています。

$41$ はオイラーの幸運数（Lucky Number）¹のひとつとされており、$n^2+n+41$ はオイラー素数²を生成します。すなわち、$n^2+n+41$ は$0 \le n \le 39$の範囲で連続して$40$個の素数を生成する式として知られています。これほどの多くの素数を生み出すには数学的な背景があるようですが、本稿ではその論理的な部分には深入りせずに、$n^2+n+p$（$p$は素数．$0 \le n \le p-2 $）という式にどれくらいの素数が含まれているのかについてほかの素数$p$について調べてみましょう。

最強の p？

$f_p(n)=n^2+n+p$ がどれだけ素数を生むかを実験的に調べて 「最強の $p$」 を探索してみましょう。ここでscore(p)を次のように定義します。

score(p) = $0\le n\le p-2$ で $f_p(n)$ が素数になる個数

このscore(p)がp-1のどれくらいの割合を占めるかという $\mathbf{ratio(p) = score(p)/(p-1)}$ が大きいほど「強い $p$」とします。次のプログラムで「p,score,ratio」をCSV形式で出力することができます。

Main.hsのソースコード

Main.hs

-- Usage:
--   ghc -O2 -threaded Main.hs -package containers
--   ./Main 100000 --csv out.csv   +RTS -N
--
-- Output CSV columns:
--   p,score,ratio

module Main (main) where

import Data.Char (isDigit)
import Data.Bits
import Data.Word
import Numeric (showFFloat)
import System.Environment (getArgs)
import System.IO (withFile, IOMode(WriteMode), hPutStrLn)

import qualified Data.Set as S

--------------------------------------------------------------------------------
-- Segmented sieve
--------------------------------------------------------------------------------
primes :: [Integer]
n##x@(m:p:y)=[n|gcd m n<2]++(n+2)##last(x:[m*p:y|p^2-3<n])
primes=2:3:5##primes

primesUpTo :: Integer -> ([Word64], S.Set Word64)
primesUpTo limit = (psW, S.fromList psW)
  where psI = takeWhile (<= limit) primes
        psW = map fromIntegral psI

--------------------------------------------------------------------------------
-- Miller–Rabin
--------------------------------------------------------------------------------
evenW :: Word64 -> Bool
evenW x = (x .&. 1) == 0
oddW :: Word64 -> Bool
oddW x = (x .&. 1) == 1
mulMod :: Word64 -> Word64 -> Word64 -> Word64
mulMod !x !y !m =
  fromInteger ((toInteger x * toInteger y) `mod` toInteger m)

mr :: Word64 -> Bool
mr n
  | n < 2     = False
  | n == 2    = True
  | evenW n   = False
  | otherwise = all (\a -> a `mod` n == 0 || witness a)
                [2,325,9375,28178,450775,9780504,1795265022]
  where
    -- (s,d) such that n-1 = d * 2^s with d odd
    (s,d) = until (oddW . snd)
                  (\(k,x) -> (k+1, x `div` 2))
                  (0, n-1)
    powMod :: Word64 -> Word64 -> Word64 -> Word64
    powMod a e m = go (a `mod` m) e 1
      where
        go !_ 0 !r = r
        go !b !e' !r =
          let !r' = if oddW e' then mulMod r b m else r
              !b' = mulMod b b m
          in go b' (e' `div` 2) r'
    witness :: Word64 -> Bool
    witness a =
      let !x0 = powMod a d n
          loop 0 !x = x == n-1
          loop k !x =
            x == n-1 || loop (k-1) (mulMod x x n)
      in x0 == 1 || loop (s-1) x0

isPrime :: S.Set Word64 -> Word64 -> Word64 -> Bool
isPrime !primeSet !limit !n
  | n <= limit = S.member n primeSet
  | otherwise  = mr n

--------------------------------------------------------------------------------
-- Scoring
--------------------------------------------------------------------------------
scoreP :: S.Set Word64 -> Word64 -> Word64 -> Word64
scoreP !primeSet !limit !p
  | p == 2 = 1
  | otherwise =
      let !maxN = p - 2
          go !n !f !acc
            | n > maxN  = acc
            | otherwise =
                let !acc' = if isPrime primeSet limit f then acc + 1 else acc
                    !f'   = f + (2*n + 2)
                in go (n + 1) f' acc'
      in go 0 p 0

--------------------------------------------------------------------------------
-- Main + CSV
--------------------------------------------------------------------------------
data Mode = PrintTop | WriteCSV FilePath deriving (Eq, Show)
parseArgs :: [String] -> (Integer, Mode)
parseArgs args =
  let limit = case [a | a <- args, all isDigit a] of
                (x:_) -> read x
                _     -> 100000
      mode  = case dropWhile (/= "--csv") args of
                (_:fp:_) -> WriteCSV fp
                _        -> PrintTop
  in (limit, mode)

main :: IO ()
main = do
  args <- getArgs
  let (limitI, mode)  = parseArgs args
      limitW          = fromIntegral limitI :: Word64
      (ps, primeSet)  = primesUpTo limitI

  case mode of
    WriteCSV fp -> do
      putStrLn $ "Writing CSV up to p <= " ++ show limitI
      withFile fp WriteMode $ \h -> do
        hPutStrLn h "p,score,ratio"
        let emit p = do
              let s      = scoreP primeSet limitW p
                  ratioD = (fromIntegral s :: Double) / fromIntegral (p - 1)
                  ratio  = showFFloat Nothing ratioD ""
              hPutStrLn h (show p ++ "," ++ show s ++ "," ++ show ratio)
        mapM_ emit ps
      putStrLn $ "Wrote CSV: " ++ fp

    PrintTop -> do
      putStrLn "Run with:"
      putStrLn "  ./Main 100000 --csv out.csv   +RTS -N"

このプログラムで$100,000$まで出力したCSVをグラフにしてみました。ratioのグラフでは$p$が大きくなると0.4より小さくなっているようです。どうもある程度大きい$p$に対しては上界があるように推測できます。

$41$は$0$付近で見えづらくなってしまっていますが$1.0$になるため「強い」ということにはなりそうです。同様に$1.0$になるのは「オイラーの幸運数¹」です。次点は、$101$ で $score(p)=68, ratio(p)=0.68$ です。$ratio(p) \ge 0.4$ となる $p$ は、$10,000$ まではいくつかありますが、$10,000 \lt p \lt 100,000$ では確認できませんでした。$10,000$ を超えたなかで強い$p$（Top10）を一覧表にしてみました。

p	score	ratio
27941	11150	0.399069434502505
21377	8398	0.392870508982036
12107	4653	0.384354865356022
22697	8695	0.383107155445894
11777	4486	0.380944293478261
55661	21155	0.380075458138699
19421	7317	0.376776519052523
14627	5434	0.371530151784493
10847	3985	0.367416559100129
10457	3829	0.366201224177506

ブニャコフスキー予想とベイトマン・ホーン予想

割合は$0.4$以上にはならないにしても$0$に近づくようには感じられません。$p$を無限にしてもある割合で素数が生成できる式となるのでしょうか。それとも割合は緩やかに$0$に近づき、$f_p(n)$ は素数を生成できなくなるのでしょうか。ここで次のscore(p)のグラフをみてみると、その下限は徐々に大きくなっているように推測できます（グラフの赤線）。$p$が無限になったとしても、この赤線が$0$に向かうようなことはないように感じられます。

この特徴に対する理論的な考察は、多項式が生成する素数の個数に対する研究です。そのような研究分野において、２つの予想を紹介します。

ブニャコフスキー予想

ウクライナ出身の数学者ブニャコフスキーが1857年に示した予想³があります。2026年2月現在、証明はされていない（つまり未解決である）が、反例も見つかっていない予想です。

ブニャコフスキー予想

整数係数を持つ2次以上の既約多項式は、自然数の引数に対して1より大きな最大公約数を持つ無限集合を生成するか、もしくは無限個の素数を生成するという予想

ベイトマン・ホーン予想

1962年にポール・T・ベイトマンとロジャー・A・ホーンが提案した、数論における非常に強力な未解決の予想⁴があります。整数係数の多項式の集合が、どれくらいの頻度で同時に素数を生成するかを予測する一般公式を与える予想です。

ベイトマン・ホーン予想

多項式セット $S = \{f_1(x), f_2(x), \dots, f_m(x) \}$ について $n \lt x$ の範囲で、すべての $f_i(x)$ が素数になるような $n$ の個数 $\pi _{S} (x)$ は、以下の漸近式で表される：

$\displaystyle \pi _{S} (x) \sim \int _{2}^{x} \frac{dt}{(\log{t})^m} $

• $D$(定数)：多項式 $f_i(x)$ の次数の積（$\displaystyle D = \prod _{i=1}^{m} \deg{f_i}$）
• $C$(定数)：多項式のセット $S$ に依存する定数（$S$ の素数生成効率を表す）
• $m$：多項式の個数

ここで、今回の多項式 $ n^2+n+p $ が無限に素数をとるか？という問題は、このブニャコフスキー予想の特別な場合です。また、ベイトマン・ホーン予想によると、二次多項式 $ n^2+n+p $ が素数を与える個数は、概ね

$ \displaystyle
\mathrm{score}(p)
\approx
\sum_{n=0}^{p-2}\frac{C_p}{\log f_p(n)}
$

と予想されます。ここで $C_p$ は局所密度（local factor）。
そこで粗い近似として

$ \displaystyle
\mathrm{score}(p)
\approx
\frac{C_p\,p}{\log(p^2)}
\approx
\frac{C_p\,p}{2\log p}
$

を使います。この式の局所密度は　$\displaystyle C_p=\prod_q \frac{1-\frac{N_q}{q}}{(1-\frac1q)^2}$という無限積ですが、各因子は 1 にかなり近く、大きい素数ほど 1 に急速に近づくため、全体として $C_p \approx 1$ としてもラフな近似ができそうです。今回は正確な$C_p$を計算せずに、$C_p \approx 1$ として粗い近似をグラフに重ねてみました（上述の図中の赤線）。
これで感覚的な下限の線が、ベイトマン・ホーン予想で予想されている個数とよく合うことも見てとれました。

参考）グラフの作成プログラム

グラフの作成(python)

draw-plot.py

import csv
import math
import matplotlib.pyplot as plt

ps = []
scores = []
ratios = []

with open("out.csv", newline="") as f:
    r = csv.DictReader(f)
    for row in r:
        ps.append(int(row["p"]))
        scores.append(int(row["score"]))
        ratios.append(float(row["ratio"]))

# --- Bateman–Horn ラフ近似 ---
bh_pred = []
for p in ps:
    if p > 2:
      bh_pred.append(p / (2 * math.log(p)))
    else:
      bh_pred.append(0)

# --- score vs p ---
plt.figure(figsize=(8, 5))
plt.scatter(ps, scores, s=1, marker=',', linewidth=0)   # 実測
plt.plot(ps, bh_pred, c='red', linewidth=1, label="Bateman–Horn (rough)")  # 理論
plt.margins(x=0, y=0)
plt.xlabel("p")
plt.ylabel("score(p)")
plt.title("score(p) for f_p(n)=n^2+n+p, n=0..p-2")
plt.tight_layout()
plt.savefig("score_vs_p.png", dpi=300)
plt.close()

# --- ratio vs p ---
plt.figure(figsize=(8, 5))
plt.scatter(ps, ratios, s=1, marker=',', edgecolors='none', linewidths=0)
plt.margins(x=0, y=0)
plt.xlabel("p")
plt.ylabel("score(p)/(p-1)")
plt.title("ratio of prime outputs")
plt.tight_layout()
plt.savefig("ratio_vs_p.png", dpi=300)
plt.close()

おわりに

いかがでしたでしょうか。本稿では、素数 $41$ の強さを考えることから始めて、多項式が生成する素数の個数について考えを進めてみました。そして『ベイトマン・ホーン予想』という現代数論において非常に一般性が高く、極めて難しい未解決問題について紹介することができました。
最後に、本稿を記載するために検証したHaskell環境を記しておきます。お手元の環境で検証する際に、動作が異なるときには参考になるかもしれません。

本稿の環境

本稿のために使用した環境は以下となります。 macOS: Tahoe 26.3 (chip: Apple M1) GHCup: 0.1.50.2 GHC: 9.12.2

ご一読いただきまして有り難うございます。
（●）（●） Happy Hacking!
/""　＿＿""＼

1. オイラーの幸運数 Wikipedia
   [2, 3, 5, 11, 17, 41] ↩ ↩²

2. オイラー素数 [41,43,47,53,61,71,83,97,113,131,151,173,197,223,251,281,313,347,383,421,461,503,547,593,641,691,743,797,853,911,971,1033,1097,1163,1231,1301,1373,1447,1523,1601] ↩

3. ブニャコフスキー予想 ↩

4. ベイトマン・ホーン予想(en) ↩