西暦13446204年には素数がない

この記事はなに
用語の約束
経緯
適当な数の総積による構成
中国剰余定理
- 全ての日に制約（367桁）
- 関係式を減らす + 乱択（~9桁）
0から探す（8桁：13446204）
なぜこんなに小さいか？

この記事はなに

グレゴリオ暦の年月日を十進法で表現された数として読むと素数になることがある（例：2022年1月3日→20220103は素数）。グレゴリオ暦が未来永劫に有効であるとして、このような素数が含まれない年は13446204年に初めて訪れる。本記事では、この値を最終的に発見するまでに試みた探索の過程を説明する。

コードは
github.com
に上げています。

バージョンは、

Rust: 1.55.0
Python: 3.8.10

用語の約束

素数日
年月日を十進法で読むと素数になる日。上述の2022年1月3日など。各年の素数日の数の列が A296008 - OEIS に登録されている。

合成数日
年月日を十進法で読むと合成数になる日。2021年12月28日は明らかに合成数日である。

無素数年
素数日を含まない年（その年の全ての日が合成数日である年）。13446204は最初の無素数年である。

「グレゴリオ暦」
歴史的にはグレゴリオ暦は1582年に導入されている。また、日付と時刻に関する国際規格であるISO8601では9999年以降の扱いに関しては定めていない。本記事で「グレゴリオ暦」（西暦）という場合、西暦0年から無期限にこの暦法が適用されると仮定する。すなわち、各月に含まれる日数や閏年・閏日に関する扱いは現状に準じる*1。なお、負数（紀元前）は考えない。

経緯

@tsujimotterさんによる記事「カレンダーの上の素数〜素数には毎年出会えるか？〜」において、年月日を十進法で読んだときに素数になる日（以下「素数日」）を含む個数の期待値が1を下回る年の見積もりとして、 ${10^{154}}$ 年が与えられている。

tsujimotter.hatenablog.com

最初の目標は、具体的に無素数年を構成することでこの評価を改良することにあった。

ところが、計算を進めるうちにこの上限がどんどん下がり、愚直な探索で最初の無素数年が見つけられることが判明、ついに13446204が実際に最初の年であることが分かった*2。

結果としては最初から愚直な探索を回せばよかったものの、この宝石探しそれ自体が楽しいものだった。そのツイートを見てくださったtsujimotterさんご本人からもコメント*3をいただき、簡単にまとめることにした。

適当な数の総積による構成

具体的に無素数年を探していく。年 ${y}$ に対して、 ${10000\, y + d}$ が全て合成数なら無素数年になる。ただし、 ${d}$ はmmdd形式の月日を十進法で読んだ数で、具体的には（カレンダーを書き下すだけだが）平年（common year）なら以下の ${D_{\rm common}}$ の要素。

${ \begin{aligned} D_{\rm common} = \{ & 101, \cdots , 131, 201, \cdots , 228, 301, \cdots , 331,\\ & 401, \cdots , 430, 501, \cdots, 531, 601, \cdots, 630,\\ & 701, \cdots, 731, 801, \cdots, 831, 901, \cdots, 930, \\ & 1001, \cdots, 1031, 1101, \cdots, 1130, 1201, \cdots, 1231 \} \end{aligned} }$

閏年（leap year）なら、閏日に対応する ${229}$ を加えて、 ${D_{\rm leap}}$ の要素を取る。

${ \begin{align} D_{\rm leap} = D_{\rm common} \cup \{ 229 \} \end{align} }$

つまり、年 ${ y }$ が平年か閏年かに応じて、

${ \begin{align} Y_{t} &= \left\{ 10000\,y + d \mid d \in D_{t} \right\} & (t = {\rm common}, {\rm leap}) \end{align} }$

で定義される ${ Y_{t} }$ の要素すべてが合成数になればよい。

アプローチのひとつは、 ${y}$ が全ての ${d}$ と2以上の公約数を持つような ${y}$ を作ることである。

1231 !（3272桁）

コード

まず、無素数年の自明な例として ${1231!}$ がある。実際に計算すると、

${ 1231! \simeq 2.7 \times 10^{3271} }$

と3272桁の無素数年が得られる。こうして無素数年が存在することが示された。

mmddの総積（1001桁）

コード

階乗による構成は明らかに無駄を含んでいる。2月30日は存在しない*4 から、積を取るとき 230 は除外してよい。

平年のmmddすべての積を

${ y = \prod_{d \in D_{\rm common}} d \simeq 8.2 \times 10^{1000} }$

によってとると、1001桁の無素数年になる。なお、平年を取ったのは、こうして得られる数が明らかに400で割り切れるから。

2か5の倍数を除いたmmddの総積（412桁）

コード

20211230, 20220102, 20220105, ...は明らかに合成数。 ${d}$ が 2 または 5 を約数に持てば、 ${10000\,y+d}$ もそれを約数に持つため。だから、積から 2 か 5 の倍数は除いてもいい。

${ E = \{ d \in D_{\rm common} \mid {\rm gcd}(d,10) = 1 \} }$

として ${E}$ の総積を取ると、

${ y = \prod_{d \in E} d \simeq 2.3 \times 10^{411} }$

と412桁の無素数年が得られる。

1より大きい最小の公約数のみ共有（167桁）

コード

積を取る中で、 ${d}$ それ自体を掛ける必要はない。たとえば、 ${d=119=7 \cdot 17}$ は、 ${y}$ との公約数として7だけを共有していればよい。

また、順に掛けていくと、 ${309=3\cdot 103}$ の約数3は、 ${111 = 3\cdot 37}$ が現れた時点で積に取り込まれているから追加する必要がない。

この方法で計算すると、

${ y \simeq 1.0 \times 10^{166} }$

が得られる。

なお、構成法から明らかなように、この ${y}$ は無平方数である。

乱択アルゴリズム（~12桁）

コード

律儀にすべての ${d}$ との1より大きい公約数が存在することを保証する必要はない。適当な割合の月日に対してのみ ${y}$ と mmdd 部分が公約数を持つように条件を課し、残りの月日に対しては素数判定（ミラー・ラビンテスト）を行って、もし素数日が存在しなければ無素数年として採択すればよい。「適当な割合」の取り方によるが、0.08程度（約ひと月）あれば無素数年がある程度効率よく得られるようになる。この方法で実際に

${ y = 8333538029817 \simeq 8.3 \times 10^{12} }$

という無素数年が得られた。かなり小さくなってきた。

中国剰余定理

ここで少しアプローチを変える。 ${10000\,y+d}$ がある素数 ${p}$ で割り切れるようにするために ${p}$ それ自体を掛ける方法では、積がどんどん膨らむためなかなか小さな値は得られない。もっと直接的に、 ${10000\,y+d \equiv 0 \mod p}$ という条件を課すほうが小さい値が得られる見込みは高いだろう。

全ての日に制約（367桁）

コード

各年月日がある（小さい）素数で割り切れるという制約を課して、その条件を満たす整数を求める。具体的には、

${ \begin{align} 10000\,y + 101 &\equiv 0 \mod 3, \\ 10000\,y + 103 &\equiv 0 \mod 7, \\ 10000\,y + 107 &\equiv 0 \mod 11, \\ &\vdots\\ 10000\,y + 1231&\equiv 0 \mod 881 \end{align} }$

を解く。これを満たす最小の正整数 ${y_0}$ は、 ${q = 3\cdot 7\cdot 11\cdots 881 \simeq 5.0 \times 10^{366} }$ 未満に唯一存在することが中国剰余定理から保証されている。

なお、閏日は条件から除いておく。もし閏日（ ${10000\,y_0+229}$ ）が素数になってしまった場合、平年の素数日として ${y_0+q}$ を取れる。

さて、実際に計算すると、

${ \begin{align} y_0 = &295382390962668722563792435480\\ &899792605764414366208510569992\\ &787810541575804048893785866617\\ &574319670744354100753850262831\\ &054931517126474879283172086016\\ &931479238343604231305828997867\\ &194049155030897841192783906031\\ &461535042719576219528156266868\\ &338484803671561776298485263002\\ &290983623162932036881744610814\\ &340176502317363191226783167702\\ &596512017606064495502847995997\\ &4119552\\ \simeq&3.0 \times 10^{366} \end{align} }$

と定まる。明らかに閏年だが、幸いにも ${10000y + 229}$ は素数である。

条件をシャッフルして

${ \begin{align} 10000\,y + 101 &\equiv 0 \mod 461, \\ 10000\,y + 103 &\equiv 0 \mod 101, \\ 10000\,y + 107 &\equiv 0 \mod 43, \\ &\vdots\\ 10000\,y + 1231&\equiv 0 \mod 19 \end{align} }$

などとすれば ${y_0}$ は変わるが、 ${q}$ を上限として ${10^{366}}$ オーダーに留まるため、あまり改善しない。

関係式を減らす + 乱択（~9桁）

コード

${10000\,y+101}$ が3の倍数なら ${10000\,y+107}$ は自動的に3の倍数になる。だから、 ${10000\,y+101}$ が3の倍数になるように条件を課したなら、 ${10000\,y+107}$ に関する条件は不要になる。

さらに、総積による方法で用いたのと同様、すべての月日に対して条件を課す必要はない。適当な閾値を定めて、解いた後でそれが無素数年であれば採択、そうでなければ再試行すればよい。

肝は「適当な閾値」の取り方にある。条件に課す素数の範囲としてもよいし、条件が課されない月日の数にとってもよい。

一例として、 ${\mod 29}$ までの制約を課すとすると、100万回程度の試行で

${ \begin{align} 10000\,y+221 &\equiv 0 \mod 3,\\ 10000\,y+1003 &\equiv 0 \mod 7,\\ 10000\,y+601 &\equiv 0 \mod 11,\\ 10000\,y+813 &\equiv 0 \mod 13,\\ 10000\,y+817 &\equiv 0 \mod 17,\\ 10000\,y+327 &\equiv 0 \mod 19,\\ 10000\,y+1027 &\equiv 0 \mod 23,\\ 10000\,y+619 &\equiv 0 \mod 29\\ \end{align} }$

という連立方程式から

${ \begin{align} y = 405956995 \simeq 4.1 \times 10^8 \end{align} }$

が見つかった。西暦4億年。太陽が赤色巨星になるより前だ。

ちなみに後で行った計算で、この405956995は122番目の無素数年であることが分かった。かなりいい線を行っている。

10の8乗？　このあたりで、無素数年がずいぶん小さい領域に存在することに気付く。

ここまでの目的は構成的に無素数年を発見することで、最初の無素数年の上限を改良することにあった。なるべく小さい無素数年を探すゲーム。ところが、10の8乗オーダーの無素数年が見つかったとなると問題は違った局面に入る。ひょっとすると、0年から順に探せば最初の無素数年が見つかるのでは？

0から探す（8桁：13446204）

コード

西暦0年から探していくだけ。特に工夫もない。

手元の1.8GHzのマシンで1秒に10万年調べられる（！）ので、実は2分程度で見つかる。10時間かかるとかならともかく……最初からこれを試すべきだったのでは？（後からでは何とでも言える）

コードと……

use primeless_year::*;
use rug::Integer;

fn main() {
    for y in 0.. {
        let y = Integer::from(y);
        if y.is_probably_primeless_year() == IsPrimelessYear::Yes {
            println!("{}", y);
            break;
        }
    }
}

出力：

13446204

こうして西暦1344万6204年が最初の無素数年であることが分かった。閏年だ。

ここまで小さいと愚直な素因数分解も簡単にできて、実際にすべての月日が合成数日であることが確かめられる。しかも64ビット整数の範囲内で。

西暦13446204年の各月日の素因数分解

この中で、

${ \begin{align} 134462040331&=328343\cdot409517,\\ 134462040403&=240197\cdot559799,\\ 134462040529&=171553\cdot783793,\\ 134462040731&=257717\cdot521743,\\ \end{align} }$

あたりは2つの素因数がどちらも6桁でいい感じだ。ちなみに ${13446204}$ 自身は

${ \begin{align} 13446204 &= 2^2\cdot 3\cdot 1120517 \end{align} }$

と素因数分解される。

なぜこんなに小さいか？

13446204は小さい。こういう類の問題で、whileループで網羅的に調べられるというあまりにも「ちょうどよい」範囲に答えがあるとは予想していなかった。

無素数年はこの後にも

${ \begin{gather} 13446204,\\ 27789755,\\ 30730255,\\ 36076657,\\ 40882247,\\ \vdots \end{gather} }$

10億以下の無素数年のリスト（485個ある）

とおおよそ1000万年に1回のペースで続き、決して最初の無素数年だけが例外的に小さいわけではない。

答えが見つかったいま、新しい問題が付け加わる。「なぜこんなに小さいか？」

OEISに「n個以上の合成数が連続する区間の最小の数」の列が登録されている。たとえば、合成数が6個以上連続する最初の区間は90,91,92,93,94,95で、その次の96も合成数だから、 ${a_6=a_7=90}$ 。

Smallest start for a run of at least n composite numbers.
A030296 - OEIS

リストに登録されているように、 ${a_{365}=10726904660\sim 10^{10}}$ 。最小の無素数年 ${y_{\rm min} =13446204}$ に対して、 ${10000\,y_{\rm min} \sim 10^{11}}$ だから、ほぼオーダーは同じ。素数定理で大雑把には見積もれるはずだがちょっと自信がない。今後の課題。