计算n的快速方法! mod m 其中 m 是素数?

Posted

技术标签:

【中文标题】计算n的快速方法! mod m 其中 m 是素数?【英文标题】:Fast way to calculate n! mod m where m is prime? 【发布时间】:2012-04-01 10:43:34 【问题描述】:

我很好奇是否有一个好的方法来做到这一点。我当前的代码是这样的:

def factorialMod(n, modulus):
    ans=1
    for i in range(1,n+1):
        ans = ans * i % modulus    
    return ans % modulus

但似乎很慢!

我也无法计算n!然后应用素数模数,因为有时 n 太大以至于 n!明确计算是不可行的。

我也遇到了http://en.wikipedia.org/wiki/Stirling%27s_approximation,想知道这里是否可以以某种方式使用它?

或者,我如何在 C++ 中创建一个递归的、记忆化的函数?

【问题讨论】:

慢有多慢?从你的伪代码,我推断你是用 Python 计算的,对吗? 任何语言,真的;就语法而言,它在 C++ 中几乎相同。我在这里选择 Python 是因为它易于阅读。不过,即使在 C++ 中,我也需要一个更快的函数。 有一种非常快速的方法可以使用不变乘法或Montgomery reduction。这两种方法都消除了模数,并允许使用循环展开技术。 您可以将模数分解为主要因数,以更容易识别为零的情况,尽管这对较大的主要因数没有帮助 - 这有多大帮助取决于您对模数的了解,如果有的话,如果素因数分解让你喜欢。 如果 ans > 模数(信用:tech.groups.yahoo.com/group/primenumbers/messages/…) 【参考方案1】:

我在 quora 上找到了以下功能: 与 f(n,m) = n!模数;

function f(n,m:int64):int64;
         begin
              if n = 1 then f:= 1
              else f:= ((n mod m)*(f(n-1,m) mod m)) mod m;
         end;

可能会使用耗时的循环并乘以存储在字符串中的大量数字。此外,它适用于任何整数 m。 我找到此功能的链接:https://www.quora.com/How-do-you-calculate-n-mod-m-where-n-is-in-the-1000s-and-m-is-a-very-large-prime-number-eg-n-1000-m-10-9+7

【讨论】:

这与实现为递归函数的朴素算法完全相同。【参考方案2】:

如果素数 m 的 n = (m - 1) 则由 http://en.wikipedia.org/wiki/Wilson's_theorem n! mod m = (m - 1)

正如已经指出的那样 n!如果 n > m

,则 mod m = 0

【讨论】:

这没有帮助。 BlueRaja-Danny-Pflughoeft 已经提到过威尔逊定理,它并没有多大作用,因为你不能指望只需要 (m-1)! 或 (m-k)!对于小 k,他的回答涵盖了但你的回答没有。【参考方案3】:

如果我们要计算M = a*(a+1) * ... * (b-1) * b (mod p),我们可以使用下面的方法,假设我们可以快速加减乘乘(mod p),得到O( sqrt(b-a) * polylog(b-a) )的运行时间复杂度。

为简单起见,假设(b-a+1) = k^2 是一个正方形。现在,我们可以将我们的产品分成 k 个部分,即M = [a*..*(a+k-1)] *...* [(b-k+1)*..*b]。此产品中的每个因子的格式为p(x)=x*..*(x+k-1),对应的x

通过使用多项式的快速乘法算法,例如Schönhage–Strassen algorithm,以分治法的方式,可以找到多项式p(x) in O( k * polylog(k) )的系数。现在,显然有一种算法可以将k点替换为O( k * polylog(k) )中的相同k次多项式,这意味着我们可以快速计算p(a), p(a+k), ..., p(b-k+1)

C. Pomerance 和 R. Crandall 在“质数”一书中描述了这种将许多点代入一个多项式的算法。最终,当您拥有这些 k 值时,您可以将它们乘以 O(k) 并获得所需的值。

请注意,我们所有的操作都是在(mod p) 进行的。 确切的运行时间是O(sqrt(b-a) * log(b-a)^2 * log(log(b-a)))

【讨论】:

H. Cormen 等人在著名的“算法介绍”一书中(在 FFT 章节中)也描述了“将许多点代入一个多项式”的算法。【参考方案4】:

n! mod m 可以在 O(n1/2 + ε) 操作中计算,而不是简单的 O(n)。这需要使用 FFT 多项式乘法,并且仅适用于非常大的 n,例如n > 104.

算法的概要和一些时间可以在这里看到:http://fredrikj.net/blog/2012/03/factorials-mod-n-and-wilsons-theorem/

【讨论】:

这是一个比接受的答案更好的答案。【参考方案5】:

n 可以任意大

好吧,n 不能任意大 - 如果n >= m,那么n! ≡ 0 (mod m) (因为m 是因子之一,根据阶乘的定义).


假设n << m 并且您需要一个精确 值,据我所知,您的算法不会变得更快。但是,如果n > m/2,您可以使用以下身份(Wilson's theorem - 谢谢@Daniel Fischer!)

将乘法次数限制在大约m-n

(米-1)! ≡ -1 (mod m) 1 * 2 * 3 * ... * (n-1) * n * (n+1) * ... * (m-2) * (m-1) ≡ -1 (mod m) 嗯! * (n+1) * ... * (m-2) * (m-1) ≡ -1 (mod m) 嗯! ≡ -[(n+1) * ... * (m-2) * (m-1)]-1 (mod m)

这为我们提供了一种简单的方法来计算 n! (mod m) 中的 m-n-1 乘法,再加上 modular inverse:

def factorialMod(n, 模数): 答案=1 如果 n modinv(ans, modulus) ans = -1*ans + 模数 返回 ans % 模数

我们可以用另一种方式改写上述等式,这可能会或可能不会执行得稍微快一些。使用以下身份:

我们可以将方程改写为

嗯! ≡ -[(n+1) * ... * (m-2) * (m-1)]-1 (mod m) 嗯! ≡ -[(n+1-m) * ... * (m-2-m) * (m-1-m)]-1 (mod m) (术语的倒序) 嗯! ≡ -[(-1) * (-2) * ... * -(m-n-2) * -(m-n-1)]-1 (mod m) 嗯! ≡ -[(1) * (2) * ... * (mn-2) * (mn-1) * (-1)(mn-1)]-1 (mod m) 嗯! ≡ [(m-n-1)!]-1 * (-1)(m-n) (mod m)

这可以用 Python 写成如下:

def factorialMod(n, 模数): 答案=1 如果 n modinv(ans, modulus) #因为m是奇素数,(-1)^(m-n) = -1如果n是偶数,+1如果n是奇数 如果 n % 2 == 0: ans = -1*ans + 模数 返回 ans % 模数

如果您不需要 精确 值,生活会轻松一些 - 您可以使用 Stirling's approximation 计算 O(log n) 时间的近似值 (使用 exponentiation by squaring ).


最后,我应该提一下,如果这是时间紧迫的,并且您正在使用 Python,请尝试切换到 C++。根据个人经验,您应该期望速度会提高一个数量级或更多,仅仅是因为这正是原生编译代码 擅长 的那种受 CPU 限制的紧密循环>(此外,无论出于何种原因,GMP 似乎比 Python 的 Bignum 更精细)

【讨论】:

"这样,当m/2 < n < m时,你只需要计算(m/2)! * (-2)^(n-m/2-1) (mod m)"那么你可以做得更好。根据威尔逊定理,(m-1)! ≡ -1 (mod m) 如果m 是素数。现在(m-1)! = n! * (m - (m-n-1)) * ... * (m - 1) ≡ (-1)^(m-n-1) * n! * (m-n-1)! (mod m),所以n! ≡ (-1)^(m-n) * ((m-n-1)!)^(-1) (mod m)。所以你需要计算(m-n-1)! mod m,找到它的模逆(O(log m) 步),并在必要时调整符号。当n 接近m/2 时差别不大,但在n > 3m/4 左右时很好。 @DanielFischer:谢谢!我已将其包含在答案中。【参考方案6】:

扩展我的评论,对于 [100, 100007] 中的所有 n,这大约需要 50% 的时间,其中 m=(117 | 1117):

Function facmod(n As Integer, m As Integer) As Integer
    Dim f As Integer = 1
    For i As Integer = 2 To n
        f = f * i
        If f > m Then
            f = f Mod m
        End If
    Next
    Return f
End Function

【讨论】:

【参考方案7】:

将我的评论扩展到答案:

是的,有更有效的方法可以做到这一点。 但它们非常混乱。

所以除非你真的需要额外的性能,否则我不建议尝试实现这些。


关键是要注意模数(本质上是一个除法)将成为瓶颈操作。幸运的是,有一些非常快速的算法可以让您多次对相同的数字进行取模。

Division by Invariant Integers using Multiplication Montgomery Reduction

这些方法速度很快,因为它们基本上消除了模数。


仅这些方法就可以给您带来适度的加速。为了真正高效,您可能需要展开循环以实现更好的 IPC:

类似这样的:

ans0 = 1
ans1 = 1
for i in range(1,(n+1) / 2):
    ans0 = ans0 * (2*i + 0) % modulus    
    ans1 = ans1 * (2*i + 1) % modulus    

return ans0 * ans1 % modulus

但考虑到奇数的迭代次数并将其与我上面链接的方法之一结合起来。

有些人可能认为循环展开应该留给编译器。我会反驳说编译器目前还不够聪明,无法展开这个特定的循环。仔细看看你就会明白为什么了。


请注意,虽然我的回答与语言无关,但它主要用于 C 或 C++。

【讨论】:

如果刚刚对 3 个热门答案投了反对票的人的评论可能会很好。 如何在 C++ 中为阶乘模型 m 进行递归 + 记忆化? @JohnSmith TBH,Memoization 可能根本没有帮助 - 没有什么可记忆的。它可能变得有用的唯一方法是尝试使用素因子分解方法并使用windowing algorithm for exponentiation by squaring。 (窗口算法是一种记忆算法。)但是从1n 的所有整数的素数分解可能会比您当前的算法慢。 好吧,在我的例子中,我是从低 n 迭代到高 n,所以这是否意味着我可以通过存储已经计算的值来节省时间?对于较大的 n,似乎只进行几次迭代而不是从 i=1 到 n 或 n/2 会节省大量时间 嗯...没有什么可以“保存”的。知道几个迭代对你其余的迭代没有帮助。【参考方案8】:

假设您选择的平台的“mod”运算符足够快,那么您主要受限于计算 n! 的速度以及可用于计算它的空间。

那么它本质上是一个两步操作:

    计算 n! (有很多快速算法,这里不再赘述) 获取结果的 mod

没有必要使事情复杂化,尤其是在速度是关键因素的情况下。一般来说,在循环中尽可能少地执行操作。

如果您需要重复计算n! mod m,那么您可能需要记住执行计算的函数的值。一如既往,这是经典的空间/时间权衡,但查找表非常快。

最后,您可以将记忆与递归(如果需要,还可以使用蹦床)结合起来,让事情变得真正快。

【讨论】:

然而,对于大的n,计算n!然后执行mod是不可行的 不可行...为什么?由于内存限制?从问题来看,速度是问题,而不是内存。如果您希望内存占用尽可能小并然后优化速度,请更新您的问题以反映这一点。 -1 计算 n!然后mod很慢,请尝试计算2000000! mod 5250307 那样。 OP 在这个问题上做得更好,你应该交错乘法和取模。 @cdeszaq:您似乎缺少的是在计算机上将两个非常大的数字(大于寄存器的大小)相乘不是O(1):它更接近O(m log m) (m = #bits)。将两个 m 位数字相乘得到 (m+m) 位,因此您的方法大约需要 m log(m) + 2m log(m) + 3m log(m) + ... + nm log(m) = nm log(m)(n+1)/2 = O(mn^2 log(m)) 操作。但是,在每次操作后取模会产生大约 2(m log (m)) + 2(m log(m)) + ...n additions... + 2(m log(m)) = 2mn log(m) = O(mn log(m)),这明显更快,即使对于较小的 n 计算 n!very 大的 n 不仅速度慢,而且完全不可能,因为数字变得如此之大,您无法再解决它们。跨度>

以上是关于计算n的快速方法! mod m 其中 m 是素数?的主要内容,如果未能解决你的问题,请参考以下文章

hdu 5446 Unknown Treasure (Lucas定理+中国剩余定理+快速乘)

1113 矩阵快速幂

51Nod 1113 矩阵快速幂

51Nod - 1113 矩阵快速幂

HDU6440 Dream(费马小定理+构造) -2018CCPC网络赛1003

51nod 1113 矩阵快速幂 如题目