14185


<https://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-
below-n/3035188#3035188> に書かれてあるprimesfrom2toという関数を使って、以下のコードを実行しました。

[code]

    import numpy
    
    def primesfrom2to(n):
        """ Input n>=6, Returns a array of primes, 2 <= p < n """
        sieve = numpy.ones(n/3 + (n%6==2), dtype=numpy.bool)
        for i in xrange(1,int(n**0.5)/3+1):
            if sieve[i]:
                k=3*i+1|1
                sieve[       k*k/3     ::2*k] = False
                sieve[k*(k-2*(i&1)+4)/3::2*k] = False
        return numpy.r_[2,3,((3*numpy.nonzero(sieve)[0][1:]+1)|1)]
    
    N = 10 ** 6
    print primesfrom2to(N)
    
    total = 0
    for i in primesfrom2to(N):
        total += i
    print total
    
[/code]

すると、実行結果は次のようになりました。

[code]

    [     2      3      5 ..., 999961 999979 999983]
    150808.py:16: RuntimeWarning: overflow encountered in int_scalars
      total += i
    -1104303641
    
[/code]

注）150808.py はファイル名です。

10^6 までの素数は求まっているのに、和が求まっていません。  
このことを解消するにはどうすればよろしいでしょうか？  
（ちなみに N = 10 ** 5 のとき素数の和は求まります。）

（追記）2015/08/08 22:25  
total = 0L  
となおすだけで上記問題は解決しました。  
なお、修正後 N = 10 ** 9 のとき素数の和は求まりますが、  
N = 10 ** 10 以上のとき素数の和は求まりませんでした。  
これは、Numpy では計算できないことなのでしょうか？


CC BY-SA 3.0
2015-08-08T11:51:46.730
13175
2015-08-09T13:24:21.000
2017-05-23T12:38:56.467
-1
5363
question
python
Numpy で素数の和を求めるには？


（追記）の  
「N = 10 ** 10 以上のとき素数の和は求まりませんでした。  
これは、Numpy では計算できないことなのでしょうか？」  
に対する回答です。

私の環境は  
OS：Windows 64bit版  
Numpy:32bit版※  
でした。  
※Pythonによるデータ分析入門1.4.1の訳注に「Windowsが64bit版であっても、64bit版のCanopy
Expressが正常に動作しない場合がありました。この場合を回避するために、翻訳時点（2013年8月現在）では、32bit版を利用することをおすすめします」とあったので、32bit版をインストールしました。

よって、Numpyの32bit版という理由で  
N = 10 ** 10 以上のとき素数の和は求まらない  
と言えそうです。

（おまけ）  
N が 10, 10^2, … ,10^9 のときの値を列挙してみた。

[code]

    import numpy
    
    def primesfrom2to(n):
        """ Input n>=6, Returns a array of primes, 2 <= p < n """
        sieve = numpy.ones(n/3 + (n%6==2), dtype=numpy.bool)
        for i in xrange(1,int(n**0.5)/3+1):
            if sieve[i]:
                k=3*i+1|1
                sieve[       k*k/3     ::2*k] = False
                sieve[k*(k-2*(i&1)+4)/3::2*k] = False
        return numpy.r_[2,3,((3*numpy.nonzero(sieve)[0][1:]+1)|1)]
    
    for i in range(1, 10):
        print primesfrom2to(10 ** i).sum(dtype=numpy.int64)
    
[/code]

実行結果  
17  
1060  
76127  
5736396  
454396537  
37550402023  
3203324994356  
279209790387276  
24739512092254535


CC BY-SA 3.0
2015-08-08T16:01:58.530
14185
2015-08-08T16:46:55.027
2015-08-08T16:46:55.027
5363
5363
13175
answer


有限桁の２進表現を使用するnp.int64等では計算できる数字の大きさに限度があります。  
任意多倍長で計算する機能はおそらくNumPy自体には含まれていません  
（あまり真面目に確認してませんが）。  
例えば代表的な多倍長のライブラリである[GMP](https://ja.wikipedia.org/wiki/GNU_Multi-
Precision_Library)のラッパーである[gmpy2](https://gmpy2.readthedocs.org/en/latest/)の[mpz](https://gmpy2.readthedocs.org/en/latest/mpz.html#examples)を使えば

[code]

    import gmpy2
    total = gmpy2.mpz(0)
    for i in primesfrom2to(N):
        total += 1
    
[/code]

として正しい値が計算できると思います。


CC BY-SA 3.0
2015-08-09T13:24:21.000
14212
2015-08-09T13:24:21.000
9572
13175
answer