题目大意
给你N和M, 要求计算
C(N, M)
如果答案10进制表示下大于12位, 只需要输出’
最高3位...最低9位
‘即可, 否然输出一个整数表示答案。
[latex] 0 \le M \le N \le 10^6[/latex]
解题报告
这算是一道处理同余和大数的技巧整合题吧…
使用lg算出最高三位, 然后最低9位使用中国剩余定理合并处理。
处理模[latex]P^K[/latex]意义下的答案时只需要把运算中每一个数的质因数分解表示法下P的指数抽出来单独处理即可。
注意处理逆元时使用欧拉定理, 不可以乱用费马小定理。
主要看代码:
代码
CRT部分略丑。
#include <cstring> #include <cstdlib> #include <cstdio> #include <algorithm> #include <climits> #include <cmath> #include <vector> using namespace std; typedef long long LL; const int MAXN = 1000010; const long long LIM = 1000000000000LL; inline int advPow (LL a, LL b, int c) { a %= c; int ret = 1; while(b) { if(b&1) ret = (LL)ret*a%c; a = (LL)a*a%c; b >>= 1; } return ret; } inline LL mult (LL a, LL b, LL mod) { a %= mod, b %= mod; LL ret = 0; while(b) { if(b&1) ret = (ret+a)%mod; a = (a+a)%mod; b >>= 1; } return ret; } LL CRT (vector<pair<int, int> > vec) { // fir -> ans, sec->mod LL prod = 1; for(int i = 0; i<(int)vec.size(); i++) prod *= vec[i].second; LL ans = 0; for(int i = 0; i<(int)vec.size(); i++) { LL M = prod/vec[i].second; int x; if(vec[i].second == 4096) x = 4096-4096/2-1; else x = vec[i].second-vec[i].second/5-1; ans = (ans+mult(M, advPow(M, x, vec[i].second)%prod*vec[i].first, prod))%prod; } return ans; } double flog[MAXN], rest; int calHi (int N, int M) { flog[0] = 0; for(int i = 1; i<=N; i++) flog[i] = flog[i-1]+log10(i); double result = flog[N]-flog[M]-flog[N-M]; rest = result; return pow(10.0, result-floor(result-2.0)); } int N, M; pair<int, int> expr[MAXN]; int solve (pair<int, int> pr) { int P = pr.first, K = pr.second; int mod = 1; for(int i = 1; i<=K; i++) mod *= P; for(int i = 1; i<=N; i++) { expr[i] = make_pair(i, 0); while(expr[i].first%P == 0) expr[i].first /= P, expr[i].second++; expr[i].first %= mod; } pair<int, int> prd(1, 0); if(N-M > N/2) M = N-M; for(int i = M+1; i<=N; i++) { prd.first = (LL)prd.first*expr[i].first%mod; prd.second += expr[i].second; } for(int i = 1; i<=N-M; i++) { //advPow(expr[i].first, mod-mod/P-1, mod)*expr[i].first%mod prd.first = (LL)prd.first*advPow(expr[i].first, mod-mod/P-1, mod)%mod; prd.second -= expr[i].second; } return (LL)prd.first*advPow(P, prd.second, mod)%mod; } int main () { scanf("%d%d", &N, &M); int hi = calHi(N, M); if(rest > 12) printf("%d...", hi); LL t = LIM; vector<pair<int, int> > vec, sol; for(LL i = 2; i*i<=t; i++) { if(t%i == 0) { pair<int, int> pr(i, 0); while(t%i == 0) {pr.second++, t/=i;} vec.push_back(pr); } } for(int i = 0; i<(int)vec.size(); i++) { LL x = 1; for(int j = 1; j<=vec[i].second; j++) x *= vec[i].first; sol.push_back(make_pair(solve(vec[i]), x)); } if(rest > 12) printf("%09I64d\n", CRT(sol)%1000000000); else printf("%I64d\n", CRT(sol)); }
Join the discussion