12109


pi(n)をn以下の素数の個数と定義します。  
pi(n)の正確な値を高速に求めるにはどのようなアルゴリズムを用いればよいでしょうか。

2×√n 程度のデータを用いて計算するコードを記しておきます。

(Ruby 2.2)

[code]

    def pi(n)
      m = Math.sqrt(n).to_i
      keys = (1..m).map{|i| n / i}
      keys += (1..keys[-1] - 1).to_a.reverse
      h = {}
      # 1を除いた個数
      keys.each{|i| h[i] = i - 1}
      # 「素数」もしくは「i以下の素数では割り切れない合成数」の個数
      (2..m).each{|i|
        if h[i] > h[i - 1] # このときiは素数
          hp = h[i - 1]
          i2 = i * i
          keys.each{|j|
            break if j < i2
            h[j] -= h[j / i] - hp # iで初めて割り切れる合成数(ちなみにi2以上）を除く
          }
        end
      }
      h[n]
    end
    
    p pi(10 ** 8)
    
[/code]

(Python 2.7)

[code]

    import math
    
    def pi(n):
        m = int(math.sqrt(n))
        keys = [n // i for i in range(1, m + 1)]
        keys += range(keys[-1] - 1, 0, -1)
        h = {i : i - 1 for i in keys}
        for i in range(2, m + 1):
            if h[i] > h[i - 1]:
                hp = h[i - 1]
                i2 = i * i
                for j in keys:
                    if j < i2: break
                    h[j] -= h[j // i] - hp
    
        return h[n]
    
    print pi(10 ** 8)
    
[/code]

（参考）  
(Ruby 2.2)  
上記コードと、Rosetta Code で紹介されている Wheel factorization を用いた方法  
（<http://rosettacode.org/wiki/Sieve_of_Eratosthenes#With_a_wheel>）とで比較してみた。

[code]

    require 'benchmark'
    puts Benchmark::CAPTION
    puts Benchmark.measure{
      p pi(10 ** 8)
    }
    puts Benchmark.measure{
      p eratosthenes2(10 ** 8).size
    }
    
[/code]

おおよそ以下のような結果になります。

[code]

          user     system      total        real
    5761455
      0.218000   0.000000   0.218000 (  0.210384)
    5761455
      8.362000   0.312000   8.674000 (  8.699177)
    
[/code]

速さで比較すると、  
上記コード＞ Wheel factorization を用いた方法（＞エラトステネスの篩）  
となっている。


CC BY-SA 3.0
2015-07-04T16:25:28.910
12033
2015-07-08T14:05:49.237
2015-07-07T13:51:52.100
5363
5363
question
ruby python アルゴリズム
n以下の素数の個数を高速に求めるには？


時間が書いてないので早いのか遅いのか分かりませんが・・・  
Cでエラトステネスの篩を使った場合１．６秒程度でカウントできます。  
~ ./a.out 100000000  
Time:1.625000

*** prime count: 5761455 ***  
~  
[エラトステネスの篩](http://www.crimsoneyes.jp/src/soe.c)  
ちなみに計算量はO(n log log n)です。


CC BY-SA 3.0
2015-07-05T09:20:35.217
12045
2015-07-05T13:14:54.393
2015-07-05T13:14:54.393
9528
9528
12033
answer


MathWorld の以下のページに method, time complexity, storage complexity
が載っているので、実装したときの速さの目安になるかと思われます。  
（<http://mathworld.wolfram.com/PrimeCountingFunction.html>）

余談ですが、上記ページには Π(10^25) (Büthe (2014)) までしか載っていませんが、OEISのA006880には Π(10^26)
まで載っています。  
（<https://oeis.org/A006880/list>）

なお、次のように改良してみました。

[code]

    def pi(n)
      m = Math.sqrt(n).to_i
      keys = (1..m).map{|i| n / i}
      keys += (1..keys[-1] - 1).to_a.reverse
      h = {}
      keys.each{|i| h[i] = i - 1}
      (2..m).each{|i|
        if h[i] > h[i - 1] # このときiは素数
          hp = h[i - 1]
          i2 = i * i
          h[n] -= h[n / i] - hp
          # h[n / j]は「i = jのときのh[n]に呼び出されるまでは計算しておく」ことにし、
          # i = jのとき、h[n / (j + 1)]から計算していくことにする
          keys[i..-1].each{|j|
            break if j < i2
            h[j] -= h[j / i] - hp
          }
        end
      }
      h[n]
    end
    
[/code]

n = 10^12 を計算させた結果

[code]

          user     system      total        real
    37607912018
    264.796000   0.078000 264.874000 (268.561594)
    
[/code]

元々のコードで計算させた結果

[code]

          user     system      total        real
    37607912018
    304.062000   0.078000 304.140000 (308.404353)
    
[/code]


CC BY-SA 3.0
2015-07-07T06:37:01.727
12109
2015-07-08T14:05:49.237
2015-07-08T14:05:49.237
5363
5363
12033
answer