【Paiza】大きい整数をエラトステネスの篩にかける

これに挑みました。

そもそもエラトステネスの篩（ふるい）って何？

指定された整数以下の全ての素数を発見するための単純なアルゴリズムである。古代ギリシアの科学者、エラトステネスが考案したとされるため、この名がついている。

上記ページには、アルゴリズムの内容が分かりやすいGIFアニメもあります。

Github Copilotさんに聞いてみた

プログラミングの勉強にあたって、基礎知識に幅広く触れ、コードリーディングをつうじてアルゴリズムをよりたくさん身につけたい、と思っているので、根性に頼らずAI様に解かせてみました。 副操縦士にメインパイロットさせとるやん、というツッコミは無しで。

konna_kanji.rb

n = gets.to_i
is_prime = Array.new(n+1, true)
is_prime[0] = false
is_prime[1] = false

(2..n).each do |i|
  if is_prime[i]
    (2..n).each do |j|
      is_prime[i*j] = false if i*j <= n
    end
  end
end

puts is_prime[n] ? "YES" : "NO"

予め配列の0、１にfalseを入れる という仕様はpaizaの出題に合わせて実装しています。

2以上の整数iの倍数を、i == nに至るまで配列へ入れていくというコードですね。
すでに配列[i]にfalseを入れてある場合はスキップする という処理は私はすぐに思いつかなかったので、さすがAIといった感じです。

試しに n = 1223 で計算した時、 i * j の部分をコンソールに出力してみたところ、

出力結果

~~~
1085
1116
1147
1178
1209
1240
1271
1302
~~~

1223がないため、nは素数。ちゃんと機能していそうですね！

…が、しかしnが大きい場合

paizaさんのところでは、テストケースとしてnに大きい整数を用意していました。ここでは仮に n=400000 とします。（どこまで外部に情報を書いていいのか分からないので、すこしぼやかしています）

上記コードだと、この時に処理が重すぎて、制限時間（paizaさんのところでRuntime Error判定に使用している一定の時間）までに間に合いませんでした。

ルートを使おう

今回のアルゴリズムに限らず、素数の判定には平方根を使うのが良いそうです。

仮に、ある自然数Nが√Nより大きい素数Pを素因数に持っているとします。
その時、N=P×M（Mは自然数）と表すことができ、
このとき P>√N より M<√N
つまり P>M となります。
（中略）
つまり、Nが√Nより大きい素数を素因数に持つならば、
√N以下の素数で調べたときにすでに消去されているはずなのです。
（ベストアンサーより引用）

数学に強いわけではないのですが、 √Nより大きい素数Pがある ときに　Nは、Pの何かしらの自然数Mの乗算で表すことができる ため、今回のアルゴリズムでは その自然数Mの倍数を先に見つけることができる 、ということでしょうか。

tatoeba.rb

def prime_factorization(number)
  prime_factors = []
  i = 2
  while i * i <= number
    while number % i == 0
      prime_factors << i
      number /= i
    end
    i += 1
  end
  prime_factors << number if number > 1
  prime_factors
end

prime_factorization(2446)
# => [2, 1223] 
# 1223は素数。今回だと、この「2」の倍数が、先に篩へ含まれることになる。
# 参考に↓
Math.sqrt(2446)
# => 49.457052075512955
# 2の12乗は4096なので、確かに49に至るまでに2446は先に網羅されているはずであろう。

計算量を抑えよう

修正したのがこちら。

refactored_program.rb

n = gets.to_i
is_prime = Array.new(n+1, true)
is_prime[0] = false
is_prime[1] = false
# count = 0 <= 計算量の調査

(2..Math.sqrt(n).ceil).each do |i| # <= 変更点
  if is_prime[i]
    (2..n).each do |j|
      # count += 1 <= 計算量の調査
      is_prime[i*j] = false if i*j <= n
    end
  end
end

puts is_prime[n] ? "YES" : "NO"

これでpaizaさんが用意した n = 400000 のケースも通りました。（Runtime Error→約1.5秒で完了という改善）

簡単に計算量を調査する

n = 2446 の時、二重ループの内側に入った際に回数 count は36,675回 でした。

ループを行った総数は Math.sqrt(n).ceil * n = 122,300回なので、内側ループ手前の if is_prime[i] によるガードによって85,625回分の処理をスキップできています。ですが、それでもまあまあ膨れ上がりますね。（参考に 2446**2 は5,982,916）

すでにチューニング案を出したということで、これ以上の詮索はやめておきます。もう少しプログラミングが上達した際は、このアルゴリズムに関してさらなる速度向上を目指したいです。