【拓展Lucas】模板

求\(C_n^m \mod p\)，写得太丑了qwq。
第一次写拓展Lucas竟然是在胡策的时候qwq写了两个半小时啊_(:з」∠)_还写挂了一个地方qwq
当然今天胡策我也是第一次写中国剩余定理(ˇˍˇ)
↑平时懒得动手的后果→_→

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;

int ipow2(int a, int b) {
	int ret = 1, w = a;
	while (b) {
		if (b & 1) ret = ret * w;
		w = w * w;
		b >>= 1;
	}
	return ret;
}

int ipow(int a, int b, int c) {
	int ret = 1, w = a;
	while (b) {
		if (b & 1) ret = 1ll * ret * w % c;
		w = 1ll * w * w % c;
		b >>= 1;
	}
	return ret;
}

int pic;

const int N = 100003;

struct data {int num, zhi;};

data work(int n, int p, int c) {
	if (n == 1 || n == 0) return (data) {1, 0};
	int ret1, las = n % pic, ret2 = 1;
	
	for (int i = 1; i <= las; ++i) {
		if (i % p == 0) continue;
		ret2 = 1ll * ret2 * i % pic;
	}
	
	if (n >= pic) {
		ret1 = ret2;
		for (int i = las + 1; i < pic; ++i) {
			if (i % p == 0) continue;
			ret1 = 1ll * ret1 * i % pic;
		}
		ret1 = ipow(ret1, n / pic, pic);
	} else
		ret1 = 1;
	
	
	int tim = n / p;
	data s = work(n / p, p, c);
	return (data) {1ll * ret1 * ret2 % pic * s.num % pic, tim + s.zhi};
}

void exgcd(int a, int b, ll &x, ll &y) {
	if (b == 0) {x = 1; y = 0;}
	else {
		exgcd(b, a % b, x, y);
		ll t = x;
		x = y;
		y = t - a / b * y;
	}
}

int inv(int a, int b) {
	ll x, y;
	exgcd(a, b, x, y);
	return (x % b + b) % b;
}

int Lucas(int n, int m, int p, int c) {
	data s1, s2, s3;
	pic = ipow2(p, c);
	s1 = work(n, p, c);
	s2 = work(m, p, c);
	s3 = work(n - m, p, c);
	
	int tim = s1.zhi - s2.zhi - s3.zhi;
	int num1 = s1.num, num2 = s2.num, num3 = s3.num;
	num2 = inv(num2, pic); num3 = inv(num3, pic);
	num1 = 1ll * num1 * num2 % pic * num3 % pic;
	for (int i = 1; i <= tim; ++i)
		num1 = 1ll * num1 * p % pic;
	return num1;
}

int prime[N], cnt = 0, cc[N], pc[N];

void Div(int p) {
	int sq = ceil(sqrt(p));
	for (int i = 2; i <= sq; ++i) {
		if (p % i == 0) {
			prime[++cnt] = i;
			pc[cnt] = 1;
			cc[cnt] = 0;
			while (p % i == 0) {++cc[cnt]; p /= i; pc[cnt] *= i;}
		}
	}
	if (p > 1) prime[++cnt] = p, cc[cnt] = 1, pc[cnt] = p;
}

int n, m, k, p, x;
int num[N];

int main() {
	int T; scanf("%d", &T);
	while (T--) {
		scanf("%d%d%d%d", &n, &m, &k, &p);
		
		cnt = 0;
		Div(p);
		
		for (int i = 1; i <= cnt; ++i)
			num[i] = Lucas(n, m, prime[i], cc[i]);
		
		x = 0;
		for (int i = 1; i <= cnt; ++i)
			(x += (1ll * (p / pc[i]) * inv((p / pc[i]) % pc[i], pc[i]) % p * num[i] % p)) %= p;
		printf("%d\n", x);
	}
	
	return 0;
}