有效地计算卢卡斯序列答案

【问题标题】：Calculating Lucas Sequences efficiently有效地计算卢卡斯序列
【发布时间】：2014-05-19 15:45:54
【问题描述】：

我正在实现 p+1 分解算法。为此，我需要计算由以下定义的卢卡斯序列的元素：

(1) x_0 = 1, x_1 = a
(2) x_n+l  =  2 * a * x_n - x_n-l

我递归地实现了它（C#），但它对于更大的索引效率低下。

static BigInteger Lucas(BigInteger a, BigInteger Q, BigInteger N)
    {
        if (Q == 0)
            return 1;
        if (Q == 1)
            return a;
        else
            return (2 * a * Lucas(a, Q - 1, N) - Lucas(a, Q - 2, N)) % N;
    }

我也知道

(3) x_2n = 2 * (x_n)^2 - 1
(4) x_2n+1 = 2 * x_n+1 * x_n - a
(5) x_k(n+1) = 2 * x_k * x_kn - x_k(n-1)

(3) 和 (4) 应该有助于计算更大的 Q。但我不确定如何。不知何故，我认为 Q 的二进制形式。

感谢任何帮助。

【问题讨论】：

对于此类任务，您可以使用记忆功能来避免多次重新计算相同的子结果。此外，研究使用基于堆的堆栈来跟踪进度而不是递归 - 这可以防止 StackOverflowExceptions，还可以通过减少堆栈的流失来提高性能。
您正在计算这个 mod N 的事实可能很有用。 N 的值是否还有其他限制？
N 是一个正的非质数。我认为解决方案在这里的第三段中：programmingpraxis.com/2010/06/04/… 但我很难理解它并将其放入代码中。

标签： c# performance algorithm

【解决方案1】：

Here 可以看到如何使用矩阵乘以矩阵来找到第 N 个斐波那契数

      n
(1 1)
(1 0)

您可以利用这种方法来计算卢卡斯数，使用矩阵（对于您的情况x_n+l = 2 * a * x_n - x_n-l）

        n
(2a -1)
(1   0)

请注意，矩阵的 N 次方可以通过 exponentiation by squaring 的 log(N) 矩阵乘法找到

【讨论】：

【解决方案2】：

(3) x_2n = 2 * (x_n)^2 - 1
(4) x_2n+1 = 2 * x_n+1 * x_n - a

当您看到2n 时，您应该认为“这可能表示偶数”，同样2n+1 可能表示“这是一个奇数”。

您可以修改x 索引，使左侧有n（以便更容易理解这与递归函数调用的对应关系），请注意舍入。

3) 2n     n
=> n      n/2

4) it is easy to see that if x = 2n+1, then n = floor(x/2)
     and similarly n+1 = ceil(x/2)

所以，对于#3，我们有：（在伪代码中）

if Q is even
   return 2 * (the function call with Q/2) - 1

对于#4：

else // following from above if
   return 2 * (the function call with floor(Q/2))
            * (the function call with ceil(Q/2)) - a

然后我们还可以合并一点memoization来防止多次计算相同参数的返回值：

保留Q 值与返回值的映射。
在函数的开头，检查Q的值是否存在于map中。如果是，则返回相应的返回值。
返回时，将Q的值和返回值添加到map中。

【讨论】：

这实际上导致了一个 O(log n) 算法，这不是很明显

【解决方案3】：

第n个卢卡斯数的值为：

Exponentiation by squaring 可用于评估函数。例如，如果 n=1000000000，则 n = 1000 * 1000^2 = 10 * 10^2 * 1000^2 = 10 * 10^2 * (10 * 10^2 )^2。通过这种方式简化可以大大减少计算次数。

【讨论】：

【解决方案4】：

你可以得到一些改进（只是百万分之一......），而无需求助于真正花哨的数学。

首先让我们让数据流更明确一点：

    static BigInteger Lucas(BigInteger a, BigInteger Q, BigInteger N)
    {
        if (Q == 0)
        {
            return 1;
        }
        else if (Q == 1)
        {
            return a;
        }
        else
        {
            BigInteger q_1 = Lucas(a, Q - 1, N);
            BigInteger q_2 = Lucas(a, Q - 2, N);
            return (2 * a * q_1 - q_2) % N;
        }
    }

不出所料，这并没有真正改变性能。

但是，它确实清楚地表明我们只需要两个先前的值来计算下一个值。这让我们可以将函数倒置为迭代版本：

    static BigInteger IterativeLucas(BigInteger a, BigInteger Q, BigInteger N)
    {
        BigInteger[] acc = new BigInteger[2];
        Action<BigInteger> push = (el) => {
                acc[1] = acc[0];
                acc[0] = el;
        };
        for (BigInteger i = 0; i <= Q; i++)
        {
            if (i == 0)
            {
                push(1);
            }
            else if (i == 1)
            {
                push(a);
            }
            else
            {
                BigInteger q_1 = acc[0];
                BigInteger q_2 = acc[1];
                push((2 * a * q_1 - q_2) % N);
            }
        }
        return acc[0];
    }

可能有一种更清晰的方式来写这个，但它确实有效。它也快得多。它的速度要快得多，测量起来有点不切实际。在我的系统上，Lucas(4000000, 47, 4000000) 大约需要 30 分钟，IterativeLucas(4000000, 47, 4000000) 大约需要 2 毫秒。我想比较48，但我没有耐心。

使用模算术的这些属性，您可以挤出更多（可能是两倍？）：

(a + b) % n = (a%n + b%n) % n
(a * b) % n = ((a%n) * (b%n)) % n

如果您应用这些，您会发现a%N 出现了几次，因此您可以通过在循环之前预先计算一次来获胜。这在a 比N 大很多时特别有用；我不确定您的应用程序中是否会发生这种情况。

可能有一些聪明的数学技术可以将这个解决方案从水中淘汰，但我认为有趣的是，只需改组一点代码就可以实现这样的改进。

【讨论】：