中国剩余定理与扩展 Lucas定理与扩展 学习笔记

中国剩余定理

问题

求同余方程组

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪x≡c1(modm1)x≡c2(modm2)x≡c3(modm3)...x≡ck(modmk)

其中满足(mi,mj)=1,1<=i≠j<=k
x的最小正（非负）整数解

结论

令M=m1∗m2∗m3∗...∗mk
则x=(∑i=1kci∗Mmi∗inv(Mmi,mi))%M

证明

a.在模M意义下x只有唯一解 （有多解那还了得）
b.令ni满足ni%Mmi=0且ni%mi=1，则N=∑i=1kci∗ni为原问题的一个解
c.根据上面的式子容易得出ni=Mmix,ni=miy+1，则Mmix=miy+1，即Mmix−miy=1
d.由于mi两两互质，所以Mmi与mi也互质，令a=Mmi,b=mi，则ax−by=1
e.可以发现我们已经将其化简成扩展欧几里得的基本形式ax+by=(a,b)，其中也可以认为x为a在模b意义下的逆元

代码

codevs 3990

#include<algorithm>#include<iostream>#include<cstring>#include<cstdio>#include<cmath>using namespace std;#define LL long longLL k,l,r,n,M,x,y,Min,ans,m[15],c[15];void exgcd(LL a,LL b,LL &x,LL &y){    if (!b) x=1,y=0;    else exgcd(b,a%b,y,x),y-=a/b*x;}int main(){    scanf("%d%lld%lld",&k,&l,&r);    M=1LL;    for (int i=1;i<=k;++i)    {        scanf("%lld%lld",&m[i],&c[i]);        M*=m[i];    }    for (int i=1;i<=k;++i)    {        LL a=M/m[i],b=m[i];        exgcd(a,b,x,y);        x=(x%b+b)%b;        if (!x) x+=b;        n+=c[i]*a*x;    }    n%=M;    if (!n) n+=M;    if (r>=n)        ans=(r-n)/M+1;    if (l>=n) ans=ans-((l-n)/M+1);    if ((l-n)%M==0) ++ans;    if (ans)    {        if (l<=n) Min=n;        else Min=n+((l-n)/M+1)*M;    }    printf("%lld\n%lld\n",ans,Min);}

扩展中国剩余定理

问题

求同余方程组

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪x≡c1(modm1)x≡c2(modm2)x≡c3(modm3)...x≡ck(modmk)

x的最小正（非负）整数解

结论

对于两个方程
x≡c1(modm1)
x≡c2(modm2)
将其合并成一个方程，有解条件为(m1,m2)|(c2−c1)
m=m1m2(m1,m2)
c=(inv(m1(m1,m2),m2(m1,m2))∗(c2−c1)(m1,m2))%m2(m1,m2)∗m1+c1
最终得出一个式子x≡c(modm1)
x=c%m即为原问题的一个解

证明

a.将两个方程写成x=m1k1+c1,x=m2k2+c2，两式联立得m1k1+c1=m2k2+c2，移项m1k1=c2−c1+m2k2，
b.根据贝祖定理，以上等式有解充要条件为(m1,m2)|(c2−c1)
c.将等式两边同除(m1,m2)得m1(m1,m2)∗k1=(c2−c1)(m1,m2)+m2(m1,m2)∗k2，即m1(m1,m2)k1≡(c2−c1)(m1,m2)(modm2(m1,m2))，进一步化简k1≡inv(m1(m1,m2),m2(m1,m2))∗(c2−c1)(m1,m2)(modm2(m1,m2))，k1=inv(m1(m1,m2),m2(m1,m2))∗(c2−c1)(m1,m2)+ym2(m1,m2)
d.将其回代x=m1k1+c1，得x=inv(m1(m1,m2),m2(m1,m2))∗(c2−c1)(m1,m2)∗m1+ym1m2(m1,m2)+c1，即x≡inv(m1(m1,m2),m2(m1,m2))∗(c2−c1)(m1,m2)∗m1+c1(modm1m2(m1,m2))
e.至此我们又将其化简成了x≡c(modm)的形式，其中m=m1m2(m1,m2)，c=(inv(m1(m1,m2),m2(m1,m2))∗(c2−c1)(m1,m2))%m2(m1,m2)∗m1+c1

代码

pku 2891

#include<algorithm>#include<iostream>#include<cstring>#include<cstdio>#include<cmath>using namespace std;#define LL long long#define N 1005int k;LL c[N],m[N],c1,c2,m1,m2,t;bool flag;LL gcd(LL a,LL b){    if (!b) return a;    else return gcd(b,a%b);}void exgcd(LL a,LL b,LL &x,LL &y){    if (!b) x=1LL,y=0LL;    else exgcd(b,a%b,y,x),y-=a/b*x;}LL inv(LL a,LL b){    LL x=0LL,y=0LL;    exgcd(a,b,x,y);    x=(x%b+b)%b;    if (!x) x+=b;    return x;}int main(){    while (~scanf("%d",&k))    {        flag=true;        for (int i=1;i<=k;++i)            scanf("%I64d%I64d",&m[i],&c[i]);        for (int i=2;i<=k;++i)        {            m1=m[i-1],m2=m[i],c1=c[i-1],c2=c[i];            t=gcd(m1,m2);            if ((c2-c1)%t!=0) {flag=false;break;}            m[i]=m1*m2/t;            c[i]=inv(m1/t,m2/t)*((c2-c1)/t)%(m2/t)*m1+c1;            c[i]=(c[i]%m[i]+m[i])%m[i];        }        if (!flag) puts("-1");        else printf("%I64d\n",c[k]);    }}

Lucas定理

问题

求Cmn%p，其中p为质数

结论

令
n=nk∗pk+nk−1∗pk−1+...+n2∗p2+n1∗p+n0
m=mk∗pk+mk−1∗pk−1+...+m2∗p2+m1∗p+m0
则Cmn=∏i=0kCmini

证明

去问卢卡斯
听说要用二项式定理什么的
我怎么可能会这种东西
记住就行了

代码

zoj 3557

#include<algorithm>#include<iostream>#include<cstring>#include<cstdio>#include<cmath>using namespace std;#define LL long longLL n,m,Mod;LL fast_pow(LL a,LL p){    LL ans=1LL;    for (;p;p>>=1,a=a*a%Mod)        if (p&1)            ans=ans*a%Mod;    return ans;}LL inv(LL x){    return fast_pow(x,Mod-2);}LL C(LL n,LL m){    if (m>n) return 0LL;    LL up=1LL,down=1LL;    for (LL i=n-m+1;i<=n;++i) up=up*i%Mod;    for (LL i=1;i<=m;++i) down=down*i%Mod;    return up*inv(down)%Mod;}LL lucas(LL n,LL m){    if (m>n) return 0LL;    LL ans=1;    for (;m;n/=Mod,m/=Mod)        ans=ans*C(n%Mod,m%Mod)%Mod;    return ans;}int main(){    while (~scanf("%lld%lld%lld",&n,&m,&Mod))        printf("%lld\n",lucas(n-m+1,m));}

扩展Lucas定理

问题

求Cmn%p

结论

令
p=p1k1∗p2k2∗...∗pqkq
列出同余方程组

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪ans≡c1(modp1k1)ans≡c2(modp2k2)...ans≡cq(modpqkq)

其中c1...cq是对于每一个Cmn%piki求出的答案
然后根据中国剩余定理合并
可见piki并不是质数，而是一个质数的幂的形式
对于如何求Cmn%piki将在证明中给出

证明

a.由于同于方程组在模p=p1k1∗p2k2∗...∗pqkq意义下有唯一解，可以证明上面做法的正确性
b.由于p=p1k1∗p2k2∗...∗pqkq实质上是将p进行质因数分解，所以piki满足两两互质，可以直接用中国剩余定理合并
c.对于如何求Cmn%piki
根据Cmn=n!m!(n−m)!，只要分别求出n!%piki,m!%piki,(n−m)!%piki就可以通过逆元求出Cmn%piki
d.对于如何求n!%piki
我们以n=19,pi=3,ki=2为例
n!=1∗2∗3∗4∗5∗6∗7∗8∗9∗10∗11∗12∗13∗14∗15∗16∗17∗18∗19
=(1∗2∗4∗5∗7∗8∗10∗11∗13∗14∗16∗17∗19)∗36∗(1∗2∗3∗4∗5∗6)
根据这个例子发现，求解n!可以分为3部分：第一部分是pi的幂的部分，也就是36即pi⌊npi⌋，可以直接求解；第二部分是一个新的阶乘，也就是6!即⌊npi⌋!，可以递归下去求解；第三部分是除前两部分之外剩下的数
考虑第三部分如何求解
发现第三部分在模piki意义下是以piki为周期的，即(1∗2∗4∗5∗7∗8)≡(10∗11∗13∗14∗16∗17)(modpiki)，所以只求piki长度的即可；但是还剩下一个孤立的19，可以发现剩下孤立的数长度不会超过piki，只需要暴力求解即可
e.最后一个问题是对于求出的m!%piki和(n−m)!%piki有可能与piki不互质，无法求逆元
所以要将m!%piki和(n−m)!%piki中质因子pi先全部除去，求出逆元后再全部乘回去
计算n!中质因子p的个数x的公式为x=⌊np⌋+⌊np2⌋+⌊np3⌋+...
递推式也可以写为f(n)=f(⌊np⌋)+⌊np⌋

代码

codeforces2015ICL,Finals,Div.1#J

#include<algorithm>#include<iostream>#include<cstring>#include<cstdio>#include<cmath>using namespace std;#define LL long longLL n,m,MOD,ans;LL fast_pow(LL a,LL p,LL Mod){    LL ans=1LL;    for (;p;p>>=1,a=a*a%Mod)        if (p&1)            ans=ans*a%Mod;    return ans;}void exgcd(LL a,LL b,LL &x,LL &y){    if (!b) x=1LL,y=0LL;    else exgcd(b,a%b,y,x),y-=a/b*x;}LL inv(LL A,LL Mod){    if (!A) return 0LL;    LL a=A,b=Mod,x=0LL,y=0LL;    exgcd(a,b,x,y);    x=((x%b)+b)%b;    if (!x) x+=b;    return x;}LL Mul(LL n,LL pi,LL pk){    if (!n) return 1LL;    LL ans=1LL;    if (n/pk)    {        for (LL i=2;i<=pk;++i)            if (i%pi) ans=ans*i%pk;        ans=fast_pow(ans,n/pk,pk);    }    for (LL i=2;i<=n%pk;++i)        if (i%pi) ans=ans*i%pk;    return ans*Mul(n/pi,pi,pk)%pk;}LL C(LL n,LL m,LL Mod,LL pi,LL pk){    if (m>n) return 0LL;    LL a=Mul(n,pi,pk),b=Mul(m,pi,pk),c=Mul(n-m,pi,pk);    LL k=0LL,ans;    for (LL i=n;i;i/=pi) k+=i/pi;    for (LL i=m;i;i/=pi) k-=i/pi;    for (LL i=n-m;i;i/=pi) k-=i/pi;    ans=a*inv(b,pk)%pk*inv(c,pk)%pk*fast_pow(pi,k,pk)%pk;    return ans*(Mod/pk)%Mod*inv(Mod/pk,pk)%Mod;}int main(){    scanf("%I64d%I64d%I64d",&n,&m,&MOD);    for (LL x=MOD,i=2;i<=MOD;++i)        if (x%i==0)        {            LL pk=1LL;            while (x%i==0) pk*=i,x/=i;            ans=(ans+C(n,m,MOD,i,pk))%MOD;        }    printf("%I64d\n",ans);}