Ruby で百世紀までの素数日リスト

Ruby で素数日を得る記事が意外と少なかったので書いてみた。
直接のきっかけは以下の記事を見たこと。
[Pythonで「素数日」を求めたら処理が長すぎて日が暮れそうになった - Qiita](https://qiita.com/yukid/items/3bb121e151903666a37e)

# 素数日とは

西暦の年月日から以下のように整数を作ったとき，その整数が素数であればその日付を素数日と呼ぶらしい。

例）2019年3月1日　→ 20190301（素数！）

見てのとおり，月と日は 2 桁にしたうえで年と並べる。

# やること

この記事では，先述の記事のように，一定期間中の素数日をすべてリストアップしてみる。
具体的には西暦 1 年から 10000 年までの百世紀間としよう（先述の記事では西暦 0〜9999 年）。
この百世紀は 3,652,427 日あり，素数日は 216,511 個あるようだ[^count]。割合は 5.9% 強。

[^count]: この個数が間違ってたら恥ずかしいなあ。検証してね。

[素数定理](https://ja.wikipedia.org/wiki/%E7%B4%A0%E6%95%B0%E5%AE%9A%E7%90%86)によれば，100001231[^lastday] までの素数はざっと 576 万個超といったところ。割合にして 5.8% 弱。
ということは，この区間中の整数で日付と解釈できる数の場合，素数率が微かに高いことになる。理由は分からん。年・月・日の「日」に奇数が多い[^dayparity]ことが効いているのだろうか？

[^lastday]: 言うまでもなく西暦百世紀最後の日に対応する整数。
[^dayparity]: 一年間の日付の中で，「日」は奇数が偶数よりもやや多い（4% ほど）。このことは，大の月（31 日までの月）は奇数が偶数よりも 1 日多く，小の月は（うるう年の  2 月を除き）同数であることを考えれば分かる。

なお，英語は苦手なので次節以降に出てくる変数名は苦し紛れに考えたものが多いので注意してくれ。

# コード

コードは以下の形で記述することに決めた。

```rb
def prime_dates(start_year, end_year)
  # start_year 年の元旦から end_year 年の大晦日までの
  # 素数日の日付を "2019-03-01" といった形式の文字列に
  # したものの配列を返す
end

IO.write "prime_dates.txt",
  prime_dates(1, 10_000).join("\n") + "\n"
```

ファイルに書き出しておけば，リファクタリングでバグが混入するのに気付きやすい。


## コード(1)

最初にこんなコードを書いた。

```rb
require "date"
require "prime"

def prime_dates(start_year, end_year)
  result = []
  Date.new(start_year, 1, 1).upto(Date.new(end_year, 12, 31)) do |date|
    if (date.year * 10_000 + date.month * 100 + date.day).prime?
      result << date.to_s
    end
  end
  result
end

IO.write "prime_dates.txt",
  prime_dates(1, 10_000).join("\n") + "\n"
```

日付と素数なので，標準添付ライブラリーの date と prime は使わせてもらう。

やっていることは，最初の日から最後の日までの日付オブジェクト（Date オブジェクト）を生成し，素数かどうかを判定して，素数なら日付を文字列化したものを空配列に放り込んで，最後にそれを返しているだけ。

Ruby の Date オブジェクトは `upto` ができたり，範囲オブジェクトの始点・終点に使えたりするところがいいよね。

今の場合，空配列を用意してそこに `<<` して最後に配列を返す，などという七面倒なことをしなくても，`select` で絞ったのをいきなり返せばいい。
ただ，あとの節でリファクタリングしていく関係で，こういう形をとっている。

さて，私の環境で実行時間は 20.4 秒（本当は有効数字 3 桁もないけど，一応小数点以下 1 桁まで書いておく）。

うーん，けっこうかかるね。

ちなみに，最後にファイルに書き出すところではたいして時間を食っていないようなので，おおむね `prime_dates` メソッドの実行時間と思ってよいようだ。

また，プロファイラーを使った調査で，どうやら実行時間の大半は素数判定に費やされているぽいことが分かった。素数判定は自分で書くより標準添付ライブラリーを素直に使ったほうが速いだろう。

## コード(2)　年・月・日の 3 重ループ

先のコードでは Date オブジェクトからメソッド呼び出しを 3 回も使って年・月・日を取り出しているのが足を引っ張ってるのかなあ？と思った。

そこで，年・月・日の 3 重ループでやってみた。
当たり前だけど，この 3 重ループで作られる年・月・日には，「2 月 30 日」など，ありえない組み合わせがある。それは弾かなければならない。
自分で判定するのは面倒。速度面を考えても [Date.valid_date?](https://docs.ruby-lang.org/ja/2.6.0/method/Date/s/valid_civil=3f.html) メソッドを使うのがよさそうだ。


```rb
require "date"
require "prime"

def prime_dates(start_year, end_year)
  result = []
  start_year.upto(end_year) do |year|
    1.upto(12) do |month|
      1.upto(31) do |day|
        next unless Date.valid_date?(year, month, day)
        if (year * 10_000 + month * 100 + day).prime?
          result << "%04d-%02d-%02d" % [year, month, day]
        end
      end
    end
  end
  result
end


IO.write "prime_dates.txt",
  prime_dates(1, 10_000).join("\n") + "\n"
```

18.7 秒。速くなったぞ！

## コード(3)　偶数日を省く

年・月・日の「日」が偶数の場合，そこから作られた数値も当然偶数となる。
偶数の素数は 2 しかあり得ないが，日付から作られる整数に 2 は無い（西暦 0 年 0 月 2 日なんて無いからね）ので，「日」が偶数の日付は全部無視してよい。

そこで，「日」のループを

```diff
-      1.upto(31) do |day|
+      1.step(31, by: 2) do |day|
```

のように変更しよう。

こうなる：

```rb
require "date"
require "prime"

def prime_dates(start_year, end_year)
  result = []
  start_year.upto(end_year) do |year|
    1.upto(12) do |month|
      1.step(31, by: 2) do |day|
        next unless Date.valid_date?(year, month, day)
        if (year * 10_000 + month * 100 + day).prime?
          result << "%04d-%02d-%02d" % [year, month, day]
        end
      end
    end
  end
  result
end

IO.write "prime_dates.txt",
  prime_dates(1, 10_000).join("\n") + "\n"
```

時間のかかる（はずの）素数判定を半分程度の回数に減らせたわけだから，劇的に速くなるのでわっ？！

いざ実行　→ 18.0 秒
あれ？　？？？？

はー，缶ビールでも飲むか・・・

先のコードと新しいコードを何度か実行したところ，確かに誤差ではなく実際に少し速くなっているようではある。しかし，少ししか違わないのはなぜか？

おそらくこういうことだ。
`prime?` メソッドが重い（時間がかかる）とはいっても，偶数なんざ瞬殺（一瞬で素数でないと判定）されちまうんである。
判定に時間がかかるのは，小さな素因数を持たない数だけなのだ。

## コード(4)　5 の倍数日も省く

先のコードでは，劇的な効果はなかったが若干の効果はあった。
であれば，同じようにしてもう少し省ける日があるのではないか。

すぐに分かるのは，年・月・日の「日」が 5 の倍数の日付も省けること。つまり，偶数に加え，「5 日」「15 日」「25 日」も省ける。

こんなふうに変えよう：

```diff
 def prime_dates(start_year, end_year)
   result = []

+  days = [*1..31].reject{ |day| day.even? || day % 5 == 0 }
   start_year.upto(end_year) do |year|
     1.upto(12) do |month|
-      1.step(31, by: 2) do |day|
+      days.each do |day|
```

全体はこうなる：

```rb
require "date"
require "prime"

def prime_dates(start_year, end_year)
  result = []
  days = [*1..31].reject{ |day| day.even? || day % 5 == 0 }
  start_year.upto(end_year) do |year|
    1.upto(12) do |month|
      days.each do |day|
        next unless Date.valid_date?(year, month, day)
        if (year * 10_000 + month * 100 + day).prime?
          result << "%04d-%02d-%02d" % [year, month, day]
        end
      end
    end
  end
  result
end

IO.write "prime_dates.txt",
  prime_dates(1, 10_000).join("\n") + "\n"
```

あまり期待できないが，いざ実行　→ 17.7 秒。
何回か試したが，やはり若干の効果はあったようだ。

## コード(5)　3 の倍数を省く

先のコードでは，素因数 2 や 5 を含む数の判定を省くため，年・月・日の「日」が素因数 2, 5 を含むものをあらかじめ排除した。

素因数 3 を含むものも同じようにできないか？
「日」だけで判定することはできない。年・月が何であるかによるからだ。

うーむ，打つ手はないのか？
はっ！（膝を打つ）

`y` 年 `m` 月 `d` 日という日付から

```rb
y * 10_000 + m * 100 + d
```

という整数を作るよね。
もし `(y * 10_000 + m * 100) % 3` が 1 なら，`d % 3` が 2 のとき，全体は 3 で割り切れるから判定が省ける！

つまり，だ。`(y * 10_000 + m * 100) % 3` の値ごとに調べるべき `d` のリストをあらかじめ作っておけば，おおむね 2/3 くらいに減らせるんじゃないか。

そこで書いたのがこのコード：

```rb
require "date"
require "prime"

def prime_dates(start_year, end_year)
  result = []
  days = [*1..31].reject{ |day| day.even? || day % 5 == 0 }
  days_by_remainder = {}
  3.times do |i|
    days_by_remainder[i] = days.reject{ |day| (i + day) % 3 == 0 }
  end

  start_year.upto(end_year) do |year|
    y = year * 10_000
    1.upto(12) do |month|
      ym = y + month * 100
      days_by_remainder[ym % 3].each do |day|
        next unless Date.valid_date?(year, month, day)
        if (ym + day).prime?
          result << "%04d-%02d-%02d" % [year, month, day]
        end
      end
    end
  end
  result
end

IO.write "prime_dates.txt",
  prime_dates(1, 10_000).join("\n") + "\n"
```

ちょっと複雑ではある。
いざ実行　→ 17.7 秒。
何回か実行したけど，先のコードとどっちが速いか判らない。

我ながらいいアイデアだと思ったのだが，残念な結果に終わった。
確かに，3 の倍数なんざ `prime?` で瞬殺なのは分かる。
しかし，少しくらい速くなってくれてもよさそうなものだ。
ややこしい処理で `days_by_remainder` を作るのは，メソッドの冒頭で 1 回実行しているだけ。
ループの中ではそれほどコストのかかることをやっていないつもりなのだが。

このあたりで力尽きたので，プロファイラーで詳しく調べる，といったことはやらなかった。

# 感想

素数日判定なんて，と思ったが，やってみたら意外と面白かった。

# 追記 2020-02-18：faster_prime で爆速化

[faster_prime](https://rubygems.org/gems/faster_prime) という gem の存在を知った。
これを使うと速くなるのではないか。
ちなみに pure Ruby のライブラリーである。

まず

```
gem install faster_prime
```

でインストールしておき，コード(5) に以下の変更を施して比較。

```diff
- require "prime"
+ require "faster_prime"
```

この記事を書いたときとは環境が違うので，変更前（prime 版）と変更後（faster_prime 版）を両方計測する。

変更前：18.9 秒
変更後：7.0 秒

マジかっ！