优化整数系数的名单及其长整型之间的转换(Optimize conversion between li

2019-07-31 13:07发布

站内文章 / 前端开发

56 0

女 | 书童

私信

我试图优化多项式实现矿井。特别是我处理多项式具有系数模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测试），所以我不能因式分解它。

（*）我试过几种方法：

使用numpy阵列，而不是一个列表，并使用乘以numpy.convolve 。它的快速对于n < 2^64但对于非常慢n > 2^64 [也我想避免使用外部库]
使用scipy.fftconvolve 。不适合工作在所有n > 2^64 。
表示系数从一开始就整数（不用每次都将它们转换）。问题是，我不知道一个简单的方法做mod x^r -1操作，而整数转换为系数的列表（这违背了使用这种表示的原因）。

Answer 1:

除非你这样做是为了学习，为什么另起炉灶？一种不同的方法将是一个Python包装写信给一些其他的多项式库或程序，如果这样的包装已经不存在。

尝试PARI / GP。这是出奇的快。我最近写了一个定制的C代码，我花了两天时间写，变成了超过两行PARI / GP脚本更快只有3次。我敢打赌，一个Python代码调用帕里将最终证明是比你在python单独实施任何更快。甚至还有从蟒蛇调用PARI模块： https://code.google.com/p/pari-python/

Answer 2:

你可以尝试使用剩余的数字系统来表示你的多项式的系数。你也将你的系数分成更小的整数因为你现在做的，但你并不需要将它们转换回到一个巨大的整数做乘法或其他操作。这应该没必要再编程工作。

剩余数系统的基本原理是使用模算术数字的唯一表示。周边的RN中的整个理论让你做的小系数您的操作。

编辑：一个简单的例子：

假设你代表你的大系数的RNS与模11和13.你的系数将全部由2个小整数（<11和<13），可以结合到原来的（大）的整数。

假设你的多项式本来是33x²+ 18X + 44。在RNS，该系数将分别是（33 MOD 11，33对13取模），（18模11,18对13取模）和（44模11，44对13取模）=>（0,7），（7,5）和（0,5）。

乘以你的多项式具有恒定然后可以通过每个小系数与常数相乘来完成，并做其模。

假设你乘以3，你的系数将成为（0,21模13）=（0,8），（21模11,15模13）=（10,2）和（0模11,15模13）= （0,2）。目前还没有必要系数转换回自己的大整数。

要检查我们的乘法工作过，我们可以在新的系数转换回其大代表。这就需要“解决”各组系数作为一个模块化系统。对于第一个系数（0.8），我们将需要解决在x mod 11 = 0，在x mod 13 = 8，这应该不是太难实现。在这个例子中可以看出，X = 99是一个有效的解决方案（模13 * 11）

然后，我们得到99x²+ 54X + 132，正确的乘多项式。与其他多项式乘法是类似（但需要你乘相互系数以成对的方式）。这同样适用于除。

对于你的使用情况，您可以根据系数想要次数或它们的大小选择你ñ。

Answer 3:

如何直接实现任意精度的整数多项式作为numpy的数组列表？

让我来解释一下：说你的多项式_ΣP A _P X _[ ^页。如果大整数的_p可表示为_p _=ΣK A _P，K 2 ^{64 10k，}则^第 k ^个 numpy的阵列将包含64位int甲_P，K在位置p。

你可以根据你的问题的结构选择密集或稀疏矩阵。

实现加法和标量运算只是一个矢量化的BIGNUM执行相同操作的问题。

AB _{=ΣP，K，P '中，k'} A _P，K乙_{P '中，k'} 2 ^{64（K + K '）X} ^{P + P'：}乘法可以如下进行处理。因此，一个天真的实现与密集阵列可能会导致数_^64（N）2调用numpy.convole或scipy.fftconvolve 。

模操作应该很容易实现，因为它是左边项的线性函数和右手项具有小系数。

编辑这里有一些更多的解释

代替表示作为多项式任意精度号码的列表的（自己表示为64位的“数字”列表），转置的表示，使得：

您的矩阵表示为数组列表
^第 k ^个数组包含^第 k ^个各系数的“数字”

如果您仅有几个系数非常大，则阵列将主要在他们0，因此它可能使用稀疏阵列是值得的。

呼叫的_P，K ^第 k ^个 ^第 p ^个系数的位数。

注意大整数表示的比喻：那里有一个很大的整数将表示为

X =Σ _的K×2 _k ^64千

您的多项式A在方法一样代表

_的 2 ^至64且_k _{_=Σpk信息A，K} X ^p _的 =Σ甲

要实现此外，您只需假装你的阵列的列表是简单的数字列表，并实施另外像往常一样为大整数（注意更换if then条件句由numpy.where ）。

要实现倍增，你会发现你需要做日志^{_64（N）2个}多项式乘法。

为贯彻落实系数模操作，是再次翻译在一个大整数模运算的简单情况。

通过与小系数的多项式取模，使用此操作的线性度：

甲MOD（X ^的R - 1）= MOD（Σ _约 2 ^至64 _的 ^A）（X ^的R - 1）
=Σ _一个 2 ^到_64（ _模式 A（X ^的R - 1））

Answer 4:

我找到了一种方法来优化转换，尽管我还是希望有人能帮助我更加完善他们，并希望找到一些其他聪明的主意。

基本上，什么是错的那些功能是，他们有某种二次内存分配行为，打包整数时，或拆包时。（见这个职位吉多·范罗苏姆对这种行为的其他例子）。

我意识到这一点之后，我已经决定放弃与分而治之原则一试，我已经取得了一些成果。我简单地划分阵列中的两个部分，分别将它们转换并最终加入的结果（稍后我会尝试使用类似的迭代版本f5在Rossum的信息[编辑：它似乎没有要快得多]）。

修改后的功能：

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塞巴斯蒂安指出AKS算法得到了许多改进，所以试图执行旧版本的算法总是会导致一个非常缓慢的程序。这并不排除这样的事实，如果你已经有了一个良好的执行AKS你可以加速这一过程提高了多项式。
如果您有兴趣系数模小n （读：字大小号），你不介意的外部依赖性，然后去numpy和使用numpy.convolve或scipy.fftconvolve的乘法。这将是比什么都可以写快得多。不幸的是如果n不是字的大小，你不能使用scipy.fftconvolve可言，也numpy.convolve变得慢如地狱。
如果你没有做模运算（上系数和多项式），则可能使用ZBDDs是一个好主意（如哈罗德指出的），即使我不答应了瞩目的成绩[尽管我认为这是真的很有趣，你应该阅读港区的纸。
如果你没有做对，然后可能使用的RNS表示系数模运算，如原产地说，是一个好主意。然后，你可以合并多个numpy阵列有效地运作。
如果你想要一个纯Python实现与系数多项式的模大n ，然后我的解决方案似乎是最快的。虽然我并没有试图实现的Python系数的阵列之间的FFT乘（这可能会更快）。