Description

【题意】
求${\mathrm{C}}_n^m \bmod{p}$ ， 其中 $\mathrm{C}$ 为组合数。

【输入格式】
一行三个整数 $n,m,p$ ，含义由题所述。

【输出格式】
一行一个整数，表示答案。

【样例输入 #1】
5 3 3

【样例输出 #1】
1

【样例输入 #2】
666 233 123456

【样例输出 #2】
61728

【提示】
对于 $50 \%$ 的数据，$1 \le m \le n \le {10}^{4}$，$2 \le p \le {10}^6$，**不保证** $p$ 是质数。
对于 $100 \%$ 的数据，$1 \le m \le n \le {10}^{18}$，$2 \le p \le {10}^6$，**不保证** $p$ 是质数。

Hint

by hansang:

/*  by；hansang 
求的是 C nm %P (P不一定是质数)，把 C nm 设为 x

把 P 化为 p1^k1p2^k2... (pi是质数)
列出式子：
x % (p1^k1) = a1
x % (p2^k2) = a2
......
用中国剩余定理合并可得答案
/
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N=1e6+10;
LL a[N], c[N]; int cnt; //a 数组表示 x (C nm ) % (pi^ki) 的值，c 数组表示 pi^ki，cnt 表示 a c 数组长度
LL q_pow(LL a, LL b, LL P) //求 a ^b %P
{
LL ans=1%P; a%=P; //防止 P 为 1，a 先%P
for(; b; b>>=1) //枚举 b 的二进制位
{
if(b&1) ans=ansa%P; //b 这个位为 1，ans 就乘这个位的 a
a=aa%P; // a=a^2
}
return ans;
}
void exgcd(LL a, LL b, LL &d, LL &x, LL &y) //扩展欧几里得
{
if(b==0) {d=a; x=1; y=0;} //找到最大公约数了，d=a
else
{
exgcd(b, a%b, d, y, x);
y-=(a/b)x;
/
设新的 x=dx，y=dy
则有；
ax +by = btx +(a%b)ty
已知；
a%b = a-[a/b]b ([]为下取整)
则；
ax +by
= btx +(a%b)ty
= btx +(a-[a/b]b)ty
= btx +aty -[a/b] bty
= aty +b*(tx-[a/b]ty)
所以；
ty=x，tx=(y+[a/b]x)
倒推回来 y 就要减 (a/b)x
/
}
}
// exgcd 和 q_pow 都为很重要的板子，请熟背
LL inv(LL a, LL P) //逆元 (运用 exgcd)
{
LL A=a, B=P, K=1, d, x, y; // ax=1 (%P) = ax+P?=1 = Ax+By=K
exgcd(A, B, d, x, y); //算的是 Ax+By=d (d为 gcd(a, b))
x=x(K/d); LL dx=abs(B/d); x=(x%dx+dx)%dx; //把 x变成最小正整数解
// x(K/d)变成 Ax+By=K 的解，dx等于当前 x 与下一组解 x 的差 (下一组解x+abs(B/d)，y-abs(A/d))
// x=(x%dx+dx)%dx 是常用的把 x变成最小正整数解的方式，如果还不懂可以手模一下
return x;
}
LL fac(LL n, LL P, LL Pk) // 算 n! 不含 P 的部分 (例如 P=3，则 fac(9)=12*(1)45*(2)78*(1))
{
LL ans=1; if(n0) return 1; //ans初始化为 1，n=0 值为 1
for(LL i=1; i<Pk; i++) if(i%P!=0) ans=(ansi)%Pk; //ans等于 1 ~ Pk 里不含 P 的数的积
ans=q_pow(ans, n/Pk, Pk); // ans^(n/Pk) % Pk， ans等于 1 ~ n/PkPk 里不含 P 的数的积
for(LL i=1; i<=n%Pk; i++) if(i%P!=0) ans=(ansi)%Pk; //ans等于 1 ~ n 里不含 P 的数的积
return ansfac(n/P, P, Pk)%Pk; //继续递归，因为有些数 = x* P^i，x 也要算进来
}
LL C(LL n, LL m, LL P, LL Pk) // C n m
{
if(n<m) return 0; // n<m 值为 0
LL f1=fac(n, P, Pk), f2=fac(m, P, Pk), f3=fac(n-m, P, Pk), sum=0;
// f1表示 n!中不含 P的部分， f2表示 m!中不含 P的部分，f3表示 (n-m)!中不含 P的部分，
for(LL i=n; i; i/=P) sum+=i/P; //sum加上 n 里面有多少个 P
for(LL i=m; i; i/=P) sum-=i/P; //sum减去 m 里面有多少个 P
for(LL i=n-m; i; i/=P) sum-=i/P; //sum减去 (n-m) 里面有多少个 P
return f1inv(f2, Pk)%Pkinv(f3, Pk)%Pkq_pow(P, sum, Pk)%Pk;
// f1乘上 f2的逆元 (%Pk) 乘上 f3的逆元 (%Pk) 再乘上 P^sum (%Pk) 就是 C n m
}
LL CRT() //中国剩余定理
{
LL m=1, ans=0; // m 表示所有 (pi^ki) 的最小公倍数 (就是 P)，ans初始化为 0
for(int i=1; i<=cnt; i++) m=c[i];
for(int i=1; i<=cnt; i++)
{
ans=(ans+a[i](m/c[i])%minv(m/c[i], c[i])%m)%m;
// (m/c[i]) 是除 c[i] 其它所有 (pi^ki) 的公倍数，加上对其它 c的取模结果不影响
// inv(m/c[i], c[i]) 是 (m/c[i]) 的逆元 (%c[i])，乘以 (m/c[i]) %c[i] 等于 1
// 再乘以 a[i] 就达成了不改变其它 c的取模结果，让 %c[i] 的结果等于 a[i]
// %m 防止溢出，加减 m 对所有 c的取模结果没影响 (因为是公倍数)
}
return ans;
}
LL exLucas(LL n, LL m, LL P) //拓展卢卡斯
{
cnt=0; //cnt 初始化
for(int i=2; i*i<=P; i++) if(P%i0) //找出组成 P 的质数
{
LL tmp=1;  while(P%i==0) {tmp*=i; P/=i;} //把 P中的 i提取到 tmp里
cnt++;  a[cnt]=C(n, m, i, tmp); c[cnt]=tmp; //存到 a c 数组里
}
if(P>1) {cnt++; a[cnt]=C(n, m, P, P); c[cnt]=P;} //当 P 还剩有质数时
return CRT();
}
int main()
{
LL n, m, P; scanf("%lld%lld%lld", &n, &m, &P);
printf("%lld\n", exLucas(n, m, P));
return 0;
}

by scy:</p>

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N=110;
LL cnt,a[N],p[N],pk[N];
void exgcd(LL a,LL b,LL &d,LL &x,LL &y)
{
	if(b==0){d=a;x=1;y=0;}
	else
	{
		exgcd(b,a%b,d,y,x);
		y-=(a/b)*x;
	}
}
LL qpow(LL a,LL b,LL P)
{
	LL ans=1%P;a%=P;
	for(;b;b>>=1)
	{
		if(b&1)ans=ans*a%P;
		a=a*a%P;
	}
	return ans;
}
LL inv(LL a,LL P)
{
	LL A=a,B=P,x,y,d;exgcd(A,B,d,x,y);
	LL dx=abs(B/d);
	x=(x%dx+dx)%dx;
	return x;
}
LL fac(LL n,LL p,LL pk)
{
	if(n==0) return 1;
	LL ans=1;for(LL i=1;i<pk;i++)if(i%p!=0) ans=ans*i%pk;
	ans=qpow(ans,n/pk,pk);
	for(int i=1;i<=n%pk;i++)if(i%p!=0) ans=ans*i%pk;
	return ans*fac(n/p,p,pk) %pk;
}
LL C(LL n,LL m,LL p,LL pk)
{
	if(n<m) return 0;
	LL f1=fac(n,p,pk),f2=fac(m,p,pk),f3=fac(n-m,p,pk);
	LL xyz=0;
	for(LL i=n;i;i/=p)xyz+=i/p;
	for(LL i=m;i;i/=p)xyz-=i/p;
	for(LL i=n-m;i;i/=p)xyz-=i/p;
	return f1*inv(f2,pk)%pk*inv(f3,pk)%pk * qpow(p,xyz,pk) %pk;
}
LL exlucas(LL n,LL m,LL P)
{
	cnt=0;LL tp=P;
	for(LL i=2;i*i<=tp;i++) if(tp%i==0)
	{
		cnt++;p[cnt]=i;pk[cnt]=1;while(tp%i==0)tp/=i,pk[cnt]*=i;
	}
	if(tp>1){cnt++;p[cnt]=pk[cnt]=tp;};
	
	for(int i=1;i<=cnt;i++) a[i]=C(n,m,p[i],pk[i]);
	
	LL x=0;
	for(int i=1;i<=cnt;i++)x=( x+ a[i] *  (P/pk[i]) % P * inv(P/pk[i], pk[i]) % P )%P;
	return x;
}
int main()
{
	LL n,m,p;scanf("%lld%lld%lld",&n,&m,&p);
	printf("%lld\n",exlucas(n,m,p));
}