如何计算(n!)%1000000009

How to calculate (n!)%1000000009

本文关键字：%1000000009 计算 何计算      更新时间：2023-10-16

我需要找到n!%1000000009。在1到20的范围内，N的类型是2^k。我使用的函数是:

#define llu unsigned long long
#define MOD 1000000009
llu mulmod(llu a,llu b) // This function calculates (a*b)%MOD caring about overflows
{
    llu x=0,y=a%MOD;
    while(b > 0)
    {
        if(b%2 == 1)
        {
            x = (x+y)%MOD;
        }
        y = (y*2)%MOD;
        b /= 2;
    }
    return (x%MOD);
}
llu fun(int n)   // This function returns answer to my query ie. n!%MOD
{
    llu ans=1;
    for(int j=1; j<=n; j++)
    {
        ans=mulmod(ans,j);
    }
    return ans;
}

我的需求是这样的，我需要调用函数'fun' n/2次。我的代码运行太慢，k值在15左右。有更快的方法吗?

编辑:我在实际计算2 *((张)C (2 ^ (k - 1) 1)) * (((2 ^ (k - 1)) !) ^ 2]我的范围2 ^ 2 ^ k (k - 1)。我的程序要求(nCr)%MOD关心溢出。

编辑:我需要一个有效的方法来找到nCr%MOD的大n。

mulmod例程可以大大提高k的速度。

1)'%'是多余的，因为(a + b)都小于n。
-评估c = a+b; if (c>=N) c-=N;
就足够了2)可同时处理多个比特;参见"俄罗斯农民算法"的优化
3) a * b实际上足够小，适合64位unsigned long long，没有溢出

由于实际问题是关于nCr mod M，因此高级优化需要使用递归式

(n+1)Cr mod M = (n+1)nCr / (n+1-r) mod M.

由于公式((nCr) mod M)*(n+1)不能被(n+1-r)整除，因此除法需要通过与模逆(n+r-1)^(-1)的乘法来实现。模逆b^(-1)是b^(M-1)，因为M是素数。(否则就 b ^(φ ()),在φ是欧拉Totient函数。)

模幂运算通常是通过重复平方来实现的，在这种情况下，每个除数需要进行约45次模乘法。

如果你可以使用递归式

nC(r+1) mod M = nCr * (n-r) / (r+1) mod M

只需要计算(r+1)^(M-1) mod M一次

既然您正在为n的多个顺序值寻找nCr，您可以使用以下命令:

    (n+1)Cr = (n+1)! / ((r!)*(n+1-r)!)
    (n+1)Cr = n!*(n+1) / ((r!)*(n-r)!*(n+1-r))
    (n+1)Cr = n! / ((r!)*(n-r)!) * (n+1)/(n+1-r)
    (n+1)Cr = nCr * (n+1)/(n+1-r)

这样可以避免为每个i显式调用阶乘函数。

此外，要保存对nCr的第一次调用，可以使用:

    nC(n-1) = n //where n in your case is 2^(k-1).

编辑:
正如Aki Suihkonen所指出的，(a/b) % m != a%m / b%m。所以上面的方法。上面的方法不能一开箱即用。有两种不同的解决方案:

  1. 1000000009是素数，这意味着a/b % m == a*c % m，其中c是b的倒数模m。你可以在这里找到如何计算它的解释，并点击Extended Euclidean Algorithm的链接了解更多关于如何计算它的信息。

  2. 另一个可能更容易的选择是认识到，由于nCr * (n+1)/(n+1-r)必须给出一个整数，因此必须可以将n+1-r == a*b写入a | nCr和b | n+1(这里的|意味着除法，如果您喜欢，您可以将其重写为nCr % a == 0)。在不失一般性的前提下，先设a = gcd(n+1-r,nCr)，再设b = (n+1-r) / a。这给出了(n+1)Cr == (nCr / a) * ((n+1) / b) % MOD。现在你的除法保证是准确的，所以你只需要计算它们，然后像以前一样进行乘法运算。EDIT根据评论，我不相信这种方法会起作用

另一个我可能会尝试的是在你的llu mulmod(llu a,llu b) 

    llu mulmod(llu a,llu b)
    {
        llu q = a * b;
        if(q < a || q < b) // Overflow!
        {
            llu x=0,y=a%MOD;
            while(b > 0)
            {
                if(b%2 == 1)
                {
                    x = (x+y)%MOD;
                }
                y = (y*2)%MOD;
                b /= 2;
            }
            return (x%MOD);
        }
        else
        {
            return q % MOD;
        }
    }
相关文章：