bzoj 3738: [Ontak2013]Kapitał

.

题意已经够简短就不概括了吧。
题解：
k最大是9，那就直接当所有都$\%10^9$吧。答案再保留k位就好。
把$10^9$拆成$2^9 \times 5^9$，假设我们知道$C_{N+M}^N\%2^9$和$C_{N+M}^N\%5^9$的值，我们就可以用中国剩余定理或扩展欧几里得求出答案。具体来说，就是设$x=C_{N+M}^N,\begin{cases}x \equiv a\pmod{2^9}\\ x \equiv b\pmod{5^9}\end{cases}$，那么$x=k\cdot 2^9\cdot 5^9+a\cdot 5^9\cdot inv(5^9,2^9)+b\cdot 2^9\cdot inv(2^9,5^9)$，模掉$10^9$就只剩后面两项了。

那么怎么求$C_{N+M}^N\%2^9$和$C_{N+M}^N\%5^9$呢？Lucas定理是行不通了，因为模数不是质数。分子分母分别求出来再除（乘逆元）的话，可能会出现这种情况：结果是$\frac{a\cdot p^\alpha}{b\cdot p^\beta}$，当$\alpha$和$\beta$都$\ge9$时上下都=0，而答案可能不等于0。所以我们要保留$\alpha$和$\beta$，而$a$和$b$就可以随便模。求$C$就可以直接用阶乘模质数的幂来搞了。阶乘模质数的幂这样来做：$$
n!\\\
=1\times2\times3\times\cdots\times n\\\
=1\times\cdots\times(p-1)\times(p+1)\times(p+2)\times\cdots\times(2p-1)\times(2p+1)\times\cdots\times n\times p^{\lfloor n/p\rfloor}\times (1\times2\times\cdots\times \lfloor n/p\rfloor)\\\
=1\times\cdots\times(p-1)\times(p+1)\times(p+2)\times\cdots\times(2p-1)\times(2p+1)\times\cdots\times n\times p^{\lfloor n/p\rfloor}\times \lfloor n/p\rfloor!

$$
前面的部分预处理，最后那部分可以递归下去搞，返回一个$a\cdot p^x$。除的时候把x相减，乘上a的逆元就好了。
好像已经做完了。
代码：

#include<bits/stdc++.h>
using namespace std;

long long n,m,k,fac1[513],fac2[1953126];
struct num
{
    long long x,y;//x*p^y
};

num times(num x,num y,long long p)
{
    return (num){x.x*y.x%p,x.y+y.y};
}
long long pow(long long x,long long y,long long p)
{
    if(y==0)
    return 1;
    long long ans=pow(x,y>>1,p);
    ans=ans*ans%p;
    if(y&1)
    ans=ans*x%p;
    return ans;
}
long long exgcd(long long a,long long b,long long&x,long long&y)
{
    if(b==0)
    {
        x=1;
        y=0;
        return a;
    }
    long long tx,ty,d=exgcd(b,a%b,tx,ty);
    x=ty;
    y=tx-(a/b)*ty;
    return d;
}
long long inv(long long a,long long p)//ax%p=1
{
    long long x,y;
    exgcd(a,p,x,y);
    x=(x%p+p)%p;
    return x;
}
num inv(num x,long long p)
{
    x.x=inv(x.x,p);
    x.y*=-1;
    return x;
}
num div(num x,num y,long long p)
{
    return times(x,inv(y,p),p);
}
num wk(long long n,long long p,long long pk,long long* fac)//n!%p^k  pk=p^k
{
    if(n==0)
    return (num){1,0};
    long long hh=pow(fac[pk],n/pk,pk)*fac[n%pk]%pk;
    return times((num){hh,n/p},wk(n/p,p,pk,fac),pk);
}
int main()
{
    scanf("%lld%lld%lld",&n,&m,&k);
    fac1[0]=1;
    for(int i=1;i<=512;i++)
    fac1[i]=fac1[i-1]*((i&1)?i:1)%512;
    fac2[0]=1;
    for(int i=1;i<=1953125;i++)
    fac2[i]=fac2[i-1]*((i%5)?i:1)%1953125;
    num a=wk(n+m,2,512,fac1);
    a=div(a,wk(n,2,512,fac1),512);
    a=div(a,wk(m,2,512,fac1),512);
    num b=wk(n+m,5,1953125,fac2);
    b=div(b,wk(n,5,1953125,fac2),1953125);
    b=div(b,wk(m,5,1953125,fac2),1953125);
    long long hh=min(a.y,b.y);
    long long ans2=a.x*pow(2,a.y-hh,512)%512*pow(inv(5,512),hh,512)%512;
    long long ans5=b.x*pow(5,b.y-hh,1953125)%1953125*pow(inv(2,1953125),hh,1953125)%1953125;
    long long ans=(ans2*1953125%1000000000*inv(1953125,512)%1000000000+ans5*512*inv(512,1953125)%1000000000)%1000000000;
    char s[10];
    sprintf(s,"%%0%lldlld",k);
    printf(s,ans%pow(10,k,1000000001));
}