洛谷P4720 【模板】扩展卢卡斯

P4720 【模板】扩展卢卡斯

题目背景

这是一道模板题。

题目描述

求

C(n,m)%P

其中 C 为组合数。

输入输出格式

输入格式：

一行三个整数 n,m,p ，含义由题所述。 

输出格式：

一行一个整数，表示答案。

输入输出样例

输入样例#1：

5 3 3

输出样例#1：

1

输入样例#2：

666 233 123456

输出样例#2： 

61728

说明

1≤m≤n≤10¹⁸,2≤p≤1000000 ，不保证 p 是质数。

sol：ExLucas模板 可以做P不是质数的组合数

具体方法简单说下：把P分成若干个质因数的幂次的乘机，分别求出那时候的组合数，用Exgcd合并起来即可

口胡一下怎么快速求n！

假定n=19，模数为3²

n!=1∗2∗3∗⋯∗19

=(1*2*4*5*7*8*10*11*13*14*16*17*19)*3⁶*6!

≡(1∗2∗4∗5∗7∗8)²∗19∗3⁶∗6!

然后对于1*2*4*5*7*8和19这两段长度不超过3²暴力求解，6!递归求解，另一团快速幂。。。

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
inline ll read()
{
    ll s=0;
    bool f=0;
    char ch=' ';
    while(!isdigit(ch))
    {
        f|=(ch=='-'); ch=getchar();
    }
    while(isdigit(ch))
    {
        s=(s<<3)+(s<<1)+(ch^48); ch=getchar();
    }
    return (f)?(-s):(s);
}
#define R(x) x=read()
inline void write(ll x)
{
    if(x<0)
    {
        putchar('-'); x=-x;
    }
    if(x<10)
    {
        putchar(x+'0');    return;
    }
    write(x/10);
    putchar((x%10)+'0');
    return;
}
#define W(x) write(x),putchar(' ')
#define Wl(x) write(x),putchar('
')
const int N=65;
ll n,m;
ll cnt=0,Mo[N],Yu[N];
ll prim[N],T[N];

inline void Exgcd(ll a,ll b,ll &X,ll &Y);
inline ll Solve();
inline ll Ksm(ll x,ll y,ll Mod);
inline ll Inv(ll Num,ll Mod);
inline ll Jiecheng(ll n,ll P,ll Pk);
inline ll C(ll n,ll m,ll P,ll Pk);
inline ll CRT(ll Num,ll Mod,ll Pk);
inline ll ExLucas(ll n,ll m,ll Mod);

inline ll Ksm(ll x,ll y,ll Mod)
{
    ll ans=1;
    while(y)
    {
        if(y&1) ans=ans*x%Mod;
        x=x*x%Mod;
        y>>=1;
    }
    return ans;
}
inline void Exgcd(ll a,ll b,ll &X,ll &Y)
{
    if(b==0)
    {
        X=1;
        Y=0;
        return;
    }
    Exgcd(b,a%b,X,Y);
    ll XX=X,YY=Y;
    X=YY;
    Y=XX-a/b*YY;
    return;
}
inline ll Inv(ll Num,ll Mod) //Num*Inv = 1(%Mod)
{
    ll a,b,X,Y;
    a=Num;
    b=Mod;
    Exgcd(a,b,X=0,Y=0);
    X=(X%b+b)%b;
    return X;
}
inline ll Jiecheng(ll n,ll P,ll Pk)
{
    if(!n) return 1;
    ll i,ans=1;
    for(i=2;i<=Pk;i++) if(i%P) ans=ans*i%Pk;
    ans=Ksm(ans,n/Pk,Pk);
    for(i=2;i<=n%Pk;i++) if(i%P) ans=ans*i%Pk;
    return ans*Jiecheng(n/P,P,Pk)%Pk;
}
inline ll C(ll n,ll m,ll P,ll Pk)
{
    ll Jn=Jiecheng(n,P,Pk);
    ll Jm=Jiecheng(m,P,Pk);
    ll Jnm=Jiecheng(n-m,P,Pk);
    ll i,oo=0;
    for(i=n;i;i/=P) oo+=i/P;
    for(i=m;i;i/=P) oo-=i/P;
    for(i=n-m;i;i/=P) oo-=i/P;
    return Jn*Ksm(P,oo,Pk)%Pk*Inv(Jm,Pk)%Pk*Inv(Jnm,Pk)%Pk%Pk;
}
inline ll ExLucas(ll n,ll m,ll Mod)
{
    ll i,ans=0,tmp=Mod;
    for(i=2;i<=sqrt(tmp);i++) if(tmp%i==0)
    {
        ll Num=1;
        while(tmp%i==0)
        {
            Num*=i;
            tmp/=i;
        }
        Mo[++cnt]=Num;
        Yu[cnt]=C(n,m,i,Num);
    }
    if(tmp>1)
    {
        Mo[++cnt]=tmp;
        Yu[cnt]=C(n,m,tmp,tmp);
    }
    return Solve();
}
inline ll Solve()
{
    memmove(prim,Mo,sizeof Mo);
    memmove(T,Yu,sizeof T);
    int i;
    ll Lcm=prim[1];
    for(i=2;i<=cnt;i++)
    {
        ll a=prim[i-1],b=prim[i],c=T[i]-T[i-1],r=1,X,Y;
        Exgcd(a,b,X=0,Y=0);
        ll tmp=b/r;
        X=(X*c%tmp+tmp)%tmp;
        prim[i]=Lcm=Lcm*prim[i];
        T[i]=X*prim[i-1]+T[i-1];
    }
    return T[cnt];
}
int main()
{
    int i,Mod;
    R(n); R(m); R(Mod);
    Wl(ExLucas(n,m,Mod));
    return 0;
}
/*
input
5 3 3
output
1

input
666 233 123456
output
61728
*/