Lucas序列实现

Question

您能否建议改进卢卡斯序列's的实现：

from mpmath.libmp import bitcount as _bitlength

def _int_tuple(*i):
    return tuple(int(_) for _ in i)

def _lucas_sequence(n, P, Q, k): 
    """Return the modular Lucas sequence (U_k, V_k, Q_k).
    Given a Lucas sequence defined by P, Q, returns the kth values for
    U and V, along with Q^k, all modulo n.
    """
    D = P*P - 4*Q
    if n < 2:
        raise ValueError("n must be >= 2")
    if k < 0:
        raise ValueError("k must be >= 0")
    if D == 0:
        raise ValueError("D must not be zero")

    if k == 0:
        return _int_tuple(0, 2, Q)
    U = 1
    V = P
    Qk = Q
    b = _bitlength(k)
    if Q == 1:
        # For strong tests
        while b > 1:
            U = (U*V) % n
            V = (V*V - 2) % n
            b -= 1
            if (k >> (b - 1)) & 1:
                t = U*D
                U = U*P + V
                if U & 1:
                    U += n
                U >>= 1
                V = V*P + t
                if V & 1:
                    V += n
                V >>= 1
    elif P == 1 and Q == -1:
        # For Selfridge parameters
        while b > 1:
            U = (U*V) % n
            if Qk == 1:
                V = (V*V - 2) % n
            else:
                V = (V*V + 2) % n
                Qk = 1
            b -= 1
            if (k >> (b-1)) & 1:
                t = U*D
                U = U + V
                if U & 1:
                    U += n
                U >>= 1
                V = V + t
                if V & 1:
                    V += n
                V >>= 1
                Qk = -1
    else:
        # The general case with any P and Q
        while b > 1:
            U = (U*V) % n
            V = (V*V - 2*Qk) % n
            Qk *= Qk
            b -= 1
            if (k >> (b - 1)) & 1:
                t = U*D
                U = U*P + V
                if U & 1:
                    U += n
                U >>= 1
                V = V*P + t
                if V & 1:
                    V += n
                V >>= 1
                Qk *= Q
            Qk %= n
    U %= n
    V %= n
    return _int_tuple(U, V, Qk)

Answer 1

if D == 0: raise ValueError("D must not be zero")

D不是一个参数，所以在异常错误消息中使用它的名称似乎有点奇怪。

此外，考虑到您所提供的引用只为D > 0定义了卢卡斯序列，您很高兴拥有D < 0似乎有点奇怪。这可能是值得补充评论说，你是支持消极的决定因素。

# For strong tests

我由此推断，你正在使用的序列卢卡斯或卢卡斯-莱默的原始性测试，但这一评论似乎有点不合适。

有三个案子似乎有点过火了。你对代码进行了剖析，看看特例有多大的不同吗？您似乎主要节省了1的一些乘法运算，而且我不认为维护工作是值得的，如果它只是函数的2%的加速比，那么无论如何也不会成为瓶颈。

B= _bitlength(k) .B> 1：.B -= 1 if (k >> (b - 1)) & 1：.

我必须对这段代码进行测试，才能确信您并没有将位向后移动(也就是说，您确实希望从最重要的位到最不重要的位到最不重要的位)。这看起来确实是正确的，但如果不给出一个解释为什么是正确的评论，那肯定太聪明了。

t = U\*D U = U\*P + V if U & 1: U += n U >>= 1 V = V\*P + t if V & 1: V += n V >>= 1

Python最伟大的特性之一，我希望有更多的语言是并行赋值。我将丢弃临时变量，并将U和V之间的并行性公开为

                U, V = U*P + V, V*P + U*D
                if U & 1:
                    U += n
                if V & 1:
                    V += n
                U, V = U >> 1, V >> 1

U %= n V %= n return \_int\_tuple(U, V, Qk)

在我看来，这并不是把太多的东西写成这样(按照格雷弗的建议)。

    return (U % n, V % n, Qk)

Answer 2

您的_int_tuple函数似乎没有任何实际用途。U和V总是整数(您将它们初始化为1，然后只对它们执行算术)。您的程序对Q的具体操作是什么--例如，字符串很难理解，但我认为它会遇到else块(Q既不是1也不是-1)，然后尝试执行Qk = Q，然后是Qk *= Qk，除了数字类型之外，几乎所有东西都失败了。

我也看不到代码的任何部分，它们会产生一个float (尽管我可能忽略了一个)。

所以我会去掉这个函数，只返回一个元组，看看它能做什么。

Answer 3

总的来说，在我看来很不错。我会进一步澄清这些评论。Q=1通常用于超强Lucas测试，P=1，Q=-1用于标准或强卢卡斯测试的大约一半情况(其余的情况有不同的Q)。关于各种测试的一些信息可以在我的伪统计页面上看到。

如果n是偶数会发生什么？例如，Lucas_5(6,1) = (1189,6726)。国防部1000应产量(189,726)。但例行公事给了我们(689,726)。这是一个工件的一半使用的方法。我的解决方案是在这种情况下使用另一种更慢的方法，因为它很少被使用(从来没有通过标准的原始性测试，但还有更多的用途)。

对于D=0，我们可以返回一个解决方案。见维基百科页面。

检查P和/或Q和/或D超过n的情况。从快速测试中可以看出，返回的Q_k值超出了k={0,1}的范围。预先进行检查和模组可能会提高性能。

考虑k.bit_length()而不是使用libmp。取决于你，但我喜欢更少的依赖。不过，我不知道两者之间的权衡。

对于Q=1情况(由特别强的Lucas测试使用)，还有另一个可以更快的方法。如果D可以倒模n，那么我们就可以用每比特2多个V_k来计算U_k，然后用逆运算来计算U_k。我不想增加更多的情况，当然，您需要进行基准测试，看看它是否更适合您的基础设施和使用。