Robert Kern wrote:
Rüdiger Werner wrote:
Hello!
Hi!
numpy questions are best answered on the numpy mailing list.
http://www.scipy.org/Mailing_Lists
Out of curiosity and to learn a little bit about the numpy package
i've tryed to implement
a vectorised version of the 'Sieve of Zakiya'.
While the code itself works fine it is astounding for me that the
numpy Version is almost 7 times slower than
the pure python version. I tryed to find out if i am doing something
wrong but wasn't able to find any answer.
The result of indexing into a numpy array is a numpy scalar object. We
do this instead of returning a float or an int because numpy supports
many more data types than just a C double or long, respectively. If I
have a uint16 array, indexing into it gives me a uint16 numpy scalar.
These are a little more complicated to set up than a regular Python
float or int, so they take more time to create.
Taking a closer look, that's only part of the problem. The larger problem is
that the numpy code wasn't vectorized enough to actually make use of numpy's
capabilities. You were paying for the overhead numpy imposes without coding in
such a way to make use of numpy's efficiencies. Here is my version that is more
than two times faster than the others. I didn't know exactly which details were
critical to the algorithm and which weren't (e.g. lndx+15), so there's some
gorpiness that can probably be removed. If you don't actually need to convert
back to a list, then you can be about 3 times faster than the others.
from numpy import arange, array, newaxis, zeros
def VSoZP5(val):
""" Vectorised Version of Sieve of Zakia"""
if val < 3 : return [2]
if val < 5 : return [2,3]
lndx = (val-1)/2
sieve = zeros(lndx+15,dtype="int32")
sieve[0:14] = 0 , 1 , 2, 0, 4, 5, 0, 7, 8, 0, 10, 0, 0, 13
_i = arange(14, lndx, 15)
if lndx > 14:
sieve[14:-14:15] = _i - 1
sieve[17:-11:15] = _i + 2
sieve[19:-9:15] = _i + 4
sieve[20:-8:15] = _i + 5
sieve[22:-6:15] = _i + 7
sieve[23:-5:15] = _i + 8
sieve[25:-3:15] = _i + 10
sieve[28::15] = _i + 13
om = array([9, 15, 18, 24, 27, 33, 42, 45], dtype="int32")
bm = array([7, 11, 13, 17, 19, 23, 29, 31], dtype="int32")
for i in xrange(2, (int(sqrt(val))-3/2)+1):
if not sieve[i]: continue
k = 30*i+45
rm = (bm * i) + om
rm2 = (rm + k * arange(0, lndx/k + 1)[:,newaxis]).flat
rm3 = rm2[rm2 < len(sieve)]
sieve[rm3] = 0
sieve = sieve[2:lndx]
primes = (2*arange(2, lndx) + 3)[sieve > 0].tolist()
primes[0:0] = [2,3,5]
return primes
--
Robert Kern
"I have come to believe that the whole world is an enigma, a harmless enigma
that is made terrible by our own mad attempt to interpret it as though it had
an underlying truth."
-- Umberto Eco
--
http://mail.python.org/mailman/listinfo/python-list
