如何使用NumPy实现多进制emirp数（可逆素数）查找？

我需要查找满足数字反转后仍是素数的素数（包括回文素数），且希望支持十进制以外的其他进制。这类素数被称为emirp。

我已经实现了一个能输出正确结果的解决方案，代码如下，现在希望用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的进制，也无法接收数组作为输入，所以在这个场景里派不上用场。我是在搜索完全无关的主题时偶然发现它的。

内容的提问来源于stack exchange，提问作者Ξένη Γήινος