[Bio] / Numeric / Demo / sieve.py Repository:
ViewVC logotype

View of /Numeric/Demo/sieve.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (download) (as text) (annotate) (vendor branch)
Thu Mar 10 23:06:43 2005 UTC (14 years, 9 months ago) by efrank
Branch: MAIN, Numeric-23-7
CVS Tags: lwc, init, HEAD
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())
	#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))


MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3