Numpy 256位计算。

Question

我有一些大数用于椭圆曲线加密。目前我使用纯Python进行所有计算，它原生支持通过使用“bignum”整数类型来处理任意大小的数字。但是，出于性能原因，我想对其中一些操作进行矢量化，并希望使用包含256位无符号整数的数组。但据我所知，使用np.longlong（或无符号版本的np.ulonglong）在numpy中整数的最大大小为64位。根据下面的帖子的一个想法是将每个256位数表示为具有四个64位整数的数组：

large = 105951305240504794066266398962584786593081186897777398483830058739006966285013
arr = np.array([large.to_bytes(32,'little')]).view('P')
arr

输出：

array([18196013122530525909, 15462736877728584896, 12869567647602165677,
       16879016735257358861], dtype=uint64)

而我可以通过以下方式将其转换回本地的Python“大数”：

int.from_bytes(arr.tobytes(), 'little')

输出：

105951305240504794066266398962584786593081186897777398483830058739006966285013

然而，在将其转换回“bignum”之前，我想执行一些操作：

arr2 = arr + 1
int.from_bytes(arr2.tobytes(), 'little')

输出：

105951305240504794072543500697971467357257258687906003363414235534976478560982

当然，这与 large + 1 不相等

---

在numpy中处理256位数是否可能？如果不行，有没有合适的解决方法？

- Kevin

---

Python支持使用int整数类型的任意精度整数。什么是bignum？另外，不要忘记语言标签。 - Mad Physicist

你将不得不自己有效地实现所有操作。那时，直接用C编写可能更快/更容易，因为你可以绕过Python的循环，而你需要比如说传递数字。 - Mad Physicist

在快速的 numpy 意义上，“向量化”意味着使用编译好的 numpy 方法。如果这些方法只编译为 int64，那么尝试绕过它们是得不到任何加速的。 - hpaulj

另外，不要使用 arr + 1，而是使用 arr + [0, 0, 0, 1]。 - Mad Physicist

@hpaulj。不过这取决于你如何做。 - Mad Physicist

Answer 1

假设你有几个四字节的数组：

a = np.array([0b0000_1000, 0b1000_1100, 0b0111_1111, 0b1111_0000], dtype=np.uint8)  # [  8, 140, 127, 240]
b = np.array([0b1111_0000, 0b1111_0000, 0b1111_0000, 0b1111_0000], dtype=np.uint8)  # [240, 240, 240, 240]

现在您想将它们相加并得到类似于以下内容的结果。

c = np.array([0b1111_1001, 0b0111_1101, 0b0111_0000, 0b1110_0000], dtype=np.uint8)  # [249, 125, 112, 224]

我故意避免了需要在最高位字节中传递的内容。 在这种情况下，结果是使用计算的。

ai = int.from_bytes(a.tobytes(), byteorder='big')
bi = int.from_bytes(b.tobytes(), byteorder='big')
c = np.frombuffer((ai + bi).to_bytes(4, 'big'), dtype=np.uint8)

那么如何在不将数据转换为 int 的情况下进行操作呢？你不能简单地使用 c = a + b：

[1111_1000, 0111_1100, 0110_1111, 1110_0000]
          +1         +1         +1

请注意，标记位置处的进位都被舍弃了。您可以使用该方法检查进位而不触发溢出：

0xFF - a < b

好消息是，每个字节只能携带一个比特位：0xFF + 0xFF + 1 == 0x1FF。如果和刚刚进位了，你也无法从加1中进位：(0xFF + 0xFF) & 0xFF == 0xFE。这里是适当加法的简单实现：

c = a.copy()
while b.sum():
    carry = np.r_[0xFF - c[1:] < b[1:], False].view(np.uint8)
    c += b
    b = carry

这个循环通常会运行1-2次。在最坏的情况下，它可能会运行N-1次，其中N=256（在您的情况下），或者N=4（在此示例中）。使用更大的数据类型将意味着您可以减少最坏情况下的迭代次数（并且可能加速最佳情况）。唯一的注意事项是您将无法使用...view(np.uint8)技巧来转换布尔掩码。但实际上这是个福音，因为您可以编写

c = a.copy()
carry = np.zeros_like(a)
cbuf = carry.view(bool)[::a.itemsize]
# On big-endian use this instead:
#cbuf = carry.view(bool)[a.itemsize - 1::a.itemsize]
cbuf = cbuf[:-1]
m = np.iinfo(a.dtype).max
while b.sum():
    c += b
    np.less(m - c[1:], b[1:], out=cbuf)
    b = carry

注意，这个版本不仅避免了重新分配进位缓冲区，而且完全不关心a和b的itemsize。

读者练习：

• 正确处理小端序（我的示例按大端序排列单词，与问题不同）
• 添加符号支持（如果您有固定宽度的二进制补码，则可以轻松实现）
• 实际上在N个维度上向量化此操作
• 实现乘法等运算

- Mad Physicist