bzoj 3738: [Ontak2013]Kapitał

题意已经够简短就不概括了吧。
题解：
k最大是9，那就直接当所有都%109吧。答案再保留k位就好。
把109拆成29×59，假设我们知道CNN+M%29和CNN+M%59的值，我们就可以用中国剩余定理或扩展欧几里得求出答案。具体来说，就是设x=CNN+M,{x≡a(mod29) x≡b(mod59)，那么x=k⋅29⋅59+a⋅59⋅inv(59,29)+b⋅29⋅inv(29,59)，模掉109就只剩后面两项了。

那么怎么求CNN+M%29和CNN+M%59呢？Lucas定理是行不通了，因为模数不是质数。分子分母分别求出来再除（乘逆元）的话，可能会出现这种情况：结果是a⋅pαb⋅pβ，当α和β都≥9时上下都=0，而答案可能不等于0。所以我们要保留α和β，而a和b就可以随便模。求C就可以直接用阶乘模质数的幂来搞了。阶乘模质数的幂这样来做：

n! =1×2×3×⋯×n =1×⋯×(p−1)×(p+1)×(p+2)×⋯×(2p−1)×(2p+1)×⋯×n×p⌊n/p⌋×(1×2×⋯×⌊n/p⌋) =1×⋯×(p−1)×(p+1)×(p+2)×⋯×(2p−1)×(2p+1)×⋯×n×p⌊n/p⌋×⌊n/p⌋!

前面的部分预处理，最后那部分可以递归下去搞，返回一个a⋅px。除的时候把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));}