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

Question

我正在尝试优化我的多项式实现.特别是我正在处理具有模数n(可能是>2^64)的多项式,并以形式x^r - 1(r是< 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 van Rossum 的这篇文章）。

当我意识到这一点后，我决定尝试一下分而治之的原则，并且我已经取得了一些成果。我只是将数组分成两部分，分别转换它们并最终加入结果（稍后我将尝试使用类似于f5Rossum 帖子中的迭代版本[编辑：它似乎并没有快得多]）。

修改后的功能：

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 测试的人来说是有用的。

• 正如 JF Sebastian 指出的，AKS 算法得到了许多改进，因此尝试实现旧版本的算法总是会导致程序非常慢。这并不排除这样一个事实：如果您已经很好地实现了 AKS，则可以加快改进多项式的速度。
• 如果您对以小数为模n（阅读：字大小的数字）的系数感兴趣并且不介意外部依赖性，那么请使用numpyornumpy.convolve进行scipy.fftconvolve乘法。它会比你写的任何东西都要快得多。不幸的是，如果不是字大小，你根本n无法使用，而且也会变得非常慢。scipy.fftconvolvenumpy.convolve
• 如果您不必进行模运算（对系数和多项式），那么使用 ZBDD 可能是一个好主意（正如 harold 所指出的），即使我不保证惊人的结果[尽管我认为它是真的很有趣，你应该阅读 Minato 的论文]。
• 如果您不必对系数进行模运算，那么使用 RNS 表示（如 Origin 所述）可能是一个好主意。然后你可以组合多个numpy数组来高效运行。
• 如果你想要一个系数模大的多项式的纯Python实现n，那么我的解决方案似乎是最快的。即使我没有尝试在 python 中的系数数组之间实现 fft 乘法（这可能更快）。