[Bio] / Numeric / Demo / sieve.py Repository: Repository Listing Bio

# View of /Numeric/Demo/sieve.py

Thu Mar 10 23:06:43 2005 UTC (14 years, 10 months ago) by efrank
Branch: MAIN, Numeric-23-7
Changes since 1.1: +0 -0 lines
```Numeric 23.7 for env migration
```

```#! /usr/bin/env python
from Numeric import *

def sieve(max):
numbers=range(max+1)
size=int(math.sqrt(max))
if size<5:
trials=[2,3]
else:
trials=sieve(size)
for i in trials:
try:
j=i*i
while 1:
numbers[j]=0
j=j+i
except IndexError:
pass
return filter(lambda x:x>1,numbers)

def factor(num,upto=1000000):
for p in sieve(upto):
if num%p==0:
print "Can divide by %d"%p

def asieve(max,atype='l'):
times = []
#times.append(time.time())
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], 'b') #atype)
#times.append(time.time())
#times.append(time.time())
#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
#times.append(time.time())
#Those entries which are nonzero are prime
a = arange(len(numbers))
#times.append(time.time())
r = repeat(a, numbers)#numbers)nonzero(numbers)
#times.append(time.time())
#times = array(times)
#print times[1:]-times[:-1]
return r

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

import time

#start = time.time()
#factor(129753308,upto=99999)
#stop = time.time()
#print "factor took %7.3f seconds"%((stop-start))
start = time.time()
asieve(999999)
afactor(129753308,upto=999999)
stop = time.time()
print "afactor took %7.3f seconds"%((stop-start))

```