从numpy数组中获取质数

4

I have a numpy array like,

nums = np.array([17, 18, 19, 20, 21, 22, 23])

如何以Pythonic的方式从数组中过滤出质数?我知道可以执行简单的筛选操作,例如:

nums[nums > 20]   #array([21, 22, 23])

有没有一种方法可以传递lambda函数来进行过滤?
期望输出:array([17, 19, 23])

2
欢迎来到StackOverflow。问题应该至少尝试解决问题,或者您可以告诉我们您尝试了什么。在这种情况下,一个好的开始是尝试检查它们是否可被任何整数整除?提示:模运算符可能会有所帮助。 - MSeifert
1
@MSeifert 的讽刺意味十分明显。楼主,gmpy 已经内置了 Miller-Rabin 素性测试。 - Goodies
找到一些代码(或编写代码!)以获取质数列表,该列表达到您的nums数组的最高限制。然后使用“for n in primes”进行迭代。 - joel goldstick
这个可能会对你有趣:https://dev59.com/snI95IYBdhLWcg3w-DH0 - MaxU - stand with Ukraine
2
@kmario23,通常当你想找到“更好的方法”时,你通常已经有了一个可行的方法,如果你将这个有效的尝试添加到你的问题中,那会很有帮助。 - Padraic Cunningham
显示剩余3条评论
5个回答

10

我会使用gmpy或者一个开发了优秀素数测试算法的第三方库来完成。米勒-拉宾素数测试通常是非常安全(而且快速!)的选择。如果你只想用慢速的方式,可以这样做:

import numpy as np
import math

def is_prime(n):
    if n % 2 == 0 and n > 2: 
        return False
    return all(n % i for i in range(3, int(math.sqrt(n)) + 1, 2))

a = np.arange(1, 10**3)
foo = np.vectorize(is_prime)
pbools = foo(a)
primes = np.extract(pbools, a)
primes  # => Output below
array([  1,   2,   3,   5,   7,  11,  13,  17,  19,  23,  29,  31,  37,
        41,  43,  47,  53,  59,  61,  67,  71,  73,  79,  83,  89,  97,
       101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163,
       167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233,
       239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311,
       313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389,
       397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
       467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563,
       569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641,
       643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727,
       733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821,
       823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907,
       911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997])
如果您想要筛选掉质数,只需在pbools变量上调用np.invert。对于任何谓词,都可以采用相同的方法。您还可以将lambda传递给vectorize。例如,假设我们只想筛选出同时满足是质数且与5的除法余数为1的数字(不管出于什么原因)。
import numpy as np
import math

def is_prime(n):
    if n % 2 == 0 and n > 2: 
        return False
    return all(n % i for i in range(3, int(math.sqrt(n)) + 1, 2))

a = np.arange(1, 10**3)
foo = np.vectorize(lambda x: (not (x + 1) % 5 or not (x - 1) % 5) and is_prime(x))
primes = a[foo(a)]  # => Shorthand.... Output below
array([  1,  11,  19,  29,  31,  41,  59,  61,  71,  79,  89, 101, 109,
   131, 139, 149, 151, 179, 181, 191, 199, 211, 229, 239, 241, 251,
   269, 271, 281, 311, 331, 349, 359, 379, 389, 401, 409, 419, 421,
   431, 439, 449, 461, 479, 491, 499, 509, 521, 541, 569, 571, 599,
   601, 619, 631, 641, 659, 661, 691, 701, 709, 719, 739, 751, 761,
   769, 809, 811, 821, 829, 839, 859, 881, 911, 919, 929, 941, 971, 991])

2
如果你关心速度和效率,我建议你使用其中一个最快的质数筛法numpy.intersect1d()函数:
import numpy as np

def primesfrom2to(n):
    # https://dev59.com/snI95IYBdhLWcg3w-DH0#3035188
    """ Input n>=6, Returns a array of primes, 2 <= p < n """
    sieve = np.ones(n//3 + (n%6==2), dtype=np.bool)
    sieve[0] = False
    for i in range(int(n**0.5)//3+1):
        if sieve[i]:
            k=3*i+1|1
            sieve[      ((k*k)//3)      ::2*k] = False
            sieve[(k*k+4*k-2*k*(i&1))//3::2*k] = False
    return np.r_[2,3,((3*np.nonzero(sieve)[0]+1)|1)]

# generate 100.000 random integers from 1 to 1.000.000.000
a1 = np.random.randint(1, 10**9, 100000)
# generate all primes that are equal or less than a1.max() 
primes = primesfrom2to(a1.max())

# print result    
print(np.intersect1d(primes, a1))

1

看起来你的问题并不是关于质数,而是关于如何将函数应用于 numpy 数组。我使用了简单的 is_odd 示例。也许你正在寻找 np.vectorize

In [34]: nums = np.array([17, 18, 19, 20, 21, 22, 23])

In [35]: def is_odd(n):
    if n % 2 == 1:
        return True
    return False
   ....: 

In [36]: is_odd_v = np.vectorize(is_odd)

In [37]: nums[is_odd_v(nums)]
Out[37]: array([17, 19, 21, 23]

如果我没记错的话,np.vectorize 主要用于方便起见,并且性能不是很好。


但这并不总是会返回质数吧?它会返回39,但它并不是一个质数。 - user6084800
3
这不会返回质数,它会返回奇数。但是如果你自己编写了is_prime函数,它将返回一个经过筛选的质数np.array。当然,我的解决方案很朴素且未经优化。 - Akavall

1

有这样的设置:

import numpy as np
import math
nums = np.array([17, 18, 19, 20, 21, 22, 23])

因此现在我们创建一个包含所有可能整数候选项的数组:

divisors = np.arange(2,int(math.sqrt(np.max(nums)))+1) # Numbers from 2 to sqrt(max(nums))
print(divisors)
# [2 3 4]

现在对数组应用模运算,但使用不同的维度,以便我们检查每个数与每个除数的情况:
print(nums[:,None] % divisors[None,:]) # Modulo operation on each element (0 means divisible)
[[1 2 1]
 [0 0 2]
 [1 1 3]
 [0 2 0]
 [1 0 1]
 [0 1 2]
 [1 2 3]]
现在我们如何得到质数...我们检查是否有零结果的行:
print(np.min(nums[:,None] % divisors[None,:], axis=1)) # Minimum of the modulo for that element
# [1 0 1 0 0 0 1]

然后创建一个掩码来索引它们:

print(nums[np.min(nums[:,None] % divisors[None,:], axis=1) > 0]) # So index them
# [17 19 23]

所以最终你所需的就是:

nums = np.array([17, 18, 19, 20, 21, 22, 23])
divisors = np.arange(2,int(math.sqrt(np.max(nums)))+1)
nums[np.min(nums[:,None] % divisors[None,:], axis=1) > 0]

所有其他的东西只是为了说明每个步骤在做什么。

这并不是简单的操作,因为它使用了将1D数组广播到2D数组中,但方法应该很清晰。如果您有任何问题,请告诉我。


如果您想进行优化,还有另一种可能性:当前除数是介于2sqrt(max(array))之间的每个数字,但您不需要测试所有这些数字。如果您有一个返回该范围内所有质数的函数,那就足够了。例如,使用@MaxU答案中的primesfrom2to,更快的可能性是:
nums = np.array([17, 18, 19, 20, 21, 22, 23])
# All prime numbers in the range from 2 to sqrt(max(nums))
divisors = primesfrom2to(int(math.sqrt(np.max(nums)))+1)
nums[np.min(nums[:,None] % divisors[None,:], axis=1) > 0]

但它使用了与之前相同的机制,但速度更快了。:-)


0
如果你真的想使用一个过滤器,你可以使用这个:
nums[[i for i in range(len(nums)) if sum([nums[i]%val==0 for val in range(2,nums[i]-1)])==0]]

这是做什么的?
我们使用质数搜索所有索引。
[i for i in range(len(nums)) if sum([nums[i]%val==0 for val in range(2,nums[i]-1)])==0]

这基本上会遍历每个值并检查它是否不能被比它自己小的任何值整除(忽略1)

[i for i in range(len(nums)) #for every index

if sum(#calculate sum of booleans

[nums[i]%val==0 for val in range(2,nums[i]-1)] # check if it is divisble by any value smaller than itself

)==0 #check if the number of divisors is zero

2
为什么会被踩?我知道这段代码很丑陋、未经优化,但它能够正常运行,并满足了提问者的问题。 - JeD

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接