Skip to content

Commit 9753572

Browse files
coellscoells
committed
day 100
1 parent 7b8880b commit 9753572

File tree

1 file changed

+216
-0
lines changed

1 file changed

+216
-0
lines changed

‎day I00 - segmented eratosthenes sieve.ipynb‎

Lines changed: 216 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,216 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": 1,
6+
"metadata": {
7+
"collapsed": true
8+
},
9+
"outputs": [],
10+
"source": [
11+
"%load_ext Cython"
12+
]
13+
},
14+
{
15+
"cell_type": "markdown",
16+
"metadata": {},
17+
"source": [
18+
"## algorithm"
19+
]
20+
},
21+
{
22+
"cell_type": "code",
23+
"execution_count": 2,
24+
"metadata": {
25+
"collapsed": true
26+
},
27+
"outputs": [],
28+
"source": [
29+
"%%cython\n",
30+
"\n",
31+
"from libc.stdlib cimport malloc, free\n",
32+
"\n",
33+
"DEF LIMIT = 1024 * 31\n",
34+
"DEF PRIME = 1024 * 4\n",
35+
"DEF SIEVE = 1024 * 32\n",
36+
"\n",
37+
"cdef inline int imin(int a, int b) nogil:\n",
38+
" return a if a < b else b\n",
39+
"\n",
40+
"cdef inline int memset(char *p, int n) nogil:\n",
41+
" cdef:\n",
42+
" short *q = <short *>p\n",
43+
" int i, j = 0\n",
44+
"\n",
45+
" for i in range((n + 1) >> 1):\n",
46+
" j += q[i]\n",
47+
" q[i] = 0x0100\n",
48+
"\n",
49+
" return j >> 8\n",
50+
"\n",
51+
"cdef int naive_sieve(char *sieve, int *primes, int *offsets, int n) nogil:\n",
52+
" cdef int i, j\n",
53+
"\n",
54+
" memset(sieve, n)\n",
55+
"\n",
56+
" for i in range(3, n, 2):\n",
57+
" if sieve[i]:\n",
58+
" j = i * i\n",
59+
" while j < n:\n",
60+
" sieve[j] = 0\n",
61+
" j += i << 1\n",
62+
"\n",
63+
" primes[0] = i\n",
64+
" offsets[0] = j\n",
65+
" primes += 1\n",
66+
" offsets += 1\n",
67+
"\n",
68+
" primes[0] = 0\n",
69+
" offsets[0] = 0\n",
70+
"\n",
71+
" return memset(sieve, n)\n",
72+
"\n",
73+
"cdef int segmented_sieve(char *sieve, int *primes, int *offsets, int k, int n) nogil:\n",
74+
" cdef int i\n",
75+
"\n",
76+
" while primes[0]:\n",
77+
" i = offsets[0] - k\n",
78+
" while i < n:\n",
79+
" sieve[i] = 0\n",
80+
" i += primes[0] << 1\n",
81+
" offsets[0] = i + k\n",
82+
"\n",
83+
" primes += 1\n",
84+
" offsets += 1\n",
85+
"\n",
86+
" return memset(sieve, n)\n",
87+
"\n",
88+
"cpdef int eratosthenes(int n) nogil:\n",
89+
" cdef:\n",
90+
" char *sieve\n",
91+
" int *primes\n",
92+
" int *offsets\n",
93+
" int k, total\n",
94+
"\n",
95+
" if n > LIMIT * LIMIT:\n",
96+
" return -1\n",
97+
"\n",
98+
" sieve = <char *>malloc(SIEVE)\n",
99+
" primes = <int *>malloc(PRIME * sizeof(int))\n",
100+
" offsets = <int *>malloc(PRIME * sizeof(int))\n",
101+
"\n",
102+
" total = naive_sieve(sieve, primes, offsets, imin(n, LIMIT))\n",
103+
"\n",
104+
" memset(sieve, SIEVE)\n",
105+
" k = LIMIT\n",
106+
" n -= LIMIT\n",
107+
"\n",
108+
" while n > 0:\n",
109+
" total += segmented_sieve(sieve, primes, offsets, k, imin(n, SIEVE))\n",
110+
" k += SIEVE\n",
111+
" n -= SIEVE\n",
112+
"\n",
113+
" free(sieve)\n",
114+
" free(primes)\n",
115+
" free(offsets)\n",
116+
"\n",
117+
" return total"
118+
]
119+
},
120+
{
121+
"cell_type": "markdown",
122+
"metadata": {},
123+
"source": [
124+
"## run"
125+
]
126+
},
127+
{
128+
"cell_type": "code",
129+
"execution_count": 3,
130+
"metadata": {},
131+
"outputs": [
132+
{
133+
"name": "stdout",
134+
"output_type": "stream",
135+
"text": [
136+
"primes below 10**1: 4\n",
137+
"primes below 10**2: 25\n",
138+
"primes below 10**3: 168\n",
139+
"primes below 10**4: 1229\n",
140+
"primes below 10**5: 9592\n",
141+
"primes below 10**6: 78498\n",
142+
"primes below 10**7: 664579\n",
143+
"primes below 10**8: 5761455\n",
144+
"primes below 10**9: 50847534\n"
145+
]
146+
}
147+
],
148+
"source": [
149+
"for i in range(1, 10):\n",
150+
" print('primes below 10**%d: %d' % (i, eratosthenes(10**i)))"
151+
]
152+
},
153+
{
154+
"cell_type": "markdown",
155+
"metadata": {},
156+
"source": [
157+
"## timeit"
158+
]
159+
},
160+
{
161+
"cell_type": "code",
162+
"execution_count": 4,
163+
"metadata": {},
164+
"outputs": [
165+
{
166+
"name": "stdout",
167+
"output_type": "stream",
168+
"text": [
169+
"10000 loops, best of 3: 56.6 µs per loop\n",
170+
"1000 loops, best of 3: 574 µs per loop\n",
171+
"100 loops, best of 3: 6.07 ms per loop\n",
172+
"10 loops, best of 3: 68.4 ms per loop\n",
173+
"1 loop, best of 3: 863 ms per loop\n"
174+
]
175+
}
176+
],
177+
"source": [
178+
"%timeit eratosthenes(1024 * 31)\n",
179+
"%timeit eratosthenes(10**6)\n",
180+
"%timeit eratosthenes(10**7)\n",
181+
"%timeit eratosthenes(10**8)\n",
182+
"%timeit eratosthenes(10**9)"
183+
]
184+
},
185+
{
186+
"cell_type": "code",
187+
"execution_count": null,
188+
"metadata": {
189+
"collapsed": true
190+
},
191+
"outputs": [],
192+
"source": []
193+
}
194+
],
195+
"metadata": {
196+
"kernelspec": {
197+
"display_name": "Python 3",
198+
"language": "python",
199+
"name": "python3"
200+
},
201+
"language_info": {
202+
"codemirror_mode": {
203+
"name": "ipython",
204+
"version": 3
205+
},
206+
"file_extension": ".py",
207+
"mimetype": "text/x-python",
208+
"name": "python",
209+
"nbconvert_exporter": "python",
210+
"pygments_lexer": "ipython3",
211+
"version": "3.6.1"
212+
}
213+
},
214+
"nbformat": 4,
215+
"nbformat_minor": 2
216+
}

0 commit comments

Comments
 (0)