如何计算(n!)%1000000009

Question

我需要找到n!%1000000009.对于范围为1到20的k,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次.对于15左右的k值,我的代码运行速度太慢了.有没有办法更快？

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

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

Answer 1

由于您正在寻找nCr多个连续值，因此n可以使用以下内容：

    (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的倒数。您可以在此处找到如何计算它的说明，并点击链接以获取有关如何计算它的更多信息。bmExtended Euclidean Algorithm

2. 另一个可能更容易的选项是认识到，由于nCr * (n+1)/(n+1-r)必须给出一个整数，因此必须可以写出n+1-r == a*bwherea | 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. 现在你的除法保证是精确的，所以你只需计算它们，然后像​​以前一样继续进行乘法。 编辑根据评论，我不相信这种方法会起作用。

我可能会尝试的另一件事是在你的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;
        }
    }

这也可以节省一些宝贵的时间。