# [PYTHON MATRIX-SIG] Another problem?

Jim Hugunin hugunin@mit.edu
Fri, 16 Aug 1996 11:14:41 -0400

```Rob.Hooft@embl-heidelberg.de wrote:
>
> While hacking some more on my "code-cracker" :-) script, I found what I think
> is another problem:
>
>    nuNumPy% python
>    Python 1.4b2 (Aug 12 1996) [C]
>    Copyright 1991-1996 Stichting Mathematisch Centrum, Amsterdam
>    >>> from Numeric import *
>    >>> nonzero([0,2,0,4])
>    1 1 3 3 3 3
>    >>>
>
> from the description I read nonzero should return [1,3] here. It looks
> like a mixture with repeat...

nonzero is implemented using repeat, I made a mistake in writing it to
give that strange mixture.  This is fixed in alpha2 which will be
released before I go home tonight.

> Furthermore, I did some more work on the script, got it working, but
> the results are very disappointing:
>
>   Can divide by 2
>   Can divide by 29
>   factor took  10.577 seconds
>   Can divide by 2
>   Can divide by 29
>   afactor took   7.916 seconds

I'm running on a faster machine than you, so my numbers are faster in
both cases, still the relative comparision is valid.

Standard python version
factor took   2.708 seconds
Using NumPy  - properly ;)
afactor took   0.322 seconds

Here's my new version of your afactor and asieve functions.  The changes
to afactor are actually not needed to get the speed improvement.  I'd be
tempted to claim that the new version of asieve is easier to understand
as well as faster, the changes to afactor are a little awkward but I
thought you might be interested in seeing the style anyway.

I doubt that my new version is more than a factor of 2 slower than a C
implementation of the same algorithm (though feel free to show me wrong
with a code sample).

def asieve(max,atype='l'):
size=int(sqrt(max))
if size<5:
known_primes=[2,3]
else:
known_primes=asieve(size)

#Initially assume all numbers are prime
numbers=ones([max+1], atype)
#0 and 1 are not prime
numbers[0:2]=0
#print trials
for i in known_primes:
#Set multiples of i to be nonprime
numbers[(i*i)::i] = 0
#Those entries which are nonzero are prime
return nonzero(numbers)

def afactor(num,upto=1000000):
#Calculate all of the primes up to upto
s = asieve(upto)
#Return just those primes that divide evenly into num
for d in take(s, nonzero(equal(num % s, 0))):
print "Can divide by %d"%d

-Jim

=================
MATRIX-SIG  - SIG on Matrix Math for Python

send messages to: matrix-sig@python.org