有沒有辦法使用 NumPy列出1 和 n 之間的超級素數（維基百科文章中素數位置的素數）。
uj5u.com熱心網友回復：
我在沒有Numpy 的情況下得到了這個，可以嗎？這是基于的代碼
阿特金篩
import math
is_prime = list()
limit = 100
for i in range(5, limit):
is_prime.append(False)
for x in range(1, int(math.sqrt(limit)) 1):
for y in range(1, int(math.sqrt(limit)) 1):
n = 4 * x ** 2 y ** 2
if n <= limit and (n % 12 == 1 or n % 12 == 5) and n <= len(is_prime):
# print "1st if"
is_prime[n] = not is_prime[n]
n = 3 * x ** 2 y ** 2
if n <= limit and n % 12 == 7:
# print "Second if"
is_prime[n] = not is_prime[n]
n = 3 * x ** 2 - y ** 2
if x > y and n <= limit and n % 12 == 11:
# print "third if"
is_prime[n] = not is_prime[n]
for n in range(5, int(math.sqrt(limit))):
if is_prime[n]:
for k in range(n ** 2, limit 1, n ** 2):
if k <= len(is_prime):
is_prime[k] = False
for n in range(5, limit):
if n < len(is_prime) and is_prime[n]:
print(n)
uj5u.com熱心網友回復：
Numpy 代碼版本
用戶2357112對Numpy Sieve的Sieve修改支持Monica
import numpy as np
def sieve(n):
'''
Sieve of Erastosenes using numpy
'''
flags = np.ones(n, dtype=bool) # Set all values to True
# Set 0, 1, and even numbers > 2 to False
flags[0] = flags[1] = False
flags[4:n:2] = False
for i in range(3, n, 2):
# Check for odd primes
if flags[i]:
flags[i*i::i] = False
return np.flatnonzero(flags) # Indexes of True are prime
def superPrimes(n):
'''
Find superprimes
'''
# Find primes up to n using sieve
p = sieve(n)
# p(i) denotes the ith prime number, the numbers in this sequence are those of the form p(p(i))
# Want to find numbers we can index in the list (i.e. has to be less than len(p))
indexes = p[p < len(p)] # base 1 indexes within range
return p[indexes-1] # Subtract 1 since primes array is base 0 indexed
測驗
print(superPrimes(200))
# Output: [ 3 5 11 17 31 41 59 67 83 109 127 157 179 191]
轉載請註明出處，本文鏈接：https://www.uj5u.com/shujuku/398141.html
上一篇：逐行條件組合陣列