HDU 4373 Mysterious For（Lucas定理、中国剩余定理）

题目链接：
HDU 4373 Mysterious For
题意：
定义两种循环：
$for(int\ a[i] = 0; a[i] < n; ++a[i]) { \\...\\ }$
$for(int\ a[i] = a[i-1]; a[i] < n; ++a[i]) {\\ ... \\}$
会有$num$层循环属于第一种循环，给出这些循环的位置，求总共的循环次数是多少？答案$mod\ 364875103$
分析：
首先我们考虑第一种循环下面跟$k-1$个第二种循环，总共是$k$层循环的情况。定义在$t$层循环到$i$的次数为$f_t(i)$次

• 对于第一层循环，循环到$i$，显然有 $$f_1(i)=1\\ 总的循环次数是：\sum{f_1(i)}=n=C_{n}^{1}$$
• 对于第二层循环，因为第二层循环属于第二种循环，想要循环到$i$,那么相应的第一层循环就要$\leq i$,也就是 $$f_2(i)=\sum_{j\leq i}f_1(j)=i\\ 总的循环次数是:\sum{f_2(i)}=\frac{(n+1)*n}{2}=C_{n+1}^{2}$$
• 对于第三层循环，同理可得： $$f_3(i)=\sum_{j\leq i}{f_2(i)}=\frac{(i+1)*i}{2}\\ 总的循环次数是：\sum{f_3(i)}=\sum{\frac{1}{2}*(i^2+i)}=\frac{(n+2)*(n+1)*(n)}{6}=C_{n+2}^{3}$$
• …
• 可以发现对于第$k$层循环总的循环次数是 $$C_{n+k-1}^{k}$$

接着我们需要计算$C_{n+k-1}^{k}\ \% 364875103的值$.比较坑爹的是$364875103$不是一个素数！但是对$364875103$进行素数检查时可以发现$364875103=97*3761599，$而这两个因子都是素数，令$mod1 = 97,mod2 = 3761599,mod=364875103$
令$m=k,n=n+k-1$,我们

$$把C_{n}^{m}\ \ \% \ \ mod的值记为x,\\ 把C_{n}^{m}\ \%\ mod1的值记为x1,\\ 把C_{n}^{m}\ \%\ mod2 的值记为x2$$
则有：
$$x \equiv \ x_1(\% \ \ mod1)\\ x \equiv \ x_2(\% \ \ mod2)$$
因为 $mod1,mod2$都是素数，所以 $x_1和x_2$都可以用 $Lucas定理求解出来。$那么剩下的问题就是求解上面的 同余方程了。
利用 中国剩余定理求解同余方程。
我们定义： $$inv1为:(mod2在模mod1域下的逆元)*mod2\\ inv2为:(mod1在模mod2域下的逆元)*mod1$$
那么答案就是:
$$x=(inv1*x_1\ + inv2 * x_2)\% \ mod$$
最后在根据乘法原理易知每段求得的答案需要累乘。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <string>
#include <algorithm>
#include <climits>
#include <cmath>
#include <ctime>
#include <cassert>
#define IOS ios_base::sync_with_stdio(0); cin.tie(0);
using namespace std;
typedef long long ll;
const int mod1 = 97;
const int mod2 = 3761599;
const ll MOD = 364875103;
const int MAX_M = 100010;

ll fac1[mod1 + 10], fac2[mod2 + 10], e1, e2;

ll quick_pow(ll a, ll b, ll mod)
{
    ll res = 1, tmp = a;
    while(b) {
        if(b & 1) res = res * tmp % mod;
        tmp = tmp * tmp % mod;
        b >>= 1;
    }
    return res;
}

void init() 
{
    fac1[0] = fac2[0] = 1;
    for(int i = 1; i < mod1; ++i) { fac1[i] = fac1[i - 1] * i % mod1; }
    for(int i = 1; i < mod2; ++i) { fac2[i] = fac2[i - 1] * i % mod2; }
    e1 = quick_pow(mod2, mod1 - 2, mod1) * mod2;
    e2 = quick_pow(mod1, mod2 - 2, mod2) * mod1;
}

ll ex_gcd(ll a, ll b, ll& x, ll& y) 
{
    if(b == 0) {
        x = 1, y = 0;
        return a;
    }
    ll r = ex_gcd(b, a % b, y, x);
    y -= a / b * x;
    return r;
}

ll inv(ll a, ll n)
{
    ll x, y;
    ll d = ex_gcd(a, n, x, y);
    if(d == 1) return (x % n + n) % n;
    else return 0;
}

ll C(int n, int m, int id, int mod)
{
    if(m > n) return 0;
    ll a, b, c;
    if(id == 1) {
        a = fac1[n], b = fac1[n - m], c = fac1[m];
    } else {
        a = fac2[n], b = fac2[n - m], c = fac2[m];
    }
    //return a * inv(b * c % mod, mod) % mod;  //也可以用扩展欧几里德求逆元
    return a * quick_pow(b * c % mod, mod - 2, mod) % mod;
}

ll Lucas(int n, int m, ll mod, int id)
{
    if(m == 0) return 1;
    return C(n % mod, m % mod, id, mod) * Lucas(n / mod, m / mod, mod, id) % mod; 
}

int main()
{
    init();
    int T, n, m, num, data[30], cases = 0;
    scanf("%d", &T);
    while(T--) {
        scanf("%d%d%d", &n, &m, &num);
        for(int i = 0; i < num; ++i) { scanf ("%d", &data[i]); }
        data[num++] = m;
        ll ans = 1;
        for(int i = 1; i < num; ++i) {
            int k = data[i] - data[i - 1];
            ll t1 = Lucas(n + k - 1, k, mod1, 1) * e1 % MOD;
            ll t2 = Lucas(n + k - 1, k, mod2, 2) * e2 % MOD;
            ans = ans * (t1 + t2) % MOD;
        }
        printf("Case #%d: %lld\n", ++cases, ans);
    }
    return 0;
}