HDU-4373 Mysterious For（中国剩余定理+卢卡斯定理+可重组合）

题意

定义两种类型的 $for$ 循环：

for(int a[i]=0;a[i]<n;a[i]++){   //类型一
    ...
}
for(int a[i]=a[i-1];a[i]<n;a[i]++){   //类型二
    ...
}

现在给定循环范围 $n$ 并将 $m$ 重循环嵌套，最里层的循环内是 $f$ 函数，已知类型一的循环有 $k$ 个并给定这 $k$ 个循环的位置。求 $f$ 函数调用了多少次，次数模以 $364875103$ 。
$1 \leq n \leq 10^6$
$1 \leq m\leq 10^5$
$1 \leq k \leq 15$

思路

首先对 $364875103$ 进行质因数分解得到 $97 \times 3761599$ ，分别对这两个质数求模再用 $CRT$ 合并即可。
把一个类型一循环和下面与它连续的类型二循环看成一个整体。这里的循环变量从上到下满足单调不减，循环的次数事实上是一个可重组合，可重组合的 $n$ 就是第一层循环的长度 $n$ ，可重组合的 $m$ 是这个整体的循环语句数。
对于一个可重的组合 $H$ ，有 $H_n^m=C_{n+m-1}^m$
证明可采用构造法，设选的 $m$ 个数为 $a_1,a_2,a_3...a_m (a_i<=a_{i+1},a_i \in [1,n])$ ，再构造一个 $b$ 数组满足 $b_i=a_i+i-1$ ，这个用 $a$ 数组 表示 $b$ 数组得 $a_1,a_2+1,a_3+2...a_m+m-1$ ，不难发现 $a_i<a_{i+1}，a_i \in [1,n+m-1]$，这满足一个普通的组合 $C_{n+m-1}^m$，证毕。
而计算 $C_n^m \% P$ 时常常因为要进行阶乘的数大于 $P$ 导致阶乘均为 $0$ ，从而无法得出正确答案。这是 $Lucas$ 定理的好处就体现出来了。它能让计算时 $n,m$ 均小于 P P <script type="math/tex" id="MathJax-Element-91">P</script> ，防止出现以上的情况。

#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#define FOR(i,x,y) for(int i=(x);i<=(y);i++)
#define DOR(i,x,y) for(int i=(x);i>=(y);i--)
typedef long long LL;
using namespace std;
int _P[2]={97,3761599};
LL fac[2][3800003];

LL Pow(LL a,LL p,LL P)
{
    LL res=1;
    while(p)
    {
        if(p&1)(res*=a)%=P;
        (a*=a)%=P;
        p>>=1; 
    }
    return res;
}
void init()
{
    FOR(i,0,1)
    {
        LL f=1;
        FOR(j,0,_P[i]-1)
        {
            fac[i][j]=f;
            (f*=j+1)%=_P[i];
        }
    }
    return;
}

struct CRT
{
    LL A,P;
    LL gcd(LL a,LL b){return b?gcd(b,a%b):a;}
    void exgcd(LL a,LL b,LL &x,LL &y)
    {
        if(!b){x=1,y=0;return;}
        exgcd(b,a%b,y,x);y-=a/b*x;
        return;
    }
    LL Exgcd(LL A,LL B,LL C)
    {
        LL k1,k2;
        exgcd(A,B,k1,k2);
        return (k1*C%B+B)%B; 
    }
    void kase_insert(CRT _)
    {
        LL res=Exgcd(P,_.P,_.A-A);
        A+=res*P;
        P=P/gcd(P,_.P)*_.P;
        return;
    }
};

LL inv(LL a,LL P){return Pow(a,P-2,P);}
LL C(LL n,LL m,int i){return m>n?0:fac[i][n]*inv(fac[i][m]*fac[i][n-m]%_P[i],_P[i])%_P[i];}
LL Lucas(LL n,LL m,int i){return m?Lucas(n/_P[i],m/_P[i],i)*C(n%_P[i],m%_P[i],i)%_P[i]:1;}
LL rep_C(LL n,LL m,int i){return Lucas(n+m-1,m,i);}

int main()
{
    init();
    int T;
    scanf("%d",&T);
    FOR(Ti,1,T)
    {
        int n,m,k,a[18];
        scanf("%d%d%d",&n,&m,&k);
        FOR(i,1,k)scanf("%d",&a[i]);
        a[k+1]=m;
        CRT sum[2];
        FOR(i,0,1)sum[i].P=_P[i],sum[i].A=1;
        FOR(i,0,1)FOR(j,1,k)(sum[i].A*=rep_C(n,a[j+1]-a[j],i))%=_P[i];
        sum[0].kase_insert(sum[1]);
        printf("Case #%d: %lld\n",Ti,sum[0].A);
    }
    return 0;
}