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

Annotation of /Numeric/Demo/sieve.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (view) (download) (as text)

1 : efrank 1.1 #! /usr/bin/env python
2 :     from Numeric import *
3 :    
4 :     def sieve(max):
5 :     numbers=range(max+1)
6 :     size=int(math.sqrt(max))
7 :     if size<5:
8 :     trials=[2,3]
9 :     else:
10 :     trials=sieve(size)
11 :     for i in trials:
12 :     try:
13 :     j=i*i
14 :     while 1:
15 :     numbers[j]=0
16 :     j=j+i
17 :     except IndexError:
18 :     pass
19 :     return filter(lambda x:x>1,numbers)
20 :    
21 :     def factor(num,upto=1000000):
22 :     for p in sieve(upto):
23 :     if num%p==0:
24 :     print "Can divide by %d"%p
25 :    
26 :     def asieve(max,atype='l'):
27 :     times = []
28 :     #times.append(time.time())
29 :     size=int(sqrt(max))
30 :     if size<5:
31 :     known_primes=[2,3]
32 :     else:
33 :     known_primes=asieve(size)
34 :    
35 :     #Initially assume all numbers are prime
36 :     numbers=ones([max+1], 'b') #atype)
37 :     #times.append(time.time())
38 :     #add(numbers, array(1, 'b'), numbers)
39 :     #times.append(time.time())
40 :     #0 and 1 are not prime
41 :     numbers[0:2]=0
42 :     #print trials
43 :     for i in known_primes:
44 :     #Set multiples of i to be nonprime
45 :     numbers[(i*i)::i] = 0
46 :     #times.append(time.time())
47 :     #Those entries which are nonzero are prime
48 :     #a = add.accumulate(ones([len(numbers)]))-1
49 :     a = arange(len(numbers))
50 :     #times.append(time.time())
51 :     r = repeat(a, numbers)#numbers)nonzero(numbers)
52 :     #times.append(time.time())
53 :     #times = array(times)
54 :     #print times[1:]-times[:-1]
55 :     return r
56 :    
57 :     def afactor(num,upto=1000000):
58 :     #Calculate all of the primes up to upto
59 :     s = asieve(upto)
60 :     #Return just those primes that divide evenly into num
61 :     for d in take(s, nonzero(equal(num % s, 0))):
62 :     print "Can divide by %d"%d
63 :    
64 :    
65 :     import time
66 :    
67 :     #start = time.time()
68 :     #factor(129753308,upto=99999)
69 :     #stop = time.time()
70 :     #print "factor took %7.3f seconds"%((stop-start))
71 :     start = time.time()
72 :     asieve(999999)
73 :     afactor(129753308,upto=999999)
74 :     stop = time.time()
75 :     print "afactor took %7.3f seconds"%((stop-start))
76 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3