hdu 4373 组合数学

http://acm.hdu.edu.cn/showproblem.php?pid=4373

哎，多校的题。感觉是个很简单的数学问题，但是比赛的时候肿么都想不出来。今天决心把它给做了。

哎，原来以为是一道简单数学题呢，还用到Lucas,还用到中国剩余定理啊，哎。。。

于是我就去研究研究了一下中国剩余定理。终于懂了。参见黑书p230页。

怎么描述中国剩余定理的解法呢，我想想。

对于两个数m，n（m和n互质），满足x = a1 （mod m)和 x= a2 （mod n）的x为：

n * （n关于模m的逆元 * a1） + m *（m关于模n的逆元 * a2）= x；

转的，他写得好：http://blog.csdn.net/cyberzhg/article/details/7877465
题意：
m个for循环嵌套，有两种形式，第一类从1开始到n，第二类从上一层循环当前数开始到n，第一层一定是第一种类型，问总的循环的次数对364875103取余的结果。

首先可以看出，每一个第一类循环都是一个新的开始，与前面的状态无关，所以可以把m个嵌套分为几个不同的部分，每一个部分由第一类循环开始，最终结果相乘就可以。
剩下的就是第二类循环的问题，假设一个m层循环，最大到n，
只有第一层：循环n次。C(n, 1)
只有前两层：循环n + (n - 1) + ... + 1 = (n + 1) * n / 2 = C(n + 1, 2)
只有前三层：
(n + 1) * n / 2 + n * (n - 1) / 2 + ... + 1
= (1 + 2 + ... + n) / 2 + (1 + 4 + ... + n^2) / 2
= (n + 1) * n / 2 / 2 + n * (n + 1) * (2n + 1) / 6 / 2
= (n + 2) * (n + 1) * n / 6 = C(n + 2, 3)
……
找规律得第m层：C(n + m - 1, m)

接下来的问题就是如何求出这些东西对364875103的模，由于n和m都非常大，暴力统计必然是不行的。
364875103 = 97 * 3761599 （比赛时候哪有心情拆数玩，这个太不厚道了）
Lucas(n, m, p) = c(n % p,m % p) * Lucas(n / p, m / p, p);
用Lucas定理分别求出两个质数的余数，然后利用中国剩余定理求出对364875103的余数。

/*
Pro: 0

Sol:

date:
*/
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <queue>
#include <set>
#include <vector>
#define MOD1 3761599
#define MOD2 97
#define MOD 364875103
using namespace std;
int t,lop_len,lop_num,fir_num,a[100009];
__int64 ans,p1[3761610],p2[110],c1,c2;
__int64 pow_mod(__int64 x, __int64 y, __int64 mod){
    __int64 res = 1;
    while(y){
        if(y & 1)
            res =( res * x ) % mod;
        x = (x * x) % mod;
        y >>= 1;
    }
    return res;
}
__int64 cm(__int64 m, __int64 n, __int64 mod, __int64 p[]){//C(m,n),m取n
    if(m < n) return 0;
    __int64 res = (p[m] * pow_mod(p[n], mod - 2, mod) ) % mod;
    return ( res * pow_mod(p[m - n] , mod - 2, mod)) % mod;
}

void init(){
    p1[0] = 1;  p2[0] = 1;
    for(int i = 1; i <= MOD1; i ++){
        p1[i] = p1[i - 1] * i % MOD1;
    }
    for(int i = 1; i <= MOD2; i ++){
        p2[i] = p2[i - 1] * i % MOD2;
    }
    c1 = MOD2 * pow_mod(MOD2, MOD1 - 2, MOD1);
    c2 = MOD1 * pow_mod(MOD1, MOD2 - 2, MOD2);
}
__int64 lucas(__int64 m, __int64 n, __int64 mod,__int64 p[]){
    __int64 res = 1;
    while(n && m && res){
        res = (res * cm(m % mod,n % mod,mod,p));
        n /= mod;
        m /= mod;
    }
    return res;
}
int main(){
//3761599  97 拆数
//    for(int i = 2; i <= 1000; i ++)
//        if(364875103 % i == 0)
//            printf("%d  %d",364875103 / i,  i);
    init();
    //求c1，c2和记录
    //c1 为mod2 * （mod2 对 mod1 的逆元），即c1 % mod1 == 1， c2 % mod2 == 0；
    //c2 为mod1 * （mod1 对 mod2 的逆元），即c2 % mod2 == 1， c1 % mod1 == 0；
    scanf("%d",&t);
    for(int ca = 1; ca <= t; ca ++){
        scanf("%d%d%d",&lop_len,&lop_num,&fir_num);
        for(int i = 0; i < fir_num; i ++)
            scanf("%d",a + i);//index of type 1, started from index 0
        a[fir_num] = lop_num;//这个，需要理解和注意
        ans = 1;
        for(int i = 0; i < fir_num; i ++){
            int tmp = (a[i + 1] - a[i]);
            __int64 m1 = lucas(tmp + lop_len - 1, tmp, MOD1, p1);
            __int64 m2 = lucas(tmp + lop_len - 1, tmp, MOD2, p2);
            __int64 mm = (m1 * c1 + m2 * c2) % MOD;
            ans = (ans * mm) % MOD;
        }
        printf("Case #%d: %I64d\n",ca,ans);
    }
	return 0;
}