[BZOJ 1300] 大数计算器 同余

题目大意

给你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));
}