题目就是求C(n,k)*H(n - k)%m
0<= k<= n <=10^9, 1 <= m <= 10^5, n != 0
其中H(n)是错排第n项。
对于C(n,k)%m可以参考我以前的文章
对于H(n)
直接套公式:
可以发现肯定在某一位会出现后面都是模完都是0
#include <cstdio> #include <map> using std::pair; #define Pair pair<long long,long long> const long long M = 100000 + 9; long long fac[M],a[100],r[100]; struct triple{long long x,y,z;triple(const long long _x,const long long _y,const long long _z):x(_x),y(_y),z(_z){}}; triple exgcd(const long long a,const long long b) { if (!b) return triple(1,0,a); triple last(exgcd(b,a%b)); return triple(last.y,last.x - a/b * last.y,last.z); } long long CRT(const long long (&a)[100],const long long (&r)[100],const long long cnt) { long long res = 0,MM = 1; for (long long i = 1; i <= cnt; ++i) MM *= a[i]; for (long long i = 1; i <= cnt; ++i) (res += exgcd(MM / a[i],a[i]).x % a[i] * r[i] % MM * (MM / a[i]) % MM) %= MM; //printf("CRT %I64d\n",res + MM); return (res + MM) % MM; } long long power(long long n,long long k,long long MOD) { n %= MOD; long long res = 1; for (; k; n = n * n % MOD,k /= 2) if (k & 1) res = res * n % MOD; return res; } long long H(const long long n,const long long MOD) { //H(n) = (3 * ... * n) - (4 * ... * n) + (5 * ... * n) + ... + (-1)^(n-1) * n + (-1)^n if (n == 0) return 1 % MOD; if (n == 1) return 0; long long res = (n%2)?(-1):1; for (long long i = n,t = res; i >= 3; --i) { if (i % MOD == 0) break; (t = - t * i) %= MOD; (res += t) %= MOD; } //printf("TEST: %I64d\n",res + MOD); return (res + MOD) % MOD; } Pair FnModP(long long n,const long long p,const long long MOD) { //Fn = n! //fac[n] = n! % p long long res = 1; long long c = 0; while (n) { (res *= power(fac[MOD],n / MOD,MOD)) %= MOD; (res *= fac[n % MOD]) %= MOD; n /= p; c += n; } return std::make_pair(c,res); } void calc_fac(const long long p,const long long MOD) { fac[0] = 1; for (long long i = 1; i <= MOD; ++i) if (i % p) fac[i] = fac[i - 1] * i % MOD; else fac[i] = fac[i - 1]; } long long C(const long long n,const long long K,const long long p,const long long MOD) { //nCK % p^c calc_fac(p,MOD); Pair a(FnModP(n,p,MOD)),b(FnModP(K,p,MOD)),c(FnModP(n - K,p,MOD)); return 1ll * power(p,a.first - b.first - c.first,MOD) * a.second % MOD * ((exgcd(1ll * b.second * c.second % MOD,MOD).x % MOD + MOD) % MOD) % MOD; } long long work(long long n,long long K,long long MOD) { long long cnt = 0,c = 0; const long long m = MOD; for (long long i = 2; i * i <= MOD; ++i) if (MOD % i == 0) { a[++cnt] = 1; for (c = 0; MOD % i == 0; MOD /= i) ++c,a[cnt] *= i; r[cnt] = C(n,K,i,a[cnt]); } if (MOD > 1) r[++cnt] = C(n,K,MOD,MOD),a[cnt] = MOD; //for (long long i = 1; i <= cnt; ++i) printf("%I64d %I64d\n",a[i],r[i]); return 1ll * CRT(a,r,cnt) * H(n - K,m) % m; } int main() { #ifndef ONLINE_JUDGE freopen("SEQN.in","r",stdin); freopen("SEQN.out","w",stdout); #endif long long T,n,K,MOD; scanf("%I64d",&T); for (long long i = 1; i <= T; ++i) { scanf("%I64d%I64d%I64d",&n,&K,&MOD); printf("Case %I64d: %I64d\n",i,(work(n,K,MOD)%MOD + MOD)%MOD); } }