Is there a way to list super primes between 1 and n in NumPy-CodePudding

Is there a way to list super primes (primes in a prime position wikipedia article) between 1 and n using NumPy.

CodePudding user response：

I got this without Numpy, is that okay? here is the code based on

Sieve of Atkin

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)

CodePudding user response：

Numpy Code Version

Sieve modification of Numpy Sieve by user2357112 supports 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

Test

print(superPrimes(200))
# Output: [  3   5  11  17  31  41  59  67  83 109 127 157 179 191]