Numpy 256 位计算答案

【问题标题】：Numpy 256-bit computingNumpy 256 位计算
【发布时间】：2021-07-01 18:25:55
【问题描述】：

我有一些用于椭圆曲线密码学的大数。到目前为止，我一直在使用“bignum”整数类型在原生支持任意大小的数字的普通 python 中进行所有计算。但是，我想对其中一些操作进行矢量化（出于性能原因），并想使用一个包含 256 位无符号整数的数组。但是AFAIK通过使用np.longlong（或np.ulonglong对于无符号版本），numpy中整数的最大大小是64位。 the following post 的一个想法是将每个 256 位数字表示为一个包含四个 64 位整数的数组：

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

输出：

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

我可以通过执行以下操作将其转换回本机 python“bignum”：

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

输出：

105951305240504794066266398962584786593081186897777398483830058739006966285013

但是，我想在将其转换回“bignum”之前进行一些操作：

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

输出：

105951305240504794072543500697971467357257258687906003363414235534976478560982

当然不等于large + 1

是否可以在 numpy 中处理 256 位数字？如果没有，是否存在任何合适的解决方法？

【问题讨论】：

Python 支持使用 int 整数类型的任意精度整数。 bignum 是什么？另外，不要忘记语言标签
您必须自己有效地实施所有操作。到那时，直接用 C 语言编写可能会更快/更容易，因为您可以绕过 python 的循环，例如，您需要这样做。结转数字。
"vectorize" 在 fast(er) numpy 意义上意味着使用已编译的 numpy 方法。如果这些方法仅针对 int64 编译，则尝试绕过它们将不会获得任何速度。
另外，不要做arr + 1，而是做arr + [0, 0, 0, 1]
@hpaulj。取决于你是怎么做的

标签： python numpy

【解决方案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 完全无关。

读者练习：

正确处理 little-endian 顺序（我的示例中的单词按 big-endian 顺序排列，与问题不同）
添加符号支持（如果您有固定宽度的二进制补码，则开箱即用）
实际上是跨 N 维向量化此操作
实现乘法等运算

【讨论】：