回答編集履歴
1
追記
answer
CHANGED
|
@@ -140,4 +140,88 @@
|
|
|
140
140
|
|
|
141
141
|
この状態で10^8までだと、手元環境ではエラトステネスの篩の倍弱の実行時間になります。
|
|
142
142
|
|
|
143
|
-
軽くプロファイリングした感じだと、一番内側のforループで大半の時間を食っているようです。
|
|
143
|
+
軽くプロファイリングした感じだと、一番内側のforループで大半の時間を食っているようです。
|
|
144
|
+
|
|
145
|
+
|
|
146
|
+
### 追記
|
|
147
|
+
最終的に、numbaのちからで区間篩の方が速くなるところまで持っていきました。numbaのjitコンパイルを通すためにコードをまただいぶ書き換えてあります。
|
|
148
|
+
```python
|
|
149
|
+
import timeit
|
|
150
|
+
import math
|
|
151
|
+
|
|
152
|
+
import numba
|
|
153
|
+
|
|
154
|
+
#旧：エラトステネスの篩
|
|
155
|
+
@numba.jit(nopython=True)
|
|
156
|
+
def eratosthenes_list(n):
|
|
157
|
+
primes = [True] * n
|
|
158
|
+
result = [2]
|
|
159
|
+
for prime in range(3, n, 2):
|
|
160
|
+
if primes[prime]:
|
|
161
|
+
result.append(prime)
|
|
162
|
+
for i in range(prime * prime, n, prime * 2):
|
|
163
|
+
primes[i] = False
|
|
164
|
+
return result
|
|
165
|
+
|
|
166
|
+
#新：区間篩
|
|
167
|
+
@numba.jit(nopython=True)
|
|
168
|
+
def primesList(n):
|
|
169
|
+
n -= 1
|
|
170
|
+
SIEVE_MAX = int(math.sqrt(n)) + 1
|
|
171
|
+
|
|
172
|
+
sieve = [True] * (SIEVE_MAX)
|
|
173
|
+
sieve[0] = sieve[1] = False
|
|
174
|
+
for i in range(2, SIEVE_MAX):
|
|
175
|
+
for j in range(i * 2, SIEVE_MAX, i):
|
|
176
|
+
sieve[j] = False
|
|
177
|
+
|
|
178
|
+
primes = []
|
|
179
|
+
for i, x in enumerate(sieve):
|
|
180
|
+
if x:
|
|
181
|
+
primes.append(i)
|
|
182
|
+
|
|
183
|
+
for i in range(1, (n // (SIEVE_MAX))+1):
|
|
184
|
+
sieve = [True] * len(sieve)
|
|
185
|
+
start_index = SIEVE_MAX * i
|
|
186
|
+
|
|
187
|
+
for p in primes:
|
|
188
|
+
if (p * p > start_index + SIEVE_MAX):
|
|
189
|
+
break
|
|
190
|
+
|
|
191
|
+
jj = (p - (start_index % p)) % p
|
|
192
|
+
for j in range(jj, SIEVE_MAX, p):
|
|
193
|
+
sieve[j] = False
|
|
194
|
+
|
|
195
|
+
tmp = n - start_index + 1
|
|
196
|
+
if tmp < SIEVE_MAX:
|
|
197
|
+
x = tmp
|
|
198
|
+
else:
|
|
199
|
+
x = SIEVE_MAX
|
|
200
|
+
for j in range(x):
|
|
201
|
+
if sieve[j]:
|
|
202
|
+
primes.append(start_index + j)
|
|
203
|
+
return primes
|
|
204
|
+
|
|
205
|
+
if True:
|
|
206
|
+
for nn in range(10, 10**3):
|
|
207
|
+
a = eratosthenes_list(nn)
|
|
208
|
+
b = primesList(nn)
|
|
209
|
+
if not a == b:
|
|
210
|
+
print(nn)
|
|
211
|
+
print("miss")
|
|
212
|
+
print(len(a), len(b))
|
|
213
|
+
print(a)
|
|
214
|
+
print(b)
|
|
215
|
+
break
|
|
216
|
+
else:
|
|
217
|
+
print("all ok")
|
|
218
|
+
|
|
219
|
+
nn = 10**8
|
|
220
|
+
print(timeit.timeit(lambda: (eratosthenes_list(nn)), number = 1))
|
|
221
|
+
print(timeit.timeit(lambda: (primesList(nn)), number = 1))
|
|
222
|
+
""" =>
|
|
223
|
+
all ok
|
|
224
|
+
1.2398602159810252
|
|
225
|
+
1.059988280001562
|
|
226
|
+
"""
|
|
227
|
+
```
|