『Lucas定理以及拓展Lucas』

<更新提示>

<第一次更新>

<正文>

Lucas定理

在『组合数学基础』中，我们已经提出了$Lucas$定理，并给出了$Lucas$定理的证明，本文仅将简单回顾，并给出代码。

$Lucas$定理：当$p$为质数时，$C_n^m\equiv C_{n\ mod\ p}^{m\ mod\ p}*C_{n/p}^{m/p}(mod\ p)$。

在计算模域组合数时，如果模数较小，那么就可以尝试使用$Lucas$定理来递归求解，其时间复杂度为$O(plog_p\min(n,m))$。

$Code:$

inline long long Lucas(long long u,long long d)
{
    if (u<=Mod&&d<=Mod)return C(u,d) % Mod;
    else return Lucas(u/Mod,d/Mod) * C(u%Mod,d%Mod) % Mod;
}

本文给出了一种递归写法，由于取模操作的存在，在组合数计算函数$C$内也需要进行一些特判，避免一些无意义式子的计算，影响答案。

$Code:$

inline long long C(long long u,long long d)
{
    if ( u>d || u<0 || d<0 )return 0LL;
    //...
}

关于组合数的计算，我们也有多种方法，常用的几种如下：

$1.$ 当$n$，$m$都比较小的时候，可以根据组合数的阶乘计算公式来预处理$n$，$m$范围以内的阶乘以及阶乘的逆元，然后根据询问$O(1)$地回答组合数。

$2.$ 当$n$的范围较大，$m$的范围较小时，由于不适合预处理，所以我们可以直接利用有关$m$组合数计算式$C_n^m=\frac{n*(n-1)*...*(n-m+1)}{m*(m-1)*...*1}$来模拟计算，同理，分母利用逆元操作即可，时间复杂度$O(m)$。

Exlucas算法

与之对应的，$Exlucas$算法是一种用来解决与$Lucas$定理形式很像但更具有一般性的问题的算法。为什么不叫$Exlucas$定理而叫$Exlucas$算法呢，是因为这个算法和$Lucas$定理没什么关系，只是一个数论算法罢了。

对于求解模域组合数$C_n^m\%p$，当$p$不一定是质数时，可以使用$Exlucas$算法求解。

对于这样一个问题，我们考虑对模数进行分解，设$p=\prod_{i=1}^kp_i^{a_i}$，答案$x=C_n^m\%p$，则可以得到$$x\equiv C_n^ m(mod\ p_i^{a_i})$$
这是一个线性同余方程组，可以使用中国剩余定理求解。

那么现在我们的问题就转换成了求解$C_n^m\%p_i^{a_i}$。将组合数写为阶乘的计算式，即$C_n^m=\frac{n!}{m!(n-m)!}$。

因为$\frac{n!}{m!(n-m)!}$与$p_i^{a_i}$不可避免的会存在公约数，所以不能直接使用$Exeuclid$等算法来求解逆元，那么我们就要考虑把$\frac{n!}{m!(n-m)!}$中有关$p_i$的项都提取掉，利用$\frac{n!}{m!(n-m)!}$是整数这一特点将分子该项的指数直接减掉分母该项的指数，再求解其余部分以及逆元，最后就能避免不互质的情况从而求解答案。

先考虑一个简单的问题，如何计算一个形如$n!$的数中含有多少个因子$p$，则和阶乘分解一题类似，若设$f(n)$表示$n!$中含有因子$p$的个数，则有：

\[f(n)=\begin{cases}0\ (n<p)\\\lfloor \frac{n}{p} \rfloor +f(\lfloor \frac{n}{p} \rfloor)\ (n\geq p)\end{cases} \]

简单地利用$for$循环求解即可。

求得阶乘中因子$p_i$的数量后，我们的问题就变成为如何求解$n!$中除去所有$p_i$后取模$p_i^{a_i}$的值。对于每一项，这又需要我们分两种情况考虑：

$1.$ 该项是除去$p_i$后得到的，我们将所有这样的项放在一起发现，发现除去$p_i$后，剩下的系数乘积也是一个阶乘，使用递归求解。
$2.$ 该项是阶乘中原本的一项，除去情况$1.$，我们发现阶乘中若干个连续的项在模$p_i^{a_i}$意义下恰好构成了循环节，其长度不超过$p_i^{a_i}$，先暴力计算一个循环节，然后快速幂即可。对于不包括在完整循环节中的项，由于其长度小于循环节，可以暴力计算。

然后就能将求得的答案通过计算逆元求解，最后，乘上提取出的$p_i$剩余的若干次方即可。

考虑一个简单的例子就能加深理解：

求解$n!\ mod\ p_i^{p_i}$，此时，$n=19$，$p_i=3$，$a_i=2$。

\[n!=1*2*3*...*19\\=(1*2*4*5*7*8*10*11*13*14*16*17*19)*(3^6*6!) \\=(1*2*4*5*7*8)^2*19*(3^6*6!) \]

可以证明，其时间复杂度与$O(plog_2p)$同级。

$Code:$

#include <bits/stdc++.h>
using namespace std;
const int SIZE=300;
long long A,B,Mod,m[SIZE],r[SIZE];
inline void input(void)
{
    scanf("%lld%lld%lld",&A,&B,&Mod);
    //计算C(A,B)%Mod
}
inline long long power(long long a,long long b,long long p)
{
    long long res = 1;
    while (b)
    {
        if (1&b)res = res * a % p;
        b >>= 1;
        a = a * a % p;
    }
    return res;
}
inline long long Exeuclid(long long a,long long &x,long long b,long long &y,long long c)
{
    if (!b){x=c/a,y=0;return a;}
    else
    {
        long long p = Exeuclid(b,x,a%b,y,c);
        long long x_ = x , y_ = y;
        x = y_ , y = x_ - a/b * y_;
        return p;
    }
}
inline long long inv(long long a,long long p)
{
    long long x_,y_;
    Exeuclid(a,x_,p,y_,1);
    return (x_+p)%p;
}
inline long long calc(long long x,long long val,long long p)
{
    //计算x!中除去val后模p意义下的值
    if (!x)return 1;
    long long res = 1 , last = x%p;
    for (long long i=1;i<=p;i++)
        if (i%val)res = res * i % p;
    res = power(res,x/p,p);
    for (long long i=1;i<=last;i++)
        if (i%val)res = res * i % p;
    return res * calc(x/val,val,p) % p;
}
inline long long C(long long d,long long u,long long val,long long Pow)
{
    long long mulup=calc(d,val,Pow) , k=0;
    long long muldown1=calc(u,val,Pow) , muldown2=calc(d-u,val,Pow);
    for (long long i=d;i;i/=val)k+=i/val;
    for (long long i=u;i;i/=val)k-=i/val;
    for (long long i=d-u;i;i/=val)k-=i/val;
    //计算剩余的val的指数
    return mulup * inv(muldown1,Pow) % Pow * inv(muldown2,Pow) % Pow * power(val,k,Pow) % Pow;
}
inline long long CRT(int cnt)
{
    long long m_ = Mod , M[SIZE] = {} , t[SIZE] = {} , res = 0;
    for (int i=1;i<=cnt;i++)
        M[i] = m_ / m[i];
    for (int i=1;i<=cnt;i++)
    {
        long long y;
        Exeuclid(M[i],t[i],m[i],y,1);
        res = (res + r[i] % m_ * M[i] % m_ * t[i] % m_) % m_;
    }
    return ( res % m_ + m_ ) % m_;
}
inline long long Exlucas(void)
{
    int temp = Mod , cnt = 0 ;
    for (long long i=2;i*i<=Mod;i++)
    {
        if (temp%i==0)
        {
            m[++cnt] = 1;
            while (temp%i==0)
                temp /= i , m[cnt] *= i;
            r[cnt] = C(A,B,i,m[cnt]);
        }
    }
    if (temp>1) m[++cnt] = temp , r[cnt] = C(A,B,temp,m[cnt]);
    return CRT(cnt);
}
int main(void)
{
    input();
    printf("%lld\n",Exlucas());
    return 0;
}

<后记>

posted @ 2019-04-24 20:34 Parsnip 阅读(470) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

『Lucas定理以及拓展Lucas』

Lucas定理

Exlucas算法

公告