#! /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())
#add(numbers, array(1, 'b'), numbers)
#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 = add.accumulate(ones([len(numbers)]))-1
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))
syntax highlighted by Code2HTML, v. 0.9.1