整数系数表及其长整数表示之间的优化转换

Question

我在试图优化我的多项式实现。特别是，我正在处理系数模n(可能是>2^64)的多项式，以及x^r - 1(r is < 2^64)形式的多项式。目前，我将系数表示为整数(*)的列表，并以最简单的方式实现了所有基本操作。

我希望指数和乘法尽可能快，为了得到这一点，我已经尝试过不同的方法。我目前的方法是把系数的列表转换成巨大的整数，乘以整数，然后解压缩系数。

问题是包装和拆包需要很长时间。

那么，有什么方法可以改善我的“打包/解压”功能吗？

def _coefs_to_long(coefs, window):
    '''Given a sequence of coefficients *coefs* and the *window* size return a
    long-integer representation of these coefficients.
    '''

    res = 0
    adder = 0
    for k in coefs:
        res += k << adder
        adder += window
    return res
    #for k in reversed(coefs): res = (res << window) + k is slower


def _long_to_coefs(long_repr, window, n):
    '''Given a long-integer representing coefficients of size *window*, return
    the list of coefficients modulo *n*.
    '''

    mask = 2**window - 1
    coefs = [0] * (long_repr.bit_length() // window + 1)
    for i in xrange(len(coefs)):
        coefs[i] = (long_repr & mask) % n
        long_repr >>= window

    # assure that the returned list is never empty, and hasn't got an extra 0.
    if not coefs:
        coefs.append(0)
    elif not coefs[-1] and len(coefs) > 1:
        coefs.pop()

    return coefs

请注意，我没有选择n，它是来自用户的输入，我的程序想要证明它的原素性(使用AKS测试)，所以我不能将它分解。

(*)我尝试过几种方法：

1. 使用numpy数组代替列表，并使用numpy.convolve进行乘法。它对于n < 2^64来说是快速的，但是对于n > 2^64来说非常慢，我也想避免使用外部库
2. 使用scipy.fftconvolve。对n > 2^64根本不起作用。
3. 从一开始就将系数表示为整数(而不是每次转换它们)。问题是，如果不将整数转换为系数列表(这就违背了使用这种表示的原因)，我不知道一种简单的方法来执行mod x^r -1操作。

Answer 1

我找到了一种优化转换的方法，尽管我仍然希望有人能帮助我更多地改进它们，并希望找到其他一些聪明的想法。

基本上，这些函数的错误之处在于，它们具有某种二次内存分配行为，在打包整数或解压整数时。(关于这类行为的另一个例子，请参阅Guido的这帖子)。

当我意识到这一点后，我决定尝试一下除法原理，并得到了一些结果。我只需将数组分成两部分，分别转换它们，并最终加入结果(稍后，我将尝试使用类似于Rossum的postedit中的f5的迭代版本:它似乎不会更快)。

修改后的职能：

def _coefs_to_long(coefs, window):
    """Given a sequence of coefficients *coefs* and the *window* size return a
    long-integer representation of these coefficients.
    """

    length = len(coefs)
    if length < 100:
        res = 0
        adder = 0
        for k in coefs:
            res += k << adder
            adder += window
        return res
    else:
        half_index = length // 2
        big_window = window * half_index
        low = _coefs_to_long(coefs[:half_index], window)
        high = _coefs_to_long(coefs[half_index:], window)
        return low + (high << big_window)


def _long_to_coefs(long_repr, window, n):
    """Given a long-integer representing coefficients of size *window*, return
    the list of coefficients modulo *n*.
    """

    win_length = long_repr.bit_length() // window
    if win_length < 256:
        mask = 2**window - 1
        coefs = [0] * (long_repr.bit_length() // window + 1)
        for i in xrange(len(coefs)):
            coefs[i] = (long_repr & mask) % n
            long_repr >>= window

        # assure that the returned list is never empty, and hasn't got an extra 0.
        if not coefs:
            coefs.append(0)
        elif not coefs[-1] and len(coefs) > 1:
            coefs.pop()

        return coefs
    else:
        half_len = win_length // 2
        low = long_repr & (((2**window) ** half_len) - 1)
        high = long_repr >> (window * half_len)
        return _long_to_coefs(low, window, n) + _long_to_coefs(high, window, n) 

其结果是：

>>> import timeit
>>> def coefs_to_long2(coefs, window):
...     if len(coefs) < 100:
...         return coefs_to_long(coefs, window)
...     else:
...         half_index = len(coefs) // 2
...         big_window = window * half_index
...         least = coefs_to_long2(coefs[:half_index], window) 
...         up = coefs_to_long2(coefs[half_index:], window)
...         return least + (up << big_window)
... 
>>> coefs = [1, 2, 3, 1024, 256] * 567
>>> # original function
>>> timeit.timeit('coefs_to_long(coefs, 11)', 'from __main__ import coefs_to_long, coefs',
...               number=1000)/1000
0.003283214092254639
>>> timeit.timeit('coefs_to_long2(coefs, 11)', 'from __main__ import coefs_to_long2, coefs',
...               number=1000)/1000
0.0007998988628387451
>>> 0.003283214092254639 / _
4.104536516782767
>>> coefs = [2**64, 2**31, 10, 107] * 567
>>> timeit.timeit('coefs_to_long(coefs, 66)', 'from __main__ import coefs_to_long, coefs',...               number=1000)/1000

0.009775240898132325
>>> 
>>> timeit.timeit('coefs_to_long2(coefs, 66)', 'from __main__ import coefs_to_long2, coefs',
...               number=1000)/1000
0.0012255229949951173
>>> 
>>> 0.009775240898132325 / _
7.97638309362875

正如您所看到的，这个版本提供了相当快的转换速度，从4到8的转换速度快了一倍(输入越大，速度越快)。第二个函数也得到了类似的结果：

>>> import timeit
>>> def long_to_coefs2(long_repr, window, n):
...     win_length = long_repr.bit_length() // window
...     if win_length < 256:
...         return long_to_coefs(long_repr, window, n)
...     else:
...         half_len = win_length // 2
...         least = long_repr & (((2**window) ** half_len) - 1)
...         up = long_repr >> (window * half_len)
...         return long_to_coefs2(least, window, n) + long_to_coefs2(up, window, n)
... 
>>> long_repr = coefs_to_long([1,2,3,1024,512, 0, 3] * 456, 13)
>>> # original function
>>> timeit.timeit('long_to_coefs(long_repr, 13, 1025)', 'from __main__ import long_to_coefs, long_repr', number=1000)/1000
0.005114212036132813
>>> timeit.timeit('long_to_coefs2(long_repr, 13, 1025)', 'from __main__ import long_to_coefs2, long_repr', number=1000)/1000
0.001701267957687378
>>> 0.005114212036132813 / _
3.006117885794327
>>> long_repr = coefs_to_long([1,2**33,3**17,1024,512, 0, 3] * 456, 40)
>>> timeit.timeit('long_to_coefs(long_repr, 13, 1025)', 'from __main__ import long_to_coefs, long_repr', number=1000)/1000
0.04037192392349243
>>> timeit.timeit('long_to_coefs2(long_repr, 13, 1025)', 'from __main__ import long_to_coefs2, long_repr', number=1000)/1000
0.005722791910171509
>>> 0.04037192392349243 / _
7.0545853417694

我试图避免在第一个函数中更多地重新分配内存，在开始和结束索引之间传递，并避免切片，但事实证明，对于小输入来说，这大大降低了函数的速度，而对于实际情况的输入，则稍微慢了一点。也许我可以尝试混合他们，尽管我不认为我会获得更好的结果。

我在最后一段时间编辑了我的问题，因此有些人给了我一些与我最近所要求的不同目标的建议。我认为重要的是要澄清一些不同来源在评论和答案中指出的结果，这样它们就可以对其他想要实现快速多项式和AKS测试的人有用。

• 正如J.F. Sebastian指出的那样，AKS算法得到了许多改进，因此尝试实现该算法的旧版本总是会导致一个非常慢的程序。这并不排除这样一个事实:如果您已经很好地实现了AKS，您可以加快它改进多项式。
• 如果您感兴趣的系数模一个小n(读:字大小的数字)，你不介意外部依赖，然后选择numpy，并使用numpy.convolve或scipy.fftconvolve的乘法。它将比你能写的任何东西都快得多。不幸的是，如果n不是字大小，那么您就根本不能使用scipy.fftconvolve，而且numpy.convolve也变得非常慢。
• 如果你不需要做模运算(对系数和多项式)，那么使用ZBDDs可能是个好主意(正如哈罗德所指出的)，尽管我没有承诺会有惊人的结果，尽管我认为它真的很有趣，你应该读一下Minato的论文。
• 如果您不需要对系数进行模块化操作，那么很可能使用RNS表示，正如Origin所述，这是一个好主意。然后，您可以组合多个numpy数组来有效地操作。
• 如果您想要一个系数模的多项式的纯python实现-一个大的n，那么我的解决方案似乎是最快的。即使我没有尝试实现python中系数数组之间的fft乘法(这可能更快)。

Answer 2

除非你这么做是为了学习，否则为什么要重新发明轮子呢？另一种方法是向其他多项式库或程序编写python包装器，如果这种包装器还不存在的话。

试试帕里/GP。令人惊讶的快。我最近写了一个定制的C代码，花了我两天的时间编写，结果只比两行PARI/GP脚本快3倍。我敢打赌，调用PARI的python代码将比单独在python中实现的代码更快。甚至还有一个模块可以从python调用PARI：https://code.google.com/p/pari-python/

Answer 3

您可以尝试使用剩余数系统来表示多项式的系数。你也会像现在一样把你的系数分解成更小的整数，但是你不需要把它们转换回一个巨大的整数来进行乘法或其他运算。这不需要太多的重新编程工作。

残数系统的基本原理是用模算法对数进行唯一表示。围绕RNS的整个理论允许您对小系数进行操作。

编辑：--一个快速示例：

假设你在一个RNS中用模11和13表示你的大系数，你的系数都由两个小整数组成(<11和<13)，它们可以组合成原始的(大)整数。

假设你的多项式原来是33x²+18x+44。RNS的相关系数分别为(33 mod 11，33 mod 13)，(18 mod 11,18 mod 13)和(44 mod 11，44 mod 13)=>(0,7)，(7,5)和(0,5)。

然后，用一个常数乘以你的多项式，就可以用这个常数乘以每个小系数，然后在它上做模运算。

假设你乘以3，你的系数将变成(0,21 mod 13)= (0 ,8)，(21 mod 11,15 mod 13)=(10,2)和(0 mod 11,15 mod 13)=(0,2)。没有必要将系数转换回它们的大整数。

为了检查我们的乘法是否有效，我们可以将新的系数转换回它们的大表示。这需要将每一组系数“求解”为一个模块系统。对于第一个系数(0, 8 )，我们需要求解x mod 11=0和need 13 =8，这应该不会太难实现。在这个例子中，您可以看到x=99是一个有效的解决方案(模块化13*11)

然后我们得到99x²+54x+132，正确的乘多项式。与其他多项式相乘是相似的(但要求你以成对的方式将系数相乘)。加法也是如此。

对于用例，可以根据所需系数的数目或它们的大小来选择n。