codeforces100633J Ceizenpok’s formula 【扩展Lucas】

题目链接

Ceizenpok’s formula

题解

求\(C_n^k \pmod m\)
扩展\(Lucas\)板题

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define LL long long int
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define BUG(s,n) for (int i = 1; i <= (n); i++) cout<<s[i]<<' '; puts("");
#define res register
using namespace std;
const int maxn = 1000005,maxm = 100005,INF = 1000000000;
inline LL read(){
	LL out = 0,flag = 1; char c = getchar();
	while (c < 48 || c > 57){if (c == '-') flag = -1; c = getchar();}
	while (c >= 48 && c <= 57){out = (out << 3) + (out << 1) + c - 48; c = getchar();}
	return out * flag;
}
LL md;
LL pi[50],pk[50],cnt,fac[maxn];
void pre(LL pi,LL pk){
	fac[0] = 1;
	for (LL i = 1; i <= pk; i++)
		if (i % pi) fac[i] = 1ll * fac[i - 1] * i % pk;
		else fac[i] = fac[i - 1];
}
void init(){
	LL x = md;
	for (res LL i = 2; i * i <= x; i++)
		if (x % i == 0){
			++cnt; pi[cnt] = i; pk[cnt] = 1;
			while (x % i == 0) pk[cnt] *= i,x /= i;
		}
	if (x - 1) ++cnt,pi[cnt] = pk[cnt] = x;
}
inline LL qpow(LL a,LL b,LL md){
	LL ans = 1;
	for (; b; b >>= 1,a = 1ll * a * a % md)
		if (b & 1) ans = 1ll * ans * a % md;
	return ans % md;
}
void exgcd(LL a,LL b,LL& d,LL& x,LL &y){
	if (!b) {x = 1; y = 0; d = a;}
	else exgcd(b,a % b,d,y,x),y -= a / b * x;
}
inline LL inv(LL a,LL P){
	if (!a) return 0;
	LL d,x,y; exgcd(a,P,d,x,y);
	x = (x % P + P) % P; if (!x) x += P;
	return x;
}
inline LL Fac(LL n,LL pi,LL pk){
	if (!n) return 1;
	int ans = qpow(fac[pk],n / pk,pk);
	return 1ll * ans * fac[n % pk] % pk * Fac(n / pi,pi,pk) % pk;
}
inline LL C(LL n,LL m,int pi,int pk){
	if (m > n) return 0;
	int a = Fac(n,pi,pk),b = Fac(m,pi,pk),c = Fac(n - m,pi,pk),ans,k = 0;
	for (res LL i = n; i; i /= pi) k  += i / pi;
	for (res LL i = m; i; i /= pi) k  -= i / pi;
	for (res LL i = n - m; i; i /= pi) k  -= i / pi;
	ans = 1ll * a * inv(b,pk) % pk * inv(c,pk) % pk * qpow(pi,k,pk) % pk;
	return 1ll * ans * (md / pk) % md * inv(md / pk,pk) % md;
}
inline int exlucas(LL n,LL m){
	if (m > n) return 0;
	LL ans = 0;
	for (res LL i = 1; i <= cnt; i++)
		pre(pi[i],pk[i]),ans = (ans + C(n,m,pi[i],pk[i])) % md;
	return ans;
}
int main(){
	LL n = read(),k = read(); md = read();
	init();
	cout << exlucas(n,k) << endl;;
	return 0;
}