8

我正在尝试优化我的多项式实现。特别是我正在处理系数模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不将整数转换为系数列表的情况下进行操作的简单方法（这违背了使用这种表示的原因）。
4

4 回答 4

2

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

尝试 PARI/GP。它出奇的快。我最近编写了一个自定义 C 代码，花了我两天时间编写，结果只比两行 PARI/GP 脚本快 3 倍。我敢打赌，调用 PARI 的 python 代码最终会比你单独在 python 中实现的更快。甚至还有一个从 python 调用 PARI 的模块： https ://code.google.com/p/pari-python/

于 2012-09-20T14:27:32.470 回答
2

您可以尝试使用残差数系统来表示多项式的系数。您还可以像现在一样将系数拆分为较小的整数，但您不需要将它们转换回大整数来进行乘法或其他操作。这应该不需要太多的重新编程工作。

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

编辑： 一个简单的例子：

假设您用模数 11 和 13 在 RNS 中表示您的大系数。您的系数将全部由 2 个小整数（<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 和 x mod 13 = 8。这应该不太难实现。在此示例中，您可以看到 x=99 是一个有效的解决方案（模 13*11）

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

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

于 2012-10-18T21:28:13.797 回答
2

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

让我解释一下：假设你的多项式是 Σ p A p X p。如果大整数 A p可以表示为 A p = Σ k A p,k 2 64 k那么第 k个numpy 数组将包含位置 p 处的 64 位 int A p,k。

您可以根据问题的结构选择密集或稀疏数组。

实现加法和标量操作只是向量化相同操作的 bignum 实现的问题。

乘法可以按如下方式处理： AB = Σ p,k,p',k' A p,k B p',k' 2 64(k+k') X p+p'。因此，使用密集数组的简单实现可能会导致 log 64 (n) 2调用numpy.convoleor scipy.fftconvolve。

模运算应该很容易实现，因为它是左侧项的线性函数，而右侧项的系数很小。

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

不是将多项式表示为任意精度数字的列表（它们本身表示为 64 位“数字”列表），而是转置表示，以便：

  • 您的多项式表示为数组列表
  • 第 k个数组包含每个系数的第 k个“数字”

如果只有少数系数非常大，那么数组中的大部分都是 0，因此使用稀疏数组可能是值得的。

将 A p,k称为第 p个系数的第k个数字。

请注意与大整数表示的类比：其中大整数将表示为

x = Σ k x k 2 64 k

您的多项式 A 的表示方式与

A = Σ k A k 2 64 k A k = Σ k A p,k X p

要实现加法，您只需假设您的数组列表是一个简单数字列表，并像往常一样为大整数实现加法（注意用 替换if then条件numpy.where）。

要实现乘法，您会发现需要进行 log 64 (n) 2多项式乘法。

为了对系数进行模运算，又是一个简单的情况，即对大整数进行模运算。

要对具有小系数的多项式取模，请使用此操作的线性：

A mod (X r - 1) = (Σ k A k 2 64 k ) mod (X r - 1)

= Σ k 2 64 k (A k mod (X r - 1))

于 2012-10-23T23:33:37.730 回答
1

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

基本上，这些函数的问题在于它们在打包整数或解包时具有某种二次内存分配行为。（有关此类行为的另一个示例，请参阅 Guido van Rossum 的这篇文章）。

在我意识到这一点后，我决定尝试使用 Divide et Impera 原则，并获得了一些结果。我只是将数组分成两部分，分别转换它们并最终加入结果（稍后我将尝试使用类似于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（阅读：字大小数）并且您不介意外部依赖关系，那么请numpy使用numpy.convolveorscipy.fftconvolve进行乘法运算。它比你能写的任何东西都要快得多。不幸的是，如果n不是单词大小，您根本无法使用scipy.fftconvolve，而且numpy.convolve会变得慢得要命。
  • 如果您不必进行模运算（在系数和多项式上），那么使用 ZBDD 可能是一个好主意（正如 harold 所指出的那样），即使我不承诺会取得惊人的结果[即使我认为这是真的很有趣，你应该阅读 Minato 的论文]。
  • 如果您不必对系数进行模运算，那么可能使用 RNS 表示，如 Origin 所述，是一个好主意。然后您可以组合多个numpy数组以高效运行。
  • 如果您想要系数模 a big 的多项式的纯 python 实现n，那么我的解决方案似乎是最快的。即使我没有尝试在python中的系数数组之间实现fft乘法（这可能更快）。
于 2012-10-21T12:02:56.390 回答