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

Question

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

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++中创建一个递归的,memoized函数？

Answer 1

n可以任意大

好吧,n不能随意大 - 如果n >= m,那么n! ? 0 (mod m) (因为m是因子的定义之一).

假设n << m你需要一个确切的值,据我所知,你的算法不能更快.但是,如果n > m/2,你可以使用以下身份(威尔逊定理 - 谢谢@Daniel Fischer!)

限制乘数的数量 m-n

(m-1)! ? -1 (mod m)
1 * 2 * 3 * ... * (n-1) * n * (n+1) * ... * (m-2) * (m-1) ? -1 (mod m)
n! * (n+1) * ... * (m-2) * (m-1) ? -1 (mod m)
n! ? -[(n+1) * ... * (m-2) * (m-1)]-1 (mod m)

这为我们提供了一种n! (mod m)在m-n-1乘法中计算的简单方法,以及模数逆:

def factorialMod(n, modulus):
    ans=1
    if n <= modulus//2:
        #calculate the factorial normally (right argument of range() is exclusive)
        for i in range(1,n+1):
            ans = (ans * i) % modulus   
    else:
        #Fancypants method for large n
        for i in range(n+1,modulus):
            ans = (ans * i) % modulus
        ans = modinv(ans, modulus)
        ans = -1*ans + modulus
    return ans % modulus

我们可以用另一种方式重新表述上述等式,这可能会或可能不会稍微快一些.使用以下标识:

我们可以将等式重新表述为

n! ? -[(n+1) * ... * (m-2) * (m-1)]-1 (mod m)
n! ? -[(n+1-m) * ... * (m-2-m) * (m-1-m)]-1 (mod m)
       (reverse order of terms)
n! ? -[(-1) * (-2) * ... * -(m-n-2) * -(m-n-1)]-1 (mod m)
n! ? -[(1) * (2) * ... * (m-n-2) * (m-n-1) * (-1)(m-n-1)]-1 (mod m)
n! ? [(m-n-1)!]-1 * (-1)(m-n) (mod m)

这可以用Python编写如下:

def factorialMod(n, modulus):
    ans=1
    if n <= modulus//2:
        #calculate the factorial normally (right argument of range() is exclusive)
        for i in range(1,n+1):
            ans = (ans * i) % modulus   
    else:
        #Fancypants method for large n
        for i in range(1,modulus-n):
            ans = (ans * i) % modulus
        ans = modinv(ans, modulus)

        #Since m is an odd-prime, (-1)^(m-n) = -1 if n is even, +1 if n is odd
        if n % 2 == 0:
            ans = -1*ans + modulus
    return ans % modulus

如果您不需要精确的值,生活会变得更容易 - 您可以使用斯特林的近似值来计算O(log n)时间的近似值(通过平方使用取幂).

最后,我应该提一下,如果这是时间关键的并且您正在使用Python,请尝试切换到C++.从个人的经验,你应该预计大约在速度或更多的订单数量级增加,仅仅是因为这正是那种CPU结合的紧密环的原生地编译代码擅长于(而且,无论出于何种原因,GMP似乎比Python的Bignum更精细调整.

Answer 2

将我的评论扩展到答案:

是的,有更有效的方法来做到这一点.但他们非常凌乱.

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

关键是要注意模数(基本上是一个除法)将成为瓶颈操作.幸运的是,有一些非常快速的算法允许您多次执行相同数量的模数.

• 使用乘法划分不变整数
• 蒙哥马利减少

这些方法很快,因为它们基本上消除了模量.

单独使用这些方法可以为您提供适度的加速.为了真正有效,您可能需要展开循环以允许更好的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++.

Answer 3

N!mod m可以用O(n ^{1/2 +ε})运算来计算,而不是天真的O(n).这需要使用FFT多项式乘法,并且仅适用于非常大的n,例如n> 10 ⁴.

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

Answer 4

如果我们想要计算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算法,以分而治之的方式,可以找到多项式的系数p(x) in O( k * polylog(k) ).现在,显然有一种算法用于替换k相同度数k的多项式中的点O( k * polylog(k) ),这意味着我们可以p(a), p(a+k), ..., p(b-k+1)快速计算.

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

请注意我们所有的操作都采取了(mod p).确切的运行时间是O(sqrt(b-a) * log(b-a)^2 * log(log(b-a))).