[ontak2013]Kapitał 解题报告

令$n!=2^a5^bc,2\not| c,5\not |c$，这样就$C_n^m$就便于计算了。a，b很好统计，至于c我们可以先求在$\mod 2^k$和$\mod 5^k$意义下的，然后扩展欧几里得合并即可。

如何求$c \mod 2^k$？
先处理出$[0,2^k-1)$中与2互质的数的阶乘，然后对于n！，与2互质的数就可以直接计算了，然后把所有2的倍数除以2，就是$\lfloor {n \over 2} \rfloor !$，就可以递归计算了。
但是这样处理出来的c还是带5的，所以我们还需要除以$5^b$.

如何扩展欧几里得合并？（也可以直接上中国剩余定理）
考虑方程组

$$\left \{ \begin{aligned} x≡a_0\mod m_0\\ x≡a_1\mod m_1 \end{aligned} \right.$$
那么就有 $x=ym_0+a_0=zm_1+a_1$，就有 $ym_0-zm_1=a_1-a_0$，如果 $(m_0,m_1)|a_1-a_0$，这个方程就有解，就可以化为一个形如 $x≡a\bmod lcm(m_0,m_1)$。
显然合法的方程组和合法的方程是一一对应的，所以我们就可以得到类似中国剩余定理的结论。这玩意儿比中国剩余定理稍微优越一点的地方在于它可以判断是否无解，而且不用要求每个模数两两互质。

有一个hack点是在大多数情况下2的个数是比5多的，最后答案是$2^{a-b}c$，但也有特殊情况下5的个数是比2多的，$5^{b-a}c$，一个很显然的例子就是$C_5^1$。
代码：

#include<cstdio>
#include<iostream>
using namespace std;
#include<cstring>
#include<cmath>
#include<iomanip>
typedef long long LL;
const LL N=1e15+5,M=1e15+5;
const int K=9;

void ex_eu(int a,int b,int &x,int &y){
    if(b){
        ex_eu(b,a%b,y,x);
        y-=a/b*x;
    }
    else x=1,y=0;
}
LL pow(LL prod,int x,int Mod){
    LL ans=1;
    for(;x;x>>=1,prod=prod*prod%Mod)
        if(x&1)
            ans=ans*prod%Mod;
    return ans;
}

const int Max_2=1<<9,Max_5=(int)(1e9)>>9;
const int Log=51;
struct FS{
    LL cnt_2,cnt_5,rest;
};
void out(FS node){
    printf("{cnt_2=%I64d,cnt_5=%I64d,rest=%I64d}\n",node.cnt_2,node.cnt_5,node.rest);
}
LL fac_2[Max_2+5];
LL fac_5[Max_5+5];
FS query(LL n,int k){

    //printf("---query(%I64d,%d)---\n",n,k);

    FS ans=(FS){0};

    LL tmp;
    LL power[Log+5];

    for(tmp=n>>1;tmp;tmp>>=1)ans.cnt_2+=tmp;
    int max_2=1<<k,phi_2=1<<k-1;
    fac_2[0]=1;
    for(int i=1;i<max_2;++i)
        if(i&1)fac_2[i]=fac_2[i-1]*i%max_2;
        else fac_2[i]=fac_2[i-1];
    LL rest_2=1;
    for(tmp=n;tmp;tmp>>=1)rest_2=rest_2*pow(fac_2[max_2-1],tmp/max_2%phi_2,max_2)%max_2*fac_2[tmp%max_2]%max_2;

    //printf("rest_2=%I64d\n",rest_2);

    for(tmp=n/5;tmp;tmp/=5)ans.cnt_5+=tmp;
    int max_5=pow(5,k),phi_5=max_5/5*4;
    fac_5[0]=1;
    for(int i=1;i<max_5;++i)
        if(i%5)fac_5[i]=fac_5[i-1]*i%max_5;
        else fac_5[i]=fac_5[i-1];
    LL rest_5=1;
    for(tmp=n;tmp;tmp/=5)rest_5=rest_5*pow(fac_5[max_5]-1,tmp/max_5%phi_5,max_5)%max_5*fac_5[tmp%max_5]%max_5;

    //printf("rest_5=%I64d\n",rest_5);

    rest_2=rest_2*pow(pow(5,ans.cnt_5%phi_2,max_2),phi_2-1,max_2)%max_2;
    rest_5=rest_5*pow(pow(2,ans.cnt_2%phi_5,max_5),phi_5-1,max_5)%max_5;
    //printf("rest_2=%I64d\n",rest_2);
    //printf("rest_5=%I64d\n",rest_5);
    int x,y;
    ex_eu(max_2,max_5,x,y);
    ans.rest=(x*(rest_5-rest_2)%max_5+max_5)%max_5*max_2+rest_2;

    //printf("ans=");
    //out(ans);

    return ans;
}
int main(){
    freopen("bzoj_3738.in","r",stdin);
    LL n,m;
    int k;
    cin>>n>>m>>k;
    FS a=query(n+m,k),b0=query(n,k),b1=query(m,k);
    a.cnt_2-=b0.cnt_2+b1.cnt_2;
    a.cnt_5-=b0.cnt_5+b1.cnt_5;
    int Mod=pow(10,k),phi=Mod/5<<1;
    a.rest=a.rest*pow(b0.rest,phi-1,Mod)%Mod*pow(b1.rest,phi-1,Mod)%Mod;

    //printf("final ans=");
    //out(a);

    cout<<setfill('0')<<setw(k)<<pow(2,a.cnt_2-a.cnt_5,Mod)*a.rest%Mod<<endl;
}

总结：
①$C_n^m$中，5的个数也有可能比2多！不要想当然。