lucas定理与扩展lucas定理(HDU 3037 以及 2015 ICL, Finals, Div. 1 J.Ceizenpok’s formula)

最新推荐文章于 2022-04-08 11:28:04 发布

hash[072]

最新推荐文章于 2022-04-08 11:28:04 发布

阅读量337

点赞数

分类专栏： M

本文链接：https://blog.csdn.net/weixin_43958964/article/details/105884869

版权

M 专栏收录该内容

32 篇文章 0 订阅

订阅专栏

之前时不时的学一点，属实浮躁，总结一下。
Lucas定理是用来求 c(n,m) mod p，p为素数的值。
表达式C(n,m)%p=C(n/p,m/p)*C(n%p,m%p)%p 。
时间复杂度： $log_p n * p$
推导：采用的是构造证明，有错误还请指出。
在这里插入图片描述举个栗子：
题目：HDU 3037
题意：n个不同的树中保存不超过m个物品（它们是相同的）总共有多少种方法。
比方说2棵树一个物品，有三种方法，不放，一个放树1，或者一个放树2。
由于n,m的值非常大，不可能n^2递推，由于m个物品是相同的，我们从m个为m个物品入手，将他们分堆，本质是x1+x2+…+xn=m，可以先分 $C_{m-1}^i$ 即i+1份正整数解，需要放到n个不同的树中，即 $C_n^{i+1}$ 种情况，之后作累加，但有可能只存放一个物品，或者两个…，所以总式子是 $\sum_{j=1}^m\sum_{i=1}^jC_{j-1}^{i-1}C_{n}^i+1$ 。如果计算复杂度将上升到 $O(n^2)$ ，显然不行。考虑本来就是求x1+x2+…+xn=m解得个数，这个式子是非负整数解，那么构造设yi=xi+1，y1+y2+…+yn=n+m，即求他的正整数解，即 $C_{n+m-1}^{n-1}$ ,n-1就是分为n份，最终也就是 $\sum_{i=0}^mC_{n+i-1}^{n-1}$ ,又因为C(n,m)=C(n-1,m-1)+C(n-1,m)[定义推导即可证明]， $\sum_{i=0}^mC_{n+i-1}^{n-1}=C_{n+m}^m$ ，又由于n,m得范围很大，p是质数，借用上述卢卡斯定理递归解决，中间带入了费马小定理(由欧拉定理推导而来(欧拉定理由最小正简化剩余来证明))。

#include<cmath>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<cstdlib>
#include<istream>
#include<vector>
#include<stack>
#include<set>
#include<map>
#include<algorithm>
#include<queue>
#define inf 0x3f3f3f3f
#define llinf 0x3f3f3f3f3f3f3f3f
#define MAX_len 50100*4
using namespace std;
typedef long long ll;
const int mod=998244353;
ll quickpow(ll a,ll n,ll p)
{
    ll res=1;
    while(n)
    {
        if(n&1)
        {
            res=(res*a)%p;
        }
        n>>=1;
        a=(a*a)%p;
    }
    return res;
}
ll comb(ll a,ll b,ll p)
{
    if(a<b)
        return 0;
    if(a==b)
        return 1;
    ll i,j,k,sum1=1,sum2=1;
    if(2*b>a)
        b=a-b;
    for(i=0;i<b;i++)
    {
        sum1=(sum1*(a-i))%p;
        sum2=(sum2*(b-i))%p;
    }
    ll temp=(sum1*quickpow(sum2,p-2,p))%p;
    return temp;
}
ll lucas(ll n,ll m,ll p)
{
    ll ans=1;
    if(m==0)
        return 1;
    while(n&&m)
    {
        ans=(ans*comb(n%p,m%p,p))%p;
        n/=p;
        m/=p;
    }
    return ans;
}
int main()
{
    int T;
    scanf("%d",&T);
    while(T--)
    {
        ll n,m,p;
        scanf("%lld %lld %lld",&n,&m,&p);
        printf("%lld\n",lucas(n+m,m,p));
    }
    return 0;
}

扩展Lucas定理：当p不是素数的时候怎么办，如上的证明过程由于 $C_p^i(mod p)=0$ 由于p是质数的关系利用费马小定理必是0，得到 $1+x^p$ ，但p不是质数时取模过程不一定是0，中途求解此式子需要利用扩展欧几里德求逆元，还需要互质的前提条件。思路及那个p分解为 $\prod_ip_i^{k_i}$ ，pi都为质数，做唯一分解，用中国剩余定理合并。
现在需要取解决 $C_n^m(mod$ $a_i)$ , $ai=p_i^{k_i}$ 。我们先对n!做操作，举个栗子：n=19,pi=3,ki=2。
$n!=(3*6*9*12*15*18)*(1*2*4*5*7*8)*(10*11*13*14*17)*19=(3^6*6!)*(1*2*4*5*7*8)^2*19$ ，分析可以看出， $\lfloor \frac{n}{p_i} \rfloor$ !* $p_i^{\lfloor \frac{n}{p_i} \rfloor}$ ，之后是以 $p_i$ 为一个周期作数，每个周期由于是mod，1和1+ $p_i^{k_i}$ mod $p_i^{k_i}$ 同余。最后直接处理一下单个的19，也就是未成循环节的部分。处理完(1 * 2 * 4 * 5 * 7 * 8)，之后递归处理6!，跟之前一样的处理方式，最后只需要再判断一下总的是p的多少次幂，到这里就解决了。
题目：Gym - 100633J
裸题 $C_n^k mod m$

#include<cmath>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<cstdlib>
#include<istream>
#include<vector>
#include<stack>
#include<set>
#include<map>
#include<algorithm>
#include<queue>
#define inf 0x3f3f3f3f
#define llinf 0x3f3f3f3f3f3f3f3f
#define MAX_len 50100*4
using namespace std;
typedef long long ll;
//const int mod=998244353;
ll tot=0;
ll prime[1000100],mprime[1000100];
void init(ll n)
{
    ll temp=n;
    for(ll i=2;i*i<=n;i++)
    {
        if(temp%i==0)
        {
            prime[tot]=i;
            mprime[tot]=1;
            while(temp%i==0)
            {
                temp/=i;
                mprime[tot]*=i;
            }
            tot++;
        }
    }
    if(temp>1)
    {
        prime[tot]=temp;
        mprime[tot]=temp;
        tot++;
    }
}
ll quickpow(ll a,ll n,ll p)
{
    ll res=1;
    while(n)
    {
        if(n&1)
        {
            res=(res*a)%p;
        }
        n>>=1;
        a=(a*a)%p;
    }
    return res;
}
ll C1(ll n,ll pi,ll pk)
{
    if(n==0)
        return 1;
    ll ans=1;
    for(ll i=2;i<=pk;i++)
    {
        if(i%pi==0)
            continue;
        ans=(ans*i)%pk;
    }
    ans=quickpow(ans,n/pk,pk);
    ll temp=n%pk;
    for(ll i=2;i<=temp;i++)
    {
        if(i%pi==0)
            continue;
        ans=(ans*i)%pk;
    }
    return ans*C1(n/pi,pi,pk)%pk;
}
ll ex_gcd(ll a,ll b,ll &x,ll &y)
{
    if(b==0)
    {
        x=1,y=0;
        return a;
    }
    ll t=ex_gcd(b,a%b,x,y);
    ll temp=y;
    y=x-(a/b)*y;
    x=temp;
    return t;
}
ll revers(ll a,ll b)
{
    ll x,y;
    ll g=ex_gcd(a,b,x,y);
    //逆元一定有解
    ll temp=(x%b+b)%b;
        return temp;
}
ll comb(ll n,ll m,ll cur,ll mod)
{
    ll pi=prime[cur],pk=mprime[cur];
    ll a,b,c;
    a=C1(n,pi,pk);
    b=C1(m,pi,pk);
    c=C1(n-m,pi,pk);
    ll am=0;
    ll temp=n;
    while(temp)
    {
        am+=temp/pi;
        temp/=pi;
    }
    temp=m;
    while(temp)
    {
        am-=temp/pi;
        temp/=pi;
    }
    temp=n-m;
    while(temp)
    {
        am-=temp/pi;
        temp/=pi;
    }
    ll ans=a*((revers(b,pk)%pk)*(revers(c,pk)%pk)%pk)*(quickpow(pi,am,pk));
    return (ans%mod)*((mod/pk)%mod)*(revers(mod/pk,pk)%mod);
}
ll exlucas(ll n,ll m,ll mod)
{
    ll ans=0;
    for(ll i=0;i<tot;i++)
    {
        ans=(ans+comb(n,m,i,mod))%mod;
    }
    return ans;
}
int main()
{
    ll n,m,mod;
    scanf("%I64d %I64d %I64d",&n,&m,&mod);
    init(mod);
    printf("%I64d\n",exlucas(n,m,mod));
    return 0;
}