如何使用NumPy查找emirp编号？

elcex8rz 于 2023-06-23 发布在其他

关注(0)|答案(2)|浏览(142)

我想找出其数字颠倒也是素数的素数，包括回文素数，我也希望在其他基中找到具有这种性质的素数，而不仅仅是十进制。质数的数字反转也是质数，称为emirp，关于它的更多信息在link上。
我已经实现了一个给出正确输出的解决方案，如下所示。我想用NumPy实现同样的事情。

import numpy as np
from itertools import cycle

def prime_wheel_sieve(n: int) -> np.ndarray:
    wheel = cycle([4, 2, 4, 2, 4, 6, 2, 6])
    primes = np.ones(n + 1, dtype=bool)
    primes[:2] = False
    for square, step in ((4, 2), (9, 6), (25, 10)):
        primes[square::step] = False

    k = 7
    while (square := k * k) <= n:
        if primes[k]:
            primes[square :: 2 * k] = False

        k += next(wheel)
    return np.where(primes)[0]

DIGITS = {
    **{n: chr(n + 48) for n in range(10)},
    **{n: chr(n + 87) for n in range(10, 36)}
}

def to_base_str(n: int, base: int) -> str:
    if not 2 <= base <= 36:
        raise ValueError

    s = []
    while n:
        n, d = divmod(n, base)
        s.append(DIGITS[d])

    return "".join(s[::-1])

def to_base(n: int, base: int) -> str | tuple:
    if base < 2:
        raise ValueError

    if base <= 36:
        return to_base_str(n, base)

    s = []
    while n:
        n, d = divmod(n, base)
        s.append(d)

    return tuple(s[::-1])

def from_base(s: str | tuple, base: int) -> int:
    if base < 2:
        raise ValueError
    if base <= 36:
        return int(s, base)
    return sum(d * base ** i for i, d in enumerate(s[::-1]))

def emirp(n: int, base: int = 10) -> np.ndarray:
    primes = prime_wheel_sieve(n)
    primes = primes[primes > base]
    elist = [to_base(p, base) for p in primes]
    eset = set(elist)
    return primes[[i for i, e in enumerate(elist) if e[::-1] in eset]]

示例：

In [90]: emirp(1000, 2)
Out[90]:
array([  3,   5,   7,  11,  13,  17,  23,  29,  31,  37,  41,  43,  47,
        53,  61,  67,  71,  73,  83,  97, 101, 107, 113, 127, 131, 151,
       163, 167, 173, 181, 193, 197, 199, 223, 227, 229, 233, 251, 257,
       263, 269, 277, 283, 307, 313, 331, 337, 349, 353, 359, 373, 383,
       409, 421, 431, 433, 443, 449, 461, 463, 479, 487, 491, 503, 509,
       521, 571, 577, 599, 601, 617, 619, 631, 643, 653, 661, 677, 683,
       691, 701, 709, 727, 739, 757, 773, 797, 821, 823, 827, 839, 853,
       857, 881, 883, 907, 911, 937, 941, 947, 953, 967])

In [91]: emirp(1000, 10)
Out[91]:
array([ 11,  13,  17,  31,  37,  71,  73,  79,  97, 101, 107, 113, 131,
       149, 151, 157, 167, 179, 181, 191, 199, 311, 313, 337, 347, 353,
       359, 373, 383, 389, 701, 709, 727, 733, 739, 743, 751, 757, 761,
       769, 787, 797, 907, 919, 929, 937, 941, 953, 967, 971, 983, 991])

In [92]: emirp(1000, 5)
Out[92]:
array([  7,  11,  13,  17,  19,  23,  29,  31,  37,  41,  47,  59,  61,
        67,  71,  73,  79,  83,  89,  97, 101, 103, 107, 109, 113, 127,
       131, 149, 151, 157, 163, 167, 191, 193, 211, 223, 227, 229, 233,
       239, 251, 257, 269, 271, 277, 281, 293, 317, 331, 337, 347, 349,
       353, 359, 367, 379, 397, 421, 431, 439, 443, 449, 457, 461, 463,
       479, 491, 503, 521, 523, 547, 563, 571, 577, 599, 601, 613, 617,
       619, 631, 701, 751, 761, 881, 911])

更新

我刚刚发现np.base_repr，它做的事情几乎和我的自定义函数完全一样，除了它不像int那样处理大于36的碱基，而且它不接受数组作为输入。所以它在这里是无用的。我之所以能找到它，是因为我在谷歌上搜索了一个与它完全无关的主题。

numpy

来源：https://stackoverflow.com/questions/76495180/how-to-find-emirp-numbers-using-numpy

2条答案

按热度按时间

y4ekin9u1#

将数字转换为字符串以反转它们的数字，然后将反转的字符串转换回整数，这样做效率不高。相反，您应该在整数空间中完成所有操作，如下所示：

def reverse_digits(value, base):
    result = 0
    while value > 0:
        right_digit = value % base
        value = value // base
        result = result * base + right_digit 
    return result

这比使用字符串要少得多。
在循环中应用上面的函数会相当慢。一个简单的解决方案是Numba，一个JIT编译器库。为了有效地在整个数组上应用elementwise，使用numba.vectorize创建一个NumPy“ufunc”：

import numba

@numba.vectorize(nopython=True)
def reverse_digits(value, base):
    # ...

你将得到一个以本机速度运行的编译函数，就像你用C写的一样。所以你可以像这样过滤整个primes的NumPy数组：

primes[np.isin(reverse_digits(primes, 10), primes)]

在我的机器上举一个简单的例子，Numba提供了大约200倍的加速。
最后，不需要反转数字：你只需要检查value == reverse_digits(value, base)。这可能允许进一步的加速。

赞(0）回复(0）举报 2023-06-23

vmjh9lq92#

余数和地板除法结果可以通过单个函数调用获得：divmod，这将消除使用//和%的计算冗余。
numba.vectorize可以进一步调整更多的标志，使其更有效。
primes[primes == reverse_digits(primes, 10)]的意思是寻找其数位反转等于其自身的素数，换句话说就是回文素数。
我打算找到在给定基数中的数字反转也是素数的素数，所以我必须通过成员检查来过滤素数，即我需要找到素数数组中的数字反转。
primes[np.isin(reverse_digits(primes, base), primes)]做我想做的。
但是np.isin的效率非常低，更好的方法是将primes数组转换为set，同时使用arr.tolist()绕过numpy，然后使用列表解析来构建掩码。

import numba
import numpy as np
from itertools import cycle

def prime_wheel_sieve(n: int) -> np.ndarray:
    wheel = cycle([4, 2, 4, 2, 4, 6, 2, 6])
    primes = np.ones(n + 1, dtype=bool)
    primes[:2] = False
    for square, step in ((4, 2), (9, 6), (25, 10)):
        primes[square::step] = False

    k = 7
    while (square := k * k) <= n:
        if primes[k]:
            primes[square :: 2 * k] = False

        k += next(wheel)

    return np.where(primes)[0]

@numba.vectorize(nopython=True, cache=True, fastmath=True, forceobj=False)
def reverse_digits(value: int | np.ndarray, base: int) -> np.ndarray:
    result = 0
    while value:
        value, remainder = divmod(value, base)
        result = result * base + remainder
    return result

def emirp(n: int, base: int = 10) -> np.ndarray:
    primes = prime_wheel_sieve(n)
    flipped = reverse_digits(primes, base).tolist()
    prime_set = set(primes.tolist())
    return primes[[x in prime_set for x in flipped]]

一位数的素数是回文素数。因此，包括它们是正确的，并将消除不必要的计算。
但实际上，np.isin对于较大的输入更好地扩展，我没有考虑过这一点，结果对于小输入Python获胜，因为有NumPy开销，但对于较大的输入Python失败，NumPy工作得更好。

def emirp(n: int, base: int = 10) -> np.ndarray:
    primes = prime_wheel_sieve(n)
    flipped = reverse_digits(primes, base)
    return primes[np.isin(flipped, primes)]

赞(0）回复(0）举报 2023-06-23

我来回答

如何使用NumPy查找emirp编号？

更新

2条答案

相关问题

热门标签

最新问答