Lucas定理&扩展Lucas定理&BSGS算法&扩展BSGS算法

（1）Lucas定理：
若p为素数：

$$C_n^m≡∏_{i=0}^kC_{ni}^{mi}(mod\,p)$$
其中 $$n=n_kp^k+n_{k−1}p^{k−1}+...+n_0$$
$$m=m_kp^k+m_{k−1}p^{k−1}+...+m_0$$
代码实现可以简便的理解为：
$$C_m^n=C_{m/p}^{n/p}×C_{m\,mod\,p}^{n\,mod\,p}(mod\,p)。$$

LL queryc(int x,int y) {
    if(y<x) return 0;
    if(x>=mod || y>=mod) 
        return 1LL*queryc(x/mod,y/mod)%mod*queryc(x%mod, y%mod)%mod;
    int ans1=1,ans2=1;  for(int i=y;i>x;i--) ans1=1LL*ans1*i%mod;
    for(int i=1;i<=y-x;i++) ans2=1LL*ans2*i%mod;
    return 1LL*ans1*Pow(ans2, mod-2)%mod;
}

(2)扩展Lucas：

若p不是素数，我们将p分解质因数，将$C_m^n$分别按照（1）中的方法求对p的质因数的模，然后用中国剩余定理合并。
例如：
当我们需要计算$C_m^n\,mod\,p$，其中$p=p_1^{q_1}×p_2^{q_2}×...×p_k^{q_k}$，我们可以求出：

$C_m^n≡a_i(mod\,p_i^{qi})(1<i<k)$

然后对于方程组：
$x≡ai(mod\,p_i^{q_i})(1<i<k)$

我们可以求出满足条件的最小的x,记为$x_0$。
那么我们有：
$C_m^n≡x_0(mod\,p)$

但是，我们发现，$pi^{qi}$并不是一个素数，它是某个素数的某次方。
下面我们介绍如何计算$C_m^n\,mod\,p^t$(t≥2,p为素数)。
计算$C_m^n\,mod\,pt$：
我们知道，$C_m^n=\frac{n!}{m!(n−m)!}$，若我们可以计算出n! mod $p^t$，我们就能计算出m! mod pt以及(n−m)! mod $p^t$。我们不妨设x=n! mod $p^t$,y=m! mod $p^t$,z=(n−m)! mod $p^t$,那么答案就是x×reverse(y,$p^t$)×reverse(z,$p^t$)(reverse(a,b)表示计算a对b的乘法逆元)。那么下面问题就转化成如何计算n! mod pt。例如p=3,t=2,n=19

n!=1×2×3×4×5×6×7×8×…×19=(1×2×4×5×7×8×…×16×17×19)×(3×6×9×12×15×18)=(1×2×4×5×7×8×…×16×17×19)×36×(1×2×3×4×5×6)
后面的部分恰好是(n/p)!，于是递归即可。前半部分是以pt为周期的(1×2×4×5×7×8)≡(10×11×13×14×16×17)(mod 9)。下面是孤立的19，可以知道孤立出来的长度不超过 pt,于是直接计算即可。对于最后剩下的36这些数我们只要计算出n!,m!,(n−m)!里含有多少个p(不妨设x,y,z)，那么x−y−z就是$C_m^n$中p的个数，直接计算就行。

LL CRT(int n,LL* a,LL* m){
    LL M=1,p=0;
    for(int i=1;i<=n;i++) M=M*m[i];
    for(int i=1;i<=n;i++){
        LL w=M/m[i],x,y;
        exgcd(w,m[i],x,y);
        p=(p+a[i]*w*x%mod)%mod;
    }
    return (p%mod+mod)%mod;
}
LL calc(LL n,LL x,LL P){
    if(!n) return 1;
    LL s=1;
    for(int i=1;i<=P;i++)if(i%x) s=s*i%P;
    s=Pow(s,n/P,P);
    for(int i=n/P*P+1;i<=n;i++)if(i%x) s=s*i%P;
    return s*calc(n/x,x,P)%P;
}
LL multilucas(LL m,LL n,LL x,LL P){
    int cnt=0;
    for(int i=m;i;i/=x) cnt+=i/x;
    for(int i=n;i;i/=x) cnt-=i/x;
    for(int i=m-n;i;i/=x) cnt-=i/x;
    return Pow(x,cnt,P)%P*calc(m,x,P)%P*inverse(calc(n,x,P),P)%P*inverse(calc(m-n,x,P),P)%P;
}
LL exlucas(LL m,LL n,LL P){
    int cnt=0;
    LL p[20],a[20];
    for(LL i=2;i*i<=P;i++){
        if(P%i==0){
            p[++cnt]=1;
            while(P%i==0) p[cnt]=p[cnt]*i,P/=i;
            a[cnt]=multilucas(m,n,i,p[cnt]);
        }
    }
    if(P>1) p[++cnt]=P,a[cnt]=multilucas(m,n,P,P);
    return CRT(cnt,a,p);
    return 0;
}

(3)Baby Step Giant Step算法
具体步骤如下：
1.先把x=i*m-j，其中m=ceil(sqrt(C))，（ceil是向上取整）。
->这样原式就变为

$$A^{i*m-j}≡B(mod\,C)$$ ，
->再变为 $$A^{j}×B≡A^{m*i} (mod\,C)$$ 。
2.枚举j(范围0-m),将A^j×B存入hash表
3.枚举i(范围1-m),从hash表中寻找第一个满足 $$A^j×B=A^{m*i} (mod\,C)$$

此时x=i*m-j即为所求。

->在网上看到的其他题解大多用的是x=i*m+j，也可以做，只是会牵扯的求逆元，所以比较麻烦。使x=i*m-j就可以轻松避免这个问题了。

至于为什么只要枚举到m呢？

->x=i*m-j 也就是x 的最大值不会超过p,那超过p的怎么办 ？
有一个公式

$$a^{k\,mod\,p-1}=a^k (mod\,p)$$
这个公式的推导需要用到费马小定理
k mod p-1可以看做 k-m（p-1） ,原式可化成 $$\frac{a^k}{(a^{p-1})^m}≡a^k (mod\, p)$$
根据费马小定理 $$a^{p-1}≡1 (mod\,p)$$ 其中p为质数 ，a,p 互质，可得 $$\frac{a^k}{1^m}≡a^k (mod\,p)$$
$$a^k≡a^k (mod p)$$ 得证。
以下是求解 $$A^x≡B(mod C)$$ 的代码：

void solve(int A,int B){
    int m=ceil(sqrt(mod));
    if(A%mod==0){
        if(!B) printf("%d\n",mod!=1);
        else printf("Math Error\n");return;
    }
    ans=1;
    mp.clear();
    for(int i=0;i<=m;i++){
        if(i==0){
            ans=B%mod;mp[ans]=i;continue;
        }
        ans=(ans*A)%mod;
        mp[ans]=i;
    }
    ans=1;
    LL p=-1,AA=Pow(A,m);
    for(int i=1;i<=m;i++){
        ans=(ans*AA)%mod;
        if(mp[ans]){
            p=i*m-mp[ans];
            p=(p%mod+mod)%mod;
            break;
        }
    }
    if(p==-1) printf("Math Error\n");
    else printf("%lld\n",p);
}

(4)扩展BSGS算法

void exBSGS(LL A,LL B){
    LL d;
    LL g=0,k=1;
    while((d=gcd(A,B))!=1){
        if(B%d){printf("Math error\n");return;}
        g++;mod/=d;B/=d;k=k*A/d%mod;
        if(B==k){printf("%lld\n",g);return;}
    }
    LL m=ceil(sqrt(mod));
    mp.clear();
    LL ans=1,AA=Pow(A,m);
    for(int i=0;i<=m;i++){
        if(i==0){ans=B%mod,mp[ans]=i;continue;}
        ans=ans*A%mod;mp[ans]=i;
    }
    bool flag=0;
    for(int i=1;i<=m;i++){
        k=k*AA%mod;
        if(mp[k]){
            printf("%lld\n",i*m-mp[k]+g);
            flag=1;
            break;
        }
    }
    if(!flag) printf("Math error\n");
}

^_^